cartodata  5.1.2
volumeutil.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 #ifndef CARTODATA_VOLUME_VOLUMEUTIL_H
35 #define CARTODATA_VOLUME_VOLUMEUTIL_H
36 
37 //--- cartodata --------------------------------------------------------------
38 #include <cartodata/volume/volumeref.h>
39 #include <cartodata/volume/functional.h>
40 //--- cartobase --------------------------------------------------------------
41 #include <cartobase/type/datatypetraits.h>
42 #include <cartobase/smart/rcptr.h>
43 #include <cartobase/containers/nditerator.h>
44 //--- std --------------------------------------------------------------------
45 #include <algorithm>
46 #include <exception>
47 //----------------------------------------------------------------------------
48 
49 namespace carto
50 {
51 
52  //==========================================================================
53  //
54  // DECLARATIONS
55  //
56  //==========================================================================
57 
58  namespace volumeutil {
59 
60  //========================================================================
61  // Generic operators
62  //========================================================================
66 
69  inline bool sameSize( const std::vector<int> & s1,
70  const std::vector<int> & s2 );
72  inline std::vector<int> maxSize( const std::vector<int> & s1,
73  const std::vector<int> & s2 );
75  inline std::vector<int> minSize( const std::vector<int> & s1,
76  const std::vector<int> & s2 );
77 
81  template <typename T, typename UnaryFunction>
82  Volume<typename UnaryFunction::result_type>
83  apply( const Volume<T> & vol, UnaryFunction func );
84  template <typename T, typename U, typename BinaryFunction>
85  Volume<typename BinaryFunction::result_type>
86  apply( const Volume<T> & vol1, const Volume<U> & vol2, BinaryFunction func );
87  template <typename T, typename UnaryFunction>
88  rc_ptr<Volume<typename UnaryFunction::result_type> >
89  apply( const rc_ptr<Volume<T> > & vol, UnaryFunction func );
90  template <typename T, typename U, typename BinaryFunction>
91  rc_ptr<Volume<typename BinaryFunction::result_type> >
92  apply( const rc_ptr<Volume<T> > & vol1, const Volume<U> & vol2, BinaryFunction func );
93 
97  template <typename T, typename UnaryFunction>
98  Volume<T> & selfApply( Volume<T> & vol, UnaryFunction func );
99  template <typename T, typename U, typename BinaryFunction>
100  Volume<T> & selfApply( Volume<T> & vol1, const Volume<U> & vol2, BinaryFunction func );
101  template <typename T, typename UnaryFunction>
102  rc_ptr<Volume<T> > & selfApply( rc_ptr<Volume<T> > & vol, UnaryFunction func );
103  template <typename T, typename U, typename BinaryFunction>
104  rc_ptr<Volume<T> > & selfApply( rc_ptr<Volume<T> > & vol1, const Volume<U> & vol2, BinaryFunction func );
106 
113  template <typename T, typename OUTP, typename UnaryFunction>
114  Volume<OUTP> &
115  applyTowards( const Volume<T> & vol, Volume<OUTP> & dst, UnaryFunction func );
116  template <typename T, typename OUTP, typename UnaryFunction>
117  Volume<OUTP> &
118  applyTowards( const T & cst, Volume<OUTP> & dst, UnaryFunction func );
119  template <typename T, typename U, typename OUTP, typename BinaryFunction>
120  Volume<OUTP> &
121  applyTowards( const Volume<T> & vol1, const Volume<U> & vol2, Volume<OUTP> & dst, BinaryFunction func );
122  template <typename T, typename U, typename OUTP, typename BinaryFunction>
123  Volume<OUTP> &
124  applyTowards( const Volume<T> & vol1, const U & cst2, Volume<OUTP> & dst, BinaryFunction func );
126 
140  template <typename OUTP, typename T, typename BinaryFunction>
141  OUTP accumulate( const Volume<T> & vol, BinaryFunction func, OUTP initial );
142  template <typename OUTP, typename T, typename BinaryFunction>
143  OUTP accumulate( const rc_ptr<Volume<T> > & vol, BinaryFunction func, OUTP initial );
145 
146  } // namespace volumeutil
147 
148  //==========================================================================
149  // Copy functions
150  //==========================================================================
151 
155  template <typename T>
156  void transfer( const Volume<T> & src, Volume<T> & dst );
157  template <typename T>
158  void transfer( const rc_ptr<Volume<T> > & src, rc_ptr<Volume<T> > & dst );
159  template <typename OUTP, typename INP>
160  void transfer( const Volume<INP> & src, Volume<OUTP> & dst );
161  template <typename OUTP, typename INP>
162  void transfer( const rc_ptr<Volume<INP> > & src, rc_ptr<Volume<OUTP> > & dst );
164 
168  template <typename T>
169  Volume<T> deepcopy( const Volume<T> & src,
170  const std::vector<int> & size = std::vector<int>() );
171  template <typename T>
172  rc_ptr<Volume<T> > deepcopy(
173  const rc_ptr<Volume<T> > & src,
174  const std::vector<int> & size = std::vector<int>() );
175  template <typename OUTP, typename INP>
176  Volume<OUTP> deepcopy( const Volume<INP> & src,
177  const std::vector<int> & size = std::vector<int>() );
178  template <typename OUTP, typename INP>
179  rc_ptr<Volume<OUTP> > deepcopy(
180  const rc_ptr<Volume<INP> > & src,
181  const std::vector<int> & size = std::vector<int>() );
183 
188  template <typename T>
189  Volume<T> copy( const Volume<T> & src );
190  template <typename T>
191  rc_ptr<Volume<T> > copy( const rc_ptr<Volume<T> > & src );
192  template <typename OUTP, typename INP>
193  Volume<OUTP> copy( const Volume<INP> & src );
194  template <typename OUTP, typename INP>
195  rc_ptr<Volume<OUTP> > copy( const rc_ptr<Volume<INP> > & src );
197 
200  template <typename T>
201  Volume<T> copyStructure( const Volume<T> & src );
202  template <typename T>
203  rc_ptr<Volume<T> > copyStructure( const rc_ptr<Volume<T> > & src );
204  template <typename OUTP, typename INP>
205  Volume<OUTP> copyStructure( const Volume<INP> & src );
206  template <typename OUTP, typename INP>
207  rc_ptr<Volume<OUTP> > copyStructure( const rc_ptr<Volume<INP> > & src );
209 
210  //==========================================================================
211  // Accumulated values: Min/Max/Sum...
212  //==========================================================================
213 
218  template <typename T>
219  T min( const Volume<T> & vol );
220  template <typename T>
221  T min( const rc_ptr<Volume<T> > & vol );
223 
228  template <typename T>
229  T max( const Volume<T> & vol );
230  template <typename T>
231  T max( const rc_ptr<Volume<T> > & vol );
233 
236  template <typename T>
237  typename DataTypeTraits<T>::LongType sum( const Volume<T> & vol );
238  template <typename OUTP, typename T>
239  OUTP sum( const Volume<T> & vol );
240  template <typename T>
241  typename DataTypeTraits<T>::LongType sum( const rc_ptr<Volume<T> > & vol );
242  template <typename OUTP, typename T>
243  OUTP sum( const rc_ptr<Volume<T> > & vol );
245 
248  template <typename T>
249  bool all( const Volume<T> & vol );
250  template <typename T>
251  bool all( const rc_ptr<Volume<T> > & vol );
253 
256  template <typename T>
257  bool any( const Volume<T> & vol );
258  template <typename T>
259  bool any( const rc_ptr<Volume<T> > & vol );
261 
262  //==========================================================================
263  // Special operators
264  //==========================================================================
265 
273  template <typename T, typename U>
274  Volume<bool> valuesIn( const Volume<T> & volume, const U & set );
275  template <typename T, typename U>
276  rc_ptr<Volume<bool> > valuesIn( const rc_ptr<Volume<T> > & volume, const U & set );
278 
286  template <typename T, typename U>
287  Volume<bool> valuesNotIn( const Volume<T> & volume, const U & set );
288  template <typename T, typename U>
289  rc_ptr<Volume<bool> > valuesNotIn( const rc_ptr<Volume<T> > & volume, const U & set );
291 
297  template <typename T, typename U>
298  void conditionalSet( Volume<T> & volume, const Volume<U> & condition, const T & value );
299  template <typename T, typename U>
300  void conditionalSet( rc_ptr<Volume<T> > & volume, const rc_ptr<Volume<U> > & condition, const T & value );
302 
307  template <typename T>
308  Volume<T> invertMinMax( const Volume<T> & volume );
309  template <typename T>
310  rc_ptr<Volume<T> > invertMinMax( const rc_ptr<Volume<T> > & volume );
312 
314  template <typename T>
315  VolumeRef<T> matrix_product( const Volume<T> & v1, const Volume<T> & v2 );
317  template <typename T>
318  VolumeRef<T> matrix_product( const rc_ptr<Volume<T> > & v1,
319  const rc_ptr<Volume<T> > & v2 );
321  template <typename T>
322  VolumeRef<T> transpose( const Volume<T> & v );
327  template <typename T>
328  VolumeRef<T> transpose( const rc_ptr<Volume<T> > & v, bool copy=false );
329 
330 
332  template <typename T>
333  inline
334  rc_ptr<Volume<T> > diag( const rc_ptr<Volume<T> > & thing )
335  {
336  if( thing->getSizeT() != 1 || thing->getSizeZ() != 1
337  || thing->getSizeY() != 1 )
338  throw std::runtime_error("diag expects a 1D (Nx1x1x1) vector as input");
339 
340  int n = thing->getSizeX();
341  // allocate the new thing
342  rc_ptr<Volume<T> > m( new Volume<T>( n, n ) );
343  // do the operations
344  for( int x = 0; x < n; x++ )
345  m->at( x, x ) = thing->at( x );
346 
347  return m;
348  }
349 
350 
351  template <typename T>
352  inline
353  rc_ptr<Volume<T> > undiag( const rc_ptr<Volume<T> > & thing )
354  {
355  if( thing->getSizeT() != 1 || thing->getSizeZ() != 1 )
356  throw std::runtime_error(
357  "undiag expects a 2D (NxMx1x1) matrix as input");
358 
359  int n = std::min( thing->getSizeX(), thing->getSizeY() );
360  // allocate the new thing
361  rc_ptr<Volume<T> > m( new Volume<T>( n, 1, 1, 1,
362  AllocatorContext::fast() ) );
363  // do the operations
364  for( int x = 0; x < n; x++ )
365  m->at( x, 0 ) = thing->at( x, x );
366 
367  return m;
368  }
369 
370  //==========================================================================
371  // Borders
372  //==========================================================================
373 
382  template <typename T>
383  void setBorders( Volume<T> & volume,
384  const typename Volume<T>::Position4Di & top,
385  const typename Volume<T>::Position4Di & bottom
386  = typename Volume<T>::Position4Di(-1, -1, -1, -1) );
389  template <typename T>
390  void setBorders( rc_ptr<Volume<T> > & volume,
391  const typename Volume<T>::Position4Di & top,
392  const typename Volume<T>::Position4Di & bottom
393  = typename Volume<T>::Position4Di(-1, -1, -1, -1) );
395 
403  template <typename T>
404  void setMinBorders( Volume<T> & volume,
405  const typename Volume<T>::Position4Di & top,
406  const typename Volume<T>::Position4Di & bottom
407  = typename Volume<T>::Position4Di(-1, -1, -1, -1) );
410  template <typename T>
411  void setMinBorders( rc_ptr<Volume<T> > & volume,
412  const typename Volume<T>::Position4Di & top,
413  const typename Volume<T>::Position4Di & bottom
414  = typename Volume<T>::Position4Di(-1, -1, -1, -1) );
416 
417 
418  template<typename T>
419  inline
420  bool isSameVolumeSize( const Volume<T> & v1, const Volume<T> & v2 )
421  {
422  std::vector<int> dim1 = v1.getSize();
423  std::vector<int> dim2 = v2.getSize();
424  int i, n1 = dim1.size(), n2 = dim2.size(), n = std::min( n1, n2 );
425  for( i=0; i<n; ++i )
426  if( dim1[i] != dim2[i] )
427  return false;
428  for( ; i<n1; ++i )
429  if( dim1[i] != 1 )
430  return false;
431  for( ; i<n2; ++i )
432  if( dim2[i] != 1 )
433  return false;
434 
435  return true;
436  }
437 
438 
439  template<typename T>
440  inline
441  bool isSameVolumeSize( const rc_ptr<Volume<T> > & v1,
442  const rc_ptr<Volume<T> > & v2 )
443  {
444  return isSameVolumeSize( *v1, *v2 );
445  }
446 
447 
448  //==========================================================================
449  //
450  // DEFINITIONS
451  //
452  //==========================================================================
453 
454  namespace volumeutil {
455 
456  //========================================================================
457  // Generic operators
458  //========================================================================
459 
460  bool sameSize( const std::vector<int> & s1, const std::vector<int> & s2 )
461  {
462  int s0 = s1.size();
463  int i;
464  if( s2.size() < s0 )
465  s0 = s2.size();
466  for( i=0; i<s0; ++i )
467  if( s1[i] != s2[i] )
468  return false;
469  for( ; i<s1.size(); ++i )
470  if( s1[i] != 1 )
471  return false;
472  for( ; i<s2.size(); ++i )
473  if( s2[i] != 1 )
474  return false;
475  return true;
476  }
477 
478 
479  std::vector<int> maxSize( const std::vector<int> & s1,
480  const std::vector<int> & s2 )
481  {
482  std::vector<int> dim( std::max( s1.size(), s2.size() ), 1 );
483  int s0 = s1.size();
484  int i;
485  if( s2.size() < s0 )
486  s0 = s2.size();
487  for( i=0; i<s0; ++i )
488  dim[i] = std::max( s1[i], s2[i] );
489  for( ; i<s1.size(); ++i )
490  dim[i] = s1[i];
491  for( ; i<s2.size(); ++i )
492  dim[i] = s2[i];
493  return dim;
494  }
495 
496 
497  std::vector<int> minSize( const std::vector<int> & s1,
498  const std::vector<int> & s2 )
499  {
500  std::vector<int> dim( std::min( s1.size(), s2.size() ), 1 );
501  int s0 = s1.size();
502  int i;
503  if( s2.size() < s0 )
504  s0 = s2.size();
505  for( i=0; i<s0; ++i )
506  dim[i] = std::min( s1[i], s2[i] );
507  for( ; i<dim.size(); ++i )
508  dim[i] = 1;
509  return dim;
510  }
511 
512 
513  //--- Volume [op] other --------------------------------------------------
514 
515  template <typename T, typename UnaryFunction>
516  inline
517  Volume<typename UnaryFunction::result_type>
518  apply( const Volume<T> & vol, UnaryFunction func )
519  {
520  typedef typename UnaryFunction::result_type OUTP;
521  Volume<OUTP> output = deepcopy<OUTP,T>(vol);
522  applyTowards( vol, output, func );
523  return output;
524  }
525 
526  template <typename T, typename U, typename BinaryFunction>
527  inline
528  Volume<typename BinaryFunction::result_type>
529  apply( const Volume<T> & vol1, const Volume<U> & vol2, BinaryFunction func )
530  {
531  typedef typename BinaryFunction::result_type OUTP;
532  Volume<OUTP> output = deepcopy<OUTP,T>(
533  vol1, maxSize( vol1.getSize(), vol2.getSize() ) );
534  applyTowards( vol1, vol2, output, func );
535  return output;
536  }
537 
538  //--- VolumeRef [op]= other ----------------------------------------------
539 
540  template <typename T, typename UnaryFunction>
541  inline
542  rc_ptr<Volume<typename UnaryFunction::result_type> >
543  apply( const rc_ptr<Volume<T> > & vol, UnaryFunction func )
544  {
545  typedef typename UnaryFunction::result_type OUTP;
546  rc_ptr<Volume<OUTP> > output = deepcopy<OUTP,T>(vol);
547  applyTowards( *vol, *output, func );
548  return output;
549  }
550 
551  template <typename T, typename U, typename BinaryFunction>
552  inline
553  rc_ptr<Volume<typename BinaryFunction::result_type> >
554  apply( const rc_ptr<Volume<T> > & vol1, const Volume<U> & vol2, BinaryFunction func )
555  {
556  typedef typename BinaryFunction::result_type OUTP;
557  rc_ptr<Volume<OUTP> > output = deepcopy<OUTP,T>(
558  vol1, maxSize( vol1->getSize(), vol2.getSize()) );
559  applyTowards( *vol1, vol2, *output, func );
560  return output;
561  }
562 
563  //--- Volume [op]= other -------------------------------------------------
564 
565  template <typename T, typename UnaryFunction>
566  inline
567  Volume<T> & selfApply( Volume<T> & vol, UnaryFunction func )
568  {
569  applyTowards( vol, vol, func );
570  return vol;
571  }
572 
573  template <typename T, typename U, typename BinaryFunction>
574  inline
575  Volume<T> & selfApply( Volume<T> & vol1, const Volume<U> & vol2, BinaryFunction func )
576  {
577  applyTowards( vol1, vol2, vol1, func );
578  return vol1;
579  }
580 
581  //--- VolumeRef [op]= other ----------------------------------------------
582 
583  template <typename T, typename UnaryFunction>
584  inline
585  rc_ptr<Volume<T> > & selfApply( rc_ptr<Volume<T> > & vol, UnaryFunction func )
586  {
587  selfApply( *vol, func );
588  return vol;
589  }
590 
591  template <typename T, typename U, typename BinaryFunction>
592  inline
593  rc_ptr<Volume<T> > & selfApply( rc_ptr<Volume<T> > & vol1, const Volume<U> & vol2, BinaryFunction func )
594  {
595  selfApply( *vol1, vol2, func );
596  return vol1;
597  }
598 
599  //--- Volume [op] other -> Output ---------------------------------------
600  template <typename T, typename OUTP, typename UnaryFunction>
601  inline
602  Volume<OUTP> &
603  applyTowards( const Volume<T> & vol, Volume<OUTP> & dst,
604  UnaryFunction func )
605  {
606  std::vector<int> size = minSize( vol.getSize(), dst.getSize() );
607  if( !sameSize( size, vol.getSize() )
608  || !sameSize( size, dst.getSize() ) )
609  {
610  // build views
611  carto::rc_ptr<Volume<T> > invvol(
612  const_cast<carto::Volume<T> *>( &vol ) );
613  std::vector<int> pos( size.size(), 0 );
614  carto::VolumeRef<T> invol( invvol, pos, size );
615  carto::rc_ptr<Volume<OUTP> > outvvol( &dst );
616  carto::VolumeRef<OUTP> outvol( outvvol, pos, size );
617  // applyTowards on views of the same size
618  applyTowards( *invol, *outvol, func );
619  // release views and restore refcount on initial volumes
620  invol.reset( 0 );
621  invvol.release();
622  outvol.reset( 0 );
623  outvvol.release();
624  // FIXME TODO: handle parts outside of the intersection
625  return dst;
626  }
627  carto::const_line_NDIterator<T> it( &vol.at( 0 ), vol.getSize(),
628  vol.getStrides() );
629  carto::line_NDIterator<OUTP> dit( &dst.at( 0 ), dst.getSize(),
630  dst.getStrides() );
631  const T *p, *pp;
632  OUTP *d;
633  for( ; !it.ended(); ++it, ++dit )
634  {
635  p = &*it;
636  d = &*dit;
637  for( pp=p + it.line_length(); p!=pp;
638  it.inc_line_ptr( p ), dit.inc_line_ptr( d ) )
639  *d = func( *p );
640  }
641  return dst;
642  }
643 
644  template <typename T, typename OUTP, typename UnaryFunction>
645  inline
646  Volume<OUTP> &
647  applyTowards( const T & cst, Volume<OUTP> & dst,
648  UnaryFunction func )
649  {
650  carto::line_NDIterator<OUTP> dit( &dst.at( 0 ), dst.getSize(),
651  dst.getStrides() );
652  OUTP *d, *dd;
653  OUTP value = func( cst );
654  for( ; !dit.ended(); ++dit )
655  {
656  d = &*dit;
657  for( dd=d + dit.line_length(); d!=dd; dit.inc_line_ptr( d ) )
658  *d = value;
659  }
660  return dst;
661  }
662 
663  template <typename T, typename U, typename OUTP, typename BinaryFunction>
664  inline
665  Volume<OUTP> &
666  applyTowards( const Volume<T> & vol1, const Volume<U> & vol2,
667  Volume<OUTP> & dst, BinaryFunction func )
668  {
669  std::vector<int> sizei = minSize( vol1.getSize(), vol2.getSize() );
670  std::vector<int> size = minSize( sizei, dst.getSize() );
671  if( !sameSize( size, vol1.getSize() )
672  || !sameSize( size, vol2.getSize() )
673  || !sameSize( size, dst.getSize() ) )
674  {
675  // build views
676  carto::rc_ptr<Volume<T> > invvol1(
677  const_cast<carto::Volume<T> *>( &vol1 ) );
678  std::vector<int> pos( size.size(), 0 );
679  carto::VolumeRef<T> invol1( invvol1, pos, size );
680  carto::rc_ptr<Volume<U> > invvol2(
681  const_cast<carto::Volume<U> *>( &vol2 ) );
682  carto::VolumeRef<U> invol2( invvol2, pos, size );
683  carto::rc_ptr<Volume<OUTP> > outvvol( &dst );
684  carto::VolumeRef<OUTP> outvol( outvvol, pos, size );
685  // applyTowards on views of the same size
686  applyTowards( *invol1, *invol2, *outvol, func );
687  // release views and restore refcount on initial volumes
688  invol1.reset( 0 );
689  invvol1.release();
690  invol2.reset( 0 );
691  invvol2.release();
692  outvol.reset( 0 );
693  outvvol.release();
694  // FIXME TODO: handle parts outside of the intersection
695  return dst;
696  }
697 
698  carto::const_line_NDIterator<T> it1( &vol1.at( 0 ), vol1.getSize(),
699  vol1.getStrides() );
700  carto::const_line_NDIterator<U> it2( &vol2.at( 0 ), vol2.getSize(),
701  vol2.getStrides() );
702  carto::line_NDIterator<OUTP> dit( &dst.at( 0 ), dst.getSize(),
703  dst.getStrides() );
704  const T *p, *pp;
705  const U *p2;
706  OUTP *d;
707  for( ; !it1.ended(); ++it1, ++it2, ++dit )
708  {
709  p = &*it1;
710  p2 = &*it2;
711  d = &*dit;
712  for( pp=p + it1.line_length(); p!=pp;
713  it1.inc_line_ptr( p ), it2.inc_line_ptr( p2 ),
714  dit.inc_line_ptr( d ) )
715  *d = func( *p, *p2 );
716  }
717  return dst;
718  }
719 
720  template <typename T, typename U, typename OUTP, typename BinaryFunction>
721  inline
722  Volume<OUTP> &
723  applyTowards( const T & cst1, const Volume<U> & vol2,
724  Volume<OUTP> & dst, BinaryFunction func )
725  {
726  std::vector<int> size = minSize( vol2.getSize(), dst.getSize() );
727  if( !sameSize( size, vol2.getSize() )
728  || !sameSize( size, dst.getSize() ) )
729  {
730  // build views
731  std::vector<int> pos( size.size(), 0 );
732  carto::rc_ptr<Volume<U> > invvol2(
733  const_cast<carto::Volume<U> *>( &vol2 ) );
734  carto::VolumeRef<U> invol2( invvol2, pos, size );
735  carto::rc_ptr<Volume<OUTP> > outvvol( &dst );
736  carto::VolumeRef<OUTP> outvol( outvvol, pos, size );
737  // applyTowards on views of the same size
738  applyTowards( cst1, *invol2, *outvol, func );
739  // release views and restore refcount on initial volumes
740  invol2.reset( 0 );
741  invvol2.release();
742  outvol.reset( 0 );
743  outvvol.release();
744  // FIXME TODO: handle parts outside of the intersection
745  return dst;
746  }
747 
748  carto::const_line_NDIterator<U> it2( &vol2.at( 0 ), vol2.getSize(),
749  vol2.getStrides() );
750  carto::line_NDIterator<OUTP> dit( &dst.at( 0 ), dst.getSize(),
751  dst.getStrides() );
752  const U *p2, *pp;
753  OUTP *d;
754  for( ; !it2.ended(); ++it2, ++dit )
755  {
756  p2 = &*it2;
757  d = &*dit;
758  for( pp=p2+ it2.line_length(); p2!=pp;
759  it2.inc_line_ptr( p2 ), dit.inc_line_ptr( d ) )
760  *d = func( cst1, *p2 );
761  }
762  return dst;
763  }
764 
765  template <typename T, typename U, typename OUTP, typename BinaryFunction>
766  inline
767  Volume<OUTP> &
768  applyTowards( const Volume<T> & vol1, const U & cst2,
769  Volume<OUTP> & dst, BinaryFunction func )
770  {
771  std::vector<int> size = minSize( vol1.getSize(), dst.getSize() );
772  if( !sameSize( size, vol1.getSize() )
773  || !sameSize( size, dst.getSize() ) )
774  {
775  // build views
776  carto::rc_ptr<Volume<T> > invvol1(
777  const_cast<carto::Volume<T> *>( &vol1 ) );
778  std::vector<int> pos( size.size(), 0 );
779  carto::VolumeRef<T> invol1( invvol1, pos, size );
780  carto::rc_ptr<Volume<OUTP> > outvvol( &dst );
781  carto::VolumeRef<OUTP> outvol( outvvol, pos, size );
782  // applyTowards on views of the same size
783  applyTowards( *invol1, cst2, *outvol, func );
784  // release views and restore refcount on initial volumes
785  invol1.reset( 0 );
786  invvol1.release();
787  outvol.reset( 0 );
788  outvvol.release();
789  // FIXME TODO: handle parts outside of the intersection
790  return dst;
791  }
792 
793  carto::const_line_NDIterator<T> it1( &vol1.at( 0 ), vol1.getSize(),
794  vol1.getStrides() );
795  carto::line_NDIterator<OUTP> dit( &dst.at( 0 ), dst.getSize(),
796  dst.getStrides() );
797  const T *p, *pp;
798  OUTP *d;
799  for( ; !it1.ended(); ++it1, ++dit )
800  {
801  p = &*it1;
802  d = &*dit;
803  for( pp=p + it1.line_length(); p!=pp;
804  it1.inc_line_ptr( p ), dit.inc_line_ptr( d ) )
805  *d = func( *p, cst2 );
806  }
807  return dst;
808  }
809 
810  //--- Volume: accumulate -------------------------------------------------
811 
812  template <typename OUTP, typename T, typename BinaryFunction>
813  inline
814  OUTP accumulate( const Volume<T> & vol, BinaryFunction func, OUTP initial )
815  {
816  carto::const_line_NDIterator<T> it( &vol.at( 0 ), vol.getSize(),
817  vol.getStrides() );
818  const T *p, *pp;
819  for( ; !it.ended(); ++it )
820  {
821  p = &*it;
822  for( pp=p + it.line_length(); p!=pp; it.inc_line_ptr( p ) )
823  initial = func( initial, *p );
824  }
825  return initial;
826  }
827 
828  //--- VolumeRef: accumulate ----------------------------------------------
829 
830  template <typename OUTP, typename T, typename BinaryFunction>
831  inline
832  OUTP accumulate( const rc_ptr<Volume<T> > & vol, BinaryFunction func, OUTP initial )
833  {
834  return accumulate<OUTP, T, BinaryFunction>( *vol, func, initial );
835  }
836 
837  } // namespace volumeutil
838 
839  //==========================================================================
840  // Copy functions
841  //==========================================================================
842 
843  //--- transfer -------------------------------------------------------------
844 
845  template <typename T>
846  inline
847  void transfer( const Volume<T> & src, Volume<T> & dst )
848  {
849  transfer<T,T>( src, dst );
850  }
851 
852  template <typename T>
853  inline
854  void transfer( const rc_ptr<Volume<T> > & src, rc_ptr<Volume<T> > & dst )
855  {
856  transfer<T,T>( src, dst );
857  }
858 
859  template <typename OUTP, typename INP>
860  inline
861  void transfer( const Volume<INP> & src, Volume<OUTP> & dst )
862  {
863  volumeutil::applyTowards( src, dst, volumeutil::identity<INP>() );
864  }
865 
866  template <typename OUTP, typename INP>
867  inline
868  void transfer( const rc_ptr<Volume<INP> > & src, rc_ptr<Volume<OUTP> > & dst )
869  {
870  transfer( *src, *dst );
871  }
872 
873  //--- deepcopy -------------------------------------------------------------
874 
875  template <typename T>
876  inline
877  Volume<T> deepcopy( const Volume<T> & src, const std::vector<int> & size )
878  {
879  return deepcopy<T,T>( src, size );
880  }
881 
882  template <typename T>
883  inline
884  rc_ptr<Volume<T> > deepcopy( const rc_ptr<Volume<T> > & src,
885  const std::vector<int> & size )
886  {
887  return deepcopy<T,T>( src, size );
888  }
889 
890  template <typename OUTP, typename INP>
891  inline
892  Volume<OUTP> deepcopy( const Volume<INP> & src,
893  const std::vector<int> & size )
894  {
895  std::vector<int> nsize = size;
896  if( nsize.empty() )
897  nsize = src.getSize();
898  Volume<OUTP> dst( nsize,
899  src.allocatorContext(),
900  src.allocatorContext().isAllocated() );
901  dst.copyHeaderFrom( src.header() );
902  if( src.allocatorContext().isAllocated()
903  && src.allocatorContext().allocatorType()
904  != AllocatorStrategy::Unallocated )
905  transfer( src, dst );
906 
907  if( src.refVolume() )
908  {
909  rc_ptr<Volume<INP> > srcparent = src.refVolume();
910  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
911  srcparent->getSize(),
912  srcparent->allocatorContext(),
913  srcparent->allocatorContext().isAllocated() ) );
914  if( srcparent->allocatorContext().isAllocated()
915  && srcparent->allocatorContext().allocatorType()
916  != AllocatorStrategy::Unallocated )
917  transfer( srcparent, dstparent );
918  dst.setRefVolume( dstparent );
919  dst.setPosInRefVolume( src.posInRefVolume() );
920  }
921 
922  rc_ptr<Volume<INP> > srcchild = src.refVolume();
923  rc_ptr<Volume<OUTP> > dstchild = dst.refVolume();
924 
925  while( srcchild->refVolume().get() )
926  {
927  rc_ptr<Volume<INP> > srcparent = srcchild->refVolume();
928  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
929  srcparent->getSize(),
930  srcparent->allocatorContext(),
931  srcparent->allocatorContext().isAllocated() ) );
932  dstparent->copyHeaderFrom( srcparent->header() );
933  if( srcparent->allocatorContext().isAllocated()
934  && srcparent->allocatorContext().allocatorType()
935  != AllocatorStrategy::Unallocated )
936  transfer( srcparent, dstparent );
937  dstchild->setRefVolume( dstparent );
938  dstchild->setPosInRefVolume( srcchild->posInRefVolume() );
939  srcchild = srcparent;
940  dstchild = dstparent;
941  }
942 
943  return dst;
944  }
945 
946  template <typename OUTP, typename INP>
947  inline
948  rc_ptr<Volume<OUTP> > deepcopy( const rc_ptr<Volume<INP> > & src,
949  const std::vector<int> & size )
950  {
951  std::vector<int> nsize = size;
952  if( nsize.empty() )
953  nsize = src->getSize();
954  rc_ptr<Volume<OUTP> > dst( new Volume<OUTP>(
955  nsize,
956  src->allocatorContext(),
957  src->allocatorContext().isAllocated() ) );
958  dst->copyHeaderFrom( src->header() );
959  if( src->allocatorContext().isAllocated()
960  && src->allocatorContext().allocatorType()
961  != AllocatorStrategy::Unallocated )
962  transfer( src, dst );
963 
964  rc_ptr<Volume<INP> > srcchild = src;
965  rc_ptr<Volume<OUTP> > dstchild = dst;
966 
967  while( srcchild->refVolume().get() )
968  {
969  rc_ptr<Volume<INP> > srcparent = srcchild->refVolume();
970  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
971  srcparent->getSize(),
972  srcparent->allocatorContext(),
973  srcparent->allocatorContext().isAllocated() ) );
974  dstparent->copyHeaderFrom( srcparent->header() );
975  if( srcparent->allocatorContext().isAllocated()
976  && srcparent->allocatorContext().allocatorType()
977  != AllocatorStrategy::Unallocated )
978  transfer( srcparent, dstparent );
979  dstchild->setRefVolume( dstparent );
980  dstchild->setPosInRefVolume( srcchild->posInRefVolume() );
981  srcchild = srcparent;
982  dstchild = dstparent;
983  }
984 
985  return dst;
986  }
987 
988 
989  //--- copy -----------------------------------------------------------------
990 
991  template <typename T>
992  inline
993  Volume<T> copy( const Volume<T> & src )
994  {
995  return copy<T,T>( src );
996  }
997 
998  template <typename T>
999  inline
1000  rc_ptr<Volume<T> > copy( const rc_ptr<Volume<T> > & src )
1001  {
1002  return copy<T,T>( src );
1003  }
1004 
1005  template <typename OUTP, typename INP>
1006  inline
1007  Volume<OUTP> copy( const Volume<INP> & src )
1008  {
1009  Volume<OUTP> dst( src.getSize() );
1010  dst.copyHeaderFrom( src.header() );
1011  transfer( src, dst );
1012  return dst;
1013  }
1014 
1015  template <typename OUTP, typename INP>
1016  inline
1017  rc_ptr<Volume<OUTP> > copy( const rc_ptr<Volume<INP> > & src )
1018  {
1019  rc_ptr<Volume<OUTP> > dst( new Volume<OUTP>( src->getSize() ) );
1020  dst->copyHeaderFrom( src->header() );
1021  transfer( src, dst );
1022  return dst;
1023  }
1024 
1025 
1026  //--- copyStructure --------------------------------------------------------
1027 
1028  template <typename T>
1029  inline
1030  Volume<T> copyStructure( const Volume<T> & src )
1031  {
1032  return copyStructure<T,T>( src );
1033  }
1034 
1035  template <typename T>
1036  inline
1037  rc_ptr<Volume<T> > copyStructure( const rc_ptr<Volume<T> > & src )
1038  {
1039  return copyStructure<T,T>( src );
1040  }
1041 
1042  template <typename OUTP, typename INP>
1043  inline
1044  Volume<OUTP> copyStructure( const Volume<INP> & src )
1045  {
1046  Volume<OUTP> dst( src.getSize(),
1047  src.allocatorContext(),
1048  src.allocatorContext().isAllocated() );
1049  dst.copyHeaderFrom( src.header() );
1050  // if( src.allocatorContext().isAllocated() )
1051  // dst.fill(0);
1052 
1053  if( src.refVolume() )
1054  {
1055  rc_ptr<Volume<INP> > srcparent = src.refVolume();
1056  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
1057  srcparent->getSize(),
1058  srcparent->allocatorContext(),
1059  srcparent->allocatorContext().isAllocated() ) );
1060  // if( srcparent->allocatorContext().isAllocated() )
1061  // dstparent->fill(0);
1062  dst.setRefVolume( dstparent );
1063  dst.setPosInRefVolume( src.posInRefVolume() );
1064 
1065  rc_ptr<Volume<INP> > srcchild = src.refVolume();
1066  rc_ptr<Volume<OUTP> > dstchild = dst.refVolume();
1067 
1068  while( srcchild->refVolume().get() )
1069  {
1070  rc_ptr<Volume<INP> > srcparent = srcchild->refVolume();
1071  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
1072  srcparent->getSize(),
1073  srcparent->allocatorContext(),
1074  srcparent->allocatorContext().isAllocated() ) );
1075  dstparent->copyHeaderFrom( srcparent->header() );
1076  // if( srcparent->allocatorContext().isAllocated() )
1077  // dstparent->fill(0);
1078  dstchild->setRefVolume( dstparent );
1079  dstchild->setPosInRefVolume( srcchild->posInRefVolume() );
1080  srcchild = srcparent;
1081  dstchild = dstparent;
1082  }
1083 
1084  }
1085 
1086  return dst;
1087  }
1088 
1089  template <typename OUTP, typename INP>
1090  inline
1091  rc_ptr<Volume<OUTP> > copyStructure( const rc_ptr<Volume<INP> > & src )
1092  {
1093  if( !src.get() )
1094  return rc_ptr<Volume<OUTP> >( (Volume<OUTP>*)0 );
1095 
1096  rc_ptr<Volume<OUTP> > dst( new Volume<OUTP>(
1097  src->getSize(),
1098  src->allocatorContext(),
1099  src->allocatorContext().isAllocated() ) );
1100  dst->copyHeaderFrom( src->header() );
1101  // if( src->allocatorContext().isAllocated() )
1102  // dst->fill(0);
1103 
1104  rc_ptr<Volume<INP> > srcchild = src;
1105  rc_ptr<Volume<OUTP> > dstchild = dst;
1106 
1107  while( srcchild->refVolume().get() )
1108  {
1109  rc_ptr<Volume<INP> > srcparent = srcchild->refVolume();
1110  rc_ptr<Volume<OUTP> > dstparent( new Volume<OUTP>(
1111  srcparent->getSize(),
1112  srcparent->allocatorContext(),
1113  srcparent->allocatorContext().isAllocated() ) );
1114  dstparent->copyHeaderFrom( srcparent->header() );
1115  // if( srcparent->allocatorContext().isAllocated() )
1116  // dstparent->fill(0);
1117  dstchild->setRefVolume( dstparent );
1118  dstchild->setPosInRefVolume( srcchild->posInRefVolume() );
1119  srcchild = srcparent;
1120  dstchild = dstparent;
1121  }
1122 
1123  return dst;
1124  }
1125 
1126 
1127  //==========================================================================
1128  // Min/Max/Sum/Any/All
1129  //==========================================================================
1130 
1131  template <typename T>
1132  inline
1133  T min( const Volume<T> & vol )
1134  {
1135  if( vol.getSizeX() == 0 && vol.getSizeY() == 0 &&
1136  vol.getSizeZ() == 0 && vol.getSizeT() == 0 )
1137  throw std::runtime_error("Cannot compute min of an empty volume");
1138  return accumulate( vol, volumeutil::select_min<T>(), vol.at(0, 0, 0, 0) );
1139  }
1140 
1141  template <typename T>
1142  inline
1143  T min( const rc_ptr<Volume<T> > & vol )
1144  {
1145  if( !vol.get() || ( vol->getSizeX() == 0 && vol->getSizeY() == 0 &&
1146  vol->getSizeZ() == 0 && vol->getSizeT() == 0 ) )
1147  throw std::runtime_error("Cannot compute min of an empty volume");
1148  return accumulate( vol, volumeutil::select_min<T>(), vol->at(0, 0, 0, 0) );
1149  }
1150 
1151  template <typename T>
1152  inline
1153  T max( const Volume<T> & vol )
1154  {
1155  if( vol.getSizeX() == 0 && vol.getSizeY() == 0 &&
1156  vol.getSizeZ() == 0 && vol.getSizeT() == 0 )
1157  throw std::runtime_error("Cannot compute max of an empty volume");
1158  return accumulate( vol, volumeutil::select_max<T>(), vol.at(0, 0, 0, 0) );
1159  }
1160 
1161  template <typename T>
1162  inline
1163  T max( const rc_ptr<Volume<T> > & vol )
1164  {
1165  if( !vol.get() || ( vol->getSizeX() == 0 && vol->getSizeY() == 0 &&
1166  vol->getSizeZ() == 0 && vol->getSizeT() == 0 ) )
1167  throw std::runtime_error("Cannot compute max of an empty volume");
1168  return accumulate( vol, volumeutil::select_max<T>(), vol->at(0, 0, 0, 0) );
1169  }
1170 
1171  template <typename T>
1172  inline
1173  typename DataTypeTraits<T>::LongType sum( const Volume<T> & vol )
1174  {
1175  return sum<typename DataTypeTraits<T>::LongType,T>( vol );
1176  }
1177 
1178 
1179  template <typename OUTP, typename T>
1180  inline
1181  OUTP sum( const Volume<T> & vol )
1182  {
1183  return accumulate( vol, volumeutil::plus<OUTP,T>(), static_cast<OUTP>(0) );
1184  }
1185 
1186  template <typename T>
1187  inline
1188  typename DataTypeTraits<T>::LongType sum( const rc_ptr<Volume<T> > & vol )
1189  {
1190  return sum<typename DataTypeTraits<T>::LongType,T>( vol );
1191  }
1192 
1193  template <typename OUTP, typename T>
1194  inline
1195  OUTP sum( const rc_ptr<Volume<T> > & vol )
1196  {
1197  if( !vol.get() )
1198  return static_cast<OUTP>(0);
1199  return accumulate( vol, volumeutil::plus<OUTP,T>(), static_cast<OUTP>(0) );
1200  }
1201 
1202  template <typename T>
1203  inline
1204  bool all( const Volume<T> & vol )
1205  {
1206  return volumeutil::accumulate( vol, volumeutil::logical_and<bool,T>(), true );
1207  }
1208 
1209  template <typename T>
1210  inline
1211  bool all( const rc_ptr<Volume<T> > & vol )
1212  {
1213  return volumeutil::accumulate( *vol, volumeutil::logical_and<bool,T>(), true );
1214  }
1215 
1216  template <typename T>
1217  inline
1218  bool any( const Volume<T> & vol )
1219  {
1220  return volumeutil::accumulate( vol, volumeutil::logical_or<bool,T>(), false );
1221  }
1222 
1223  template <typename T>
1224  inline
1225  bool any( const rc_ptr<Volume<T> > & vol )
1226  {
1227  return volumeutil::accumulate( *vol, volumeutil::logical_or<bool,T>(), false );
1228  }
1229 
1230  //==========================================================================
1231  // Special operators
1232  //==========================================================================
1233 
1234  namespace internal {
1235  template <typename T, typename U>
1236  struct inSet: public std::binary_function<T, U, bool>
1237  {
1238  bool operator() ( const T & x, const U & y )
1239  {
1240  for( typename U::const_iterator k = y.begin(); k != y.end(); ++k )
1241  if( x == *k )
1242  return true;
1243  return false;
1244  }
1245  };
1246  } // namespace internal
1247 
1248  template <typename T, typename U>
1249  inline
1250  Volume<bool> valuesIn( const Volume<T> & volume, const U & set )
1251  {
1252  Volume<bool> output = copyStructure<bool, T>( volume );
1253  volumeutil::applyTowards( volume, output, std::bind2nd( internal::inSet<T,U>(), set ) );
1254  return output;
1255  }
1256 
1257  template <typename T, typename U>
1258  inline
1259  rc_ptr<Volume<bool> > valuesIn( const rc_ptr<Volume<T> > & volume, const U & set )
1260  {
1261  rc_ptr<Volume<bool> > output = copyStructure<bool, T>( volume );
1262  volumeutil::applyTowards( *volume, *output, std::bind2nd( internal::inSet<T,U>(), set ) );
1263  return output;
1264  }
1265 
1266  namespace internal {
1267  template <typename T, typename U>
1268  struct notInSet: public std::binary_function<T, U, bool>
1269  {
1270  bool operator() ( const T & x, const U & y )
1271  {
1272  for( typename U::const_iterator k = y.begin(); k != y.end(); ++k )
1273  if( x == *k )
1274  return false;
1275  return true;
1276  }
1277  };
1278  } // namespace internal
1279 
1280  template <typename T, typename U>
1281  inline
1282  Volume<bool> valuesNotIn( const Volume<T> & volume, const U & set )
1283  {
1284  Volume<bool> output = copyStructure<bool, T>( volume );
1285  volumeutil::applyTowards( volume, output, std::bind2nd( internal::notInSet<T,U>(), set ) );
1286  return output;
1287  }
1288 
1289  template <typename T, typename U>
1290  inline
1291  rc_ptr<Volume<bool> > valuesNotIn( const rc_ptr<Volume<T> > & volume, const U & set )
1292  {
1293  rc_ptr<Volume<bool> > output = copyStructure<bool, T>( volume );
1294  volumeutil::applyTowards( *volume, *output, std::bind2nd( internal::notInSet<T,U>(), set ) );
1295  return output;
1296  }
1297 
1298  namespace internal {
1299  template <typename T, typename U>
1300  struct changeIf: public std::binary_function<T, U, T>
1301  {
1302  changeIf( const T & value ): _value(value) {}
1303  bool operator() ( const T & x, const U & y )
1304  {
1305  return y ? _value : x;
1306  }
1307  private:
1308  T _value;
1309  };
1310  } // namespace internal
1311 
1312  template <typename T, typename U>
1313  inline
1314  void conditionalSet( Volume<T> & volume, const Volume<U> & condition, const T & value )
1315  {
1316  volumeutil::selfApply( volume, condition, internal::changeIf<T,U>(value) );
1317  }
1318 
1319  template <typename T, typename U>
1320  inline
1321  void conditionalSet( rc_ptr<Volume<T> > & volume, const rc_ptr<Volume<U> > & condition, const T & value )
1322  {
1323  volumeutil::selfApply( volume, *condition, internal::changeIf<T,U>(value) );
1324  }
1325 
1326  namespace internal {
1327  template <typename T>
1328  struct invMinMax: public std::unary_function<T, T>
1329  {
1330  invMinMax( const T & min, const T & max ):
1331  _min(min),
1332  _max(max)
1333  {}
1334 
1335  T operator() ( const T & x )
1336  {
1337  return max - x + min;
1338  }
1339  private:
1340  T _min;
1341  T _max;
1342  };
1343  } // namespace internal
1344 
1345  template <typename T>
1346  inline
1347  Volume<T> invertMinMax( const Volume<T> & volume )
1348  {
1349  T min = min( volume );
1350  T max = max( volume );
1351  volumeutil::apply( volume, internal::invMinMax<T>(min, max) );
1352  }
1353 
1354  template <typename T>
1355  inline
1356  rc_ptr<Volume<T> > invertMinMax( const rc_ptr<Volume<T> > & volume )
1357  {
1358  T min = min( volume );
1359  T max = max( volume );
1360  volumeutil::apply( volume, internal::invMinMax<T>(min, max) );
1361  }
1362 
1363  //==========================================================================
1364  // Borders
1365  //==========================================================================
1366 
1367  template <typename T>
1368  inline
1369  void setBorders( Volume<T> & volume,
1370  const typename Volume<T>::Position4Di & top,
1371  const typename Volume<T>::Position4Di & bottom )
1372  {
1373  const typename Volume<T>::Position4Di & bot =
1374  ( bottom == typename Volume<T>::Position4Di(-1, -1, -1, -1) ?
1375  top : bottom );
1376  std::vector<int> b = volume.getBorders();
1377  if( b[0] != top[0] || b[1] != bot[0] ||
1378  b[2] != top[1] || b[3] != bot[1] ||
1379  b[4] != top[2] || b[5] != bot[2] ||
1380  b[6] != top[3] || b[7] != bot[3] )
1381  {
1382  rc_ptr<Volume<T> > parent( new Volume<T>(
1383  volume.getSizeX() + top[0] + bot[0],
1384  volume.getSizeY() + top[1] + bot[1],
1385  volume.getSizeZ() + top[2] + bot[2],
1386  volume.getSizeT() + top[3] + bot[3] ) );
1387  parent->copyHeaderFrom( volume.header() );
1388  typename Volume<T>::Position4Di size( volume.getSizeX(),
1389  volume.getSizeY(),
1390  volume.getSizeZ(),
1391  volume.getSizeT() );
1392  Volume<T> view( parent, top, size );
1393  transfer( volume, view );
1394  volume.reallocate( volume.getSizeX(), volume.getSizeY(),
1395  volume.getSizeZ(), volume.getSizeT(),
1396  true, volume.allocatorContext(), false );
1397  volume.setRefVolume( parent );
1398  volume.setPosInRefVolume( top );
1399  }
1400  }
1401 
1402  template <typename T>
1403  inline
1404  void setBorders( rc_ptr<Volume<T> > & volume,
1405  const typename Volume<T>::Position4Di & top,
1406  const typename Volume<T>::Position4Di & bottom )
1407  {
1408  const typename Volume<T>::Position4Di & bot =
1409  ( bottom == typename Volume<T>::Position4Di(-1, -1, -1, -1) ?
1410  top : bottom );
1411  std::vector<int> b = volume->getBorders();
1412  if( b[0] != top[0] || b[1] != bot[0] ||
1413  b[2] != top[1] || b[3] != bot[1] ||
1414  b[4] != top[2] || b[5] != bot[2] ||
1415  b[6] != top[3] || b[7] != bot[3] )
1416  {
1417  rc_ptr<Volume<T> > parent( new Volume<T>(
1418  volume->getSizeX() + top[0] + bot[0],
1419  volume->getSizeY() + top[1] + bot[1],
1420  volume->getSizeZ() + top[2] + bot[2],
1421  volume->getSizeT() + top[3] + bot[3] ) );
1422  parent->copyHeaderFrom( volume->header() );
1423  typename Volume<T>::Position4Di size( volume->getSizeX(),
1424  volume->getSizeY(),
1425  volume->getSizeZ(),
1426  volume->getSizeT() );
1427  rc_ptr<Volume<T> > view( new Volume<T>( parent, top, size ) );
1428  transfer( volume, view );
1429  volume->reallocate( volume->getSizeX(), volume->getSizeY(),
1430  volume->getSizeZ(), volume->getSizeT(),
1431  true, volume->allocatorContext(), false );
1432  volume->setRefVolume( parent );
1433  volume->setPosInRefVolume( top );
1434  }
1435  }
1436 
1437  template <typename T>
1438  inline
1439  void setMinBorders( Volume<T> & volume,
1440  const typename Volume<T>::Position4Di & top,
1441  const typename Volume<T>::Position4Di & bottom )
1442  {
1443  const typename Volume<T>::Position4Di & bot =
1444  ( bottom == typename Volume<T>::Position4Di(-1, -1, -1, -1) ?
1445  top : bottom );
1446  std::vector<int> b = volume.getBorders();
1447  if( b[0] < top[0] || b[1] < bot[0] ||
1448  b[2] < top[1] || b[3] < bot[1] ||
1449  b[4] < top[2] || b[5] < bot[2] ||
1450  b[6] < top[3] || b[7] < bot[3] )
1451  {
1452  rc_ptr<Volume<T> > parent( new Volume<T>(
1453  volume.getSizeX() + top[0] + bot[0],
1454  volume.getSizeY() + top[1] + bot[1],
1455  volume.getSizeZ() + top[2] + bot[2],
1456  volume.getSizeT() + top[3] + bot[3] ) );
1457  parent->copyHeaderFrom( volume.header() );
1458  typename Volume<T>::Position4Di size( volume.getSizeX(),
1459  volume.getSizeY(),
1460  volume.getSizeZ(),
1461  volume.getSizeT() );
1462  Volume<T> view( parent, top, size );
1463  transfer( volume, view );
1464  volume.reallocate( volume.getSizeX(), volume.getSizeY(),
1465  volume.getSizeZ(), volume.getSizeT(),
1466  true, volume.allocatorContext(), false );
1467  volume.setRefVolume( parent );
1468  volume.setPosInRefVolume( top );
1469  }
1470  }
1471 
1472  template <typename T>
1473  inline
1474  void setMinBorders( rc_ptr<Volume<T> > & volume,
1475  const typename Volume<T>::Position4Di & top,
1476  const typename Volume<T>::Position4Di & bottom )
1477  {
1478  const typename Volume<T>::Position4Di & bot =
1479  ( bottom == typename Volume<T>::Position4Di(-1, -1, -1, -1) ?
1480  top : bottom );
1481  std::vector<int> b = volume->getBorders();
1482  if( b[0] < top[0] || b[1] < bot[0] ||
1483  b[2] < top[1] || b[3] < bot[1] ||
1484  b[4] < top[2] || b[5] < bot[2] ||
1485  b[6] < top[3] || b[7] < bot[3] )
1486  {
1487  rc_ptr<Volume<T> > parent( new Volume<T>(
1488  volume->getSizeX() + top[0] + bot[0],
1489  volume->getSizeY() + top[1] + bot[1],
1490  volume->getSizeZ() + top[2] + bot[2],
1491  volume->getSizeT() + top[3] + bot[3] ) );
1492  parent->copyHeaderFrom( volume->header() );
1493  typename Volume<T>::Position4Di size( volume->getSizeX(),
1494  volume->getSizeY(),
1495  volume->getSizeZ(),
1496  volume->getSizeT() );
1497  rc_ptr<Volume<T> > view( new Volume<T>( parent, top, size ) );
1498  transfer( volume, view );
1499  volume->reallocate( volume->getSizeX(), volume->getSizeY(),
1500  volume->getSizeZ(), volume->getSizeT(),
1501  true, volume->allocatorContext(), false );
1502  volume->setRefVolume( parent );
1503  volume->setPosInRefVolume( top );
1504  }
1505  }
1506 
1507  //==========================================================================
1508  //
1509  // OLD CLASSES / METHODS
1510  //
1511  //==========================================================================
1512  template <typename T> class VolumeRef;
1513 
1514 
1520  template <typename T, class BinaryFunction>
1521  class UnaryFromConstantBinaryFunctor
1522  {
1523  public:
1524  inline UnaryFromConstantBinaryFunctor( const T & x, BinaryFunction func )
1525  : value( x ), f( func ) {}
1526  inline T operator () ( const T & y ) const
1527  { return f( y, value ); }
1528  T value;
1529  BinaryFunction f;
1530  };
1531 
1532 
1537  template <typename T, class BinaryFunction>
1538  class UnaryFromConstantBinaryFunctor2
1539  {
1540  public:
1541  inline UnaryFromConstantBinaryFunctor2( const T & x, BinaryFunction func )
1542  : value( x ), f( func ) {}
1543  inline T operator () ( const T & y ) const
1544  { return f( value, y ); }
1545  T value;
1546  BinaryFunction f;
1547  };
1552  template <typename T, typename U>
1553  class Scaler
1554  {
1555  public:
1556  inline Scaler( U x ) : scale( x ) {}
1557  inline T operator () ( const T & x ) const
1558  { return (T) ( x * scale ); }
1559  U scale;
1560  };
1561 
1562  template<> inline cfloat
1563  Scaler<cfloat, double>::operator () ( const cfloat & x ) const
1564  {
1565  return x * (float) scale;
1566  }
1567 
1568  template<> inline cfloat
1569  Scaler<cfloat, long>::operator () ( const cfloat & x ) const
1570  {
1571  return x * (float) scale;
1572  }
1573 
1574  template<> inline cdouble
1575  Scaler<cdouble, float>::operator () ( const cdouble & x ) const
1576  {
1577  return x * (double) scale;
1578  }
1579 
1580  template<> inline cdouble
1581  Scaler<cdouble, long>::operator () ( const cdouble & x ) const
1582  {
1583  return x * (double) scale;
1584  }
1585 
1590  template <typename T, typename U>
1591  class Divider
1592  {
1593  public:
1594  inline Divider( U x ) : divisor( x ) {}
1595  inline T operator () ( const T & x ) const
1596  { return (T) ( x / divisor ); }
1597  U divisor;
1598  };
1599 
1600  template<> inline cfloat
1601  Divider<cfloat, double>::operator () ( const cfloat & x ) const
1602  {
1603  return x * (float) ( 1. / divisor );
1604  }
1605 
1606  template<> inline cfloat
1607  Divider<cfloat, long>::operator () ( const cfloat & x ) const
1608  {
1609  return x * (float) ( 1. / (double) divisor );
1610  }
1611 
1612  template<> inline cdouble
1613  Divider<cdouble, double>::operator () ( const cdouble & x ) const
1614  {
1615  return x * ( 1. / divisor );
1616  }
1617 
1618  template<> inline cdouble
1619  Divider<cdouble, float>::operator () ( const cdouble & x ) const
1620  {
1621  return x * (double) ( 1. / divisor );
1622  }
1623 
1624  template<> inline cdouble
1625  Divider<cdouble, long>::operator () ( const cdouble & x ) const
1626  {
1627  return x * (double) ( 1. / divisor );
1628  }
1629 
1630  template <typename T, bool Scalar>
1631  class VolumeUtilBase {
1632  };
1633 
1635  template <typename T>
1636  class VolumeUtilBase<T, true> {
1637  public :
1640  static T min( const Volume<T> & o );
1641 
1644  static T max( const Volume<T> & o );
1645  };
1646 
1648  template <typename T>
1649  class VolumeUtilBase<T, false> {
1650  };
1651 
1653  template <typename T>
1654  class VolumeUtil: public carto::VolumeUtilBase<T, carto::DataTypeTraits<T>::is_scalar>
1655  {
1656  public:
1658  template <class UnaryFunction> static
1659  VolumeRef<T> apply( UnaryFunction f, const VolumeRef<T> & o );
1661 // template <class BinaryFunction> static
1662 // VolumeRef<T> apply( BinaryFunction f, const VolumeRef<T> & o1,
1663 // const VolumeRef<T> & o2 );
1666  template <class UnaryFunction> static
1667  void selfApply( UnaryFunction f, VolumeRef<T> & o );
1670  template <class BinaryFunction> static
1671  void selfApply( BinaryFunction f, VolumeRef<T> & o1,
1672  const VolumeRef<T> & o2 );
1683  template <class BinaryFunction> static
1684  T accumulate( BinaryFunction f, const VolumeRef<T> & o2, T initial )
1685  { return accumulate( f, *o2, initial ); }
1686  template <class BinaryFunction> static
1687  T accumulate( BinaryFunction f, const Volume<T> & o2, T initial );
1688  };
1689 
1690 
1691  template <typename T>
1692  inline
1693  void sort( Volume<T>& thing, bool ascending = true )
1694  {
1695  int i, j, n = thing.getSizeX();
1696  T tmp;
1697 
1698  if( ascending )
1699  {
1700  for( j=1; j<n; ++j )
1701  {
1702  tmp = thing.at( j );
1703  i = j - 1;
1704  while( i >= 0L && thing.at( i ) > tmp )
1705  {
1706  thing.at( i + 1 ) = thing.at( i );
1707  i--;
1708  }
1709  thing.at( i + 1 ) = tmp;
1710  }
1711  }
1712  else
1713  {
1714  for( j=1; j<n; ++j )
1715  {
1716  tmp = thing.at( j );
1717  i = j - 1;
1718  while( i >= 0L && thing.at( i ) < tmp )
1719  {
1720  thing.at( i + 1 ) = thing.at( i );
1721  i--;
1722  }
1723  thing.at( i + 1 ) = tmp;
1724  }
1725  }
1726  }
1727 
1728 
1729  template <typename T>
1730  inline
1731  void sort( rc_ptr<Volume<T> >& thing, bool ascending = true )
1732  {
1733  sort( *thing, ascending );
1734  }
1735 
1736 
1737 } // namespace carto
1738 
1739 #endif // CARTODATA_VOLUME_VOLUMEUTIL_H
1740 
carto::Divider
Divider functor.
Definition: volumeutil.h:1592
carto::Divider::operator()
T operator()(const T &x) const
Definition: volumeutil.h:1595
carto::Headered::header
const PropertySet & header() const
carto::NDIterator_base::ended
bool ended() const
carto::Scaler
Scaler functor.
Definition: volumeutil.h:1554
carto::Scaler::operator()
T operator()(const T &x) const
Definition: volumeutil.h:1557
carto::UnaryFromConstantBinaryFunctor2::UnaryFromConstantBinaryFunctor2
UnaryFromConstantBinaryFunctor2(const T &x, BinaryFunction func)
Definition: volumeutil.h:1541
carto::UnaryFromConstantBinaryFunctor::UnaryFromConstantBinaryFunctor
UnaryFromConstantBinaryFunctor(const T &x, BinaryFunction func)
Definition: volumeutil.h:1524
carto::VolumeProxy::getSizeY
int getSizeY() const
Definition: volumeproxy.h:99
carto::VolumeProxy::getSizeX
int getSizeX() const
Definition: volumeproxy.h:89
carto::VolumeProxy::getSize
std::vector< int > getSize() const
get the 4 dimensions in a vector
Definition: volumeproxy.h:129
carto::VolumeProxy::getSizeZ
int getSizeZ() const
Definition: volumeproxy.h:109
carto::VolumeProxy::getSizeT
int getSizeT() const
Definition: volumeproxy.h:119
carto::VolumeProxy::copyHeaderFrom
virtual void copyHeaderFrom(const PropertySet &other)
copy properties from other to this, avoiding forbidden properties like size.
Definition: volumeproxy.h:139
carto::VolumeRef
Convenient handle for a Volume - this is normally the entry point for all volumes handling.
Definition: volumeref.h:60
carto::VolumeUtil
Volume utilities classes.
Definition: volumeutil.h:1655
carto::VolumeUtil::accumulate
static T accumulate(BinaryFunction f, const VolumeRef< T > &o2, T initial)
Apply a binary function to each voxel of the volume, with a "current" value as second argument.
Definition: volumeutil.h:1684
carto::VolumeUtil::selfApply
static void selfApply(UnaryFunction f, VolumeRef< T > &o)
applies a binary function to each voxel of a pair of volumes
Definition: volumeutil_d.h:299
carto::VolumeUtil::apply
static VolumeRef< T > apply(UnaryFunction f, const VolumeRef< T > &o)
applies a unary function to each voxel of a volume
Definition: volumeutil_d.h:74
carto::Volume
N-D Volume main class.
Definition: volumebase.h:119
carto::Volume::allocatorContext
const AllocatorContext & allocatorContext() const
returns volume's AllocatorContext
Definition: volumebase_d.h:521
carto::Volume::setPosInRefVolume
void setPosInRefVolume(const Position4Di &pos)
Set position in parent volume.
Definition: volumebase_d.h:677
carto::Volume::reallocate
virtual void reallocate(int sizeX=1, int sizeY=1, int sizeZ=1, int sizeT=1, bool keepcontents=false, const AllocatorContext &allocatorContext=AllocatorContext(), bool allocate=true, const std::vector< size_t > *strides=0)
allows resizing and changing allocator
Definition: volumebase_d.h:1165
carto::Volume::posInRefVolume
const Position & posInRefVolume() const
Get position in parent volume.
Definition: volumebase_d.h:533
carto::Volume::getStrides
std::vector< size_t > getStrides() const
Get strides for the volume.
Definition: volumebase_d.h:726
carto::Volume::setRefVolume
void setRefVolume(const rc_ptr< Volume< T > > &refvol)
Set parent volume.
Definition: volumebase_d.h:699
carto::Volume::refVolume
rc_ptr< Volume< T > > refVolume() const
Get parent volume.
Definition: volumebase_d.h:527
carto::Volume::getBorders
std::vector< int > getBorders() const
Get borders for the volume.
Definition: volumebase_d.h:707
carto::Volume::at
const T & at(long x, long y=0, long z=0, long t=0) const
Definition: volumebase_d_inline.h:48
carto::const_line_NDIterator
carto::const_line_NDIterator::line_length
long line_length() const
carto::const_line_NDIterator::inc_line_ptr
void inc_line_ptr(const T *&p) const
carto::line_NDIterator
carto::line_NDIterator::line_length
long line_length() const
carto::line_NDIterator::inc_line_ptr
void inc_line_ptr(T *&p) const
carto::rc_ptr
carto::rc_ptr::reset
void reset(T *p=NULL)
carto::rc_ptr::get
T * get() const
carto::ref::release
T * release()
carto::volumeutil::applyTowards
Volume< OUTP > & applyTowards(const Volume< T > &vol, Volume< OUTP > &dst, UnaryFunction func)
Apply a function to all the elements of a volume (already allocated output version)
Definition: volumeutil.h:603
carto::volumeutil::accumulate
OUTP accumulate(const Volume< T > &vol, BinaryFunction func, OUTP initial)
Accumulation over a volume.
Definition: volumeutil.h:814
carto::volumeutil::accumulate
OUTP accumulate(const rc_ptr< Volume< T > > &vol, BinaryFunction func, OUTP initial)
Definition: volumeutil.h:832
carto::volumeutil::minSize
std::vector< int > minSize(const std::vector< int > &s1, const std::vector< int > &s2)
return the smallest dimensions between s1 and s2
Definition: volumeutil.h:497
carto::volumeutil::selfApply
Volume< T > & selfApply(Volume< T > &vol, UnaryFunction func)
Apply a function to all the elements of a volume (in place version)
Definition: volumeutil.h:567
carto::volumeutil::maxSize
std::vector< int > maxSize(const std::vector< int > &s1, const std::vector< int > &s2)
return the largest dimensions between s1 and s2
Definition: volumeutil.h:479
carto::volumeutil::apply
Volume< typename UnaryFunction::result_type > apply(const Volume< T > &vol, UnaryFunction func)
Apply a function to all the elements of a volume (create output version)
Definition: volumeutil.h:518
carto::volumeutil::sameSize
bool sameSize(const std::vector< int > &s1, const std::vector< int > &s2)
Used by the actual Volume and VolumeRef operators It allows to keep the loops in one place and to spe...
Definition: volumeutil.h:460
carto::any
bool any(const Volume< T > &vol)
Returns true if at least one value compares to true.
Definition: volumeutil.h:1218
carto::valuesIn
Volume< bool > valuesIn(const Volume< T > &volume, const U &set)
Find values contained in a set.
Definition: volumeutil.h:1250
carto::invertMinMax
Volume< T > invertMinMax(const Volume< T > &volume)
Invert a volume between its min and max values.
Definition: volumeutil.h:1347
carto::isSameVolumeSize
bool isSameVolumeSize(const Volume< T > &v1, const Volume< T > &v2)
Definition: volumeutil.h:420
carto::deepcopy
Volume< T > deepcopy(const Volume< T > &src, const std::vector< int > &size=std::vector< int >())
Performs a copy of the data (not only a reference copy) The whole view hierarchy is fully duplicated.
Definition: volumeutil.h:877
carto::conditionalSet
void conditionalSet(Volume< T > &volume, const Volume< U > &condition, const T &value)
Conditionally fill a volume.
Definition: volumeutil.h:1314
carto::transfer
void transfer(const Volume< T > &src, Volume< T > &dst)
Transfer data from an allocated volume to another Allocated sizes must be equal (same number of voxel...
Definition: volumeutil.h:847
carto::copyStructure
Volume< T > copyStructure(const Volume< T > &src)
Performs a copy of the view structure without transfering the data.
Definition: volumeutil.h:1030
carto::min
T min(const rc_ptr< Volume< T > > &vol)
Definition: volumeutil.h:1143
carto::copy
Volume< T > copy(const Volume< T > &src)
Performs a copy of the data (not only a reference copy) Only the data inside the view is copied.
Definition: volumeutil.h:993
carto::min
T min(const Volume< T > &vol)
Returns the minimum value of the volume.
Definition: volumeutil.h:1133
carto::max
T max(const Volume< T > &vol)
Returns the maximum value of the volume.
Definition: volumeutil.h:1153
carto::sum
DataTypeTraits< T >::LongType sum(const Volume< T > &vol)
Returns the sum of the volume values.
Definition: volumeutil.h:1173
carto::sort
void sort(Volume< T > &thing, bool ascending=true)
Definition: volumeutil.h:1693
carto::matrix_product
VolumeRef< T > matrix_product(const Volume< T > &v1, const Volume< T > &v2)
matrix product
Definition: volumeutil_d.h:384
carto::transpose
VolumeRef< T > transpose(const Volume< T > &v)
transpose
Definition: volumeutil_d.h:416
carto::all
bool all(const Volume< T > &vol)
Returns true if all values compare to true.
Definition: volumeutil.h:1204
carto::undiag
rc_ptr< Volume< T > > undiag(const rc_ptr< Volume< T > > &thing)
Definition: volumeutil.h:353
carto::setBorders
void setBorders(Volume< T > &volume, const typename Volume< T >::Position4Di &top, const typename Volume< T >::Position4Di &bottom=typename Volume< T >::Position4Di(-1, -1, -1, -1))
Set border width.
Definition: volumeutil.h:1369
carto::setMinBorders
void setMinBorders(Volume< T > &volume, const typename Volume< T >::Position4Di &top, const typename Volume< T >::Position4Di &bottom=typename Volume< T >::Position4Di(-1, -1, -1, -1))
Set border width.
Definition: volumeutil.h:1439
carto::diag
rc_ptr< Volume< T > > diag(const rc_ptr< Volume< T > > &thing)
Build a diagonal matrix from a line (1D) vector.
Definition: volumeutil.h:334
carto::valuesNotIn
Volume< bool > valuesNotIn(const Volume< T > &volume, const U &set)
Find values not contained in a set.
Definition: volumeutil.h:1282
carto::max
T max(const rc_ptr< Volume< T > > &vol)
Definition: volumeutil.h:1163
carto::internal::changeIf::operator()
bool operator()(const T &x, const U &y)
Definition: volumeutil.h:1303
carto::internal::changeIf::changeIf
changeIf(const T &value)
Definition: volumeutil.h:1302
carto::internal::inSet::operator()
bool operator()(const T &x, const U &y)
Definition: volumeutil.h:1238
carto::internal::invMinMax::invMinMax
invMinMax(const T &min, const T &max)
Definition: volumeutil.h:1330
carto::internal::notInSet::operator()
bool operator()(const T &x, const U &y)
Definition: volumeutil.h:1270
carto::volumeutil::logical_and
Definition: functional.h:924
carto::volumeutil::logical_or
Definition: functional.h:933
carto::volumeutil::select_max
Definition: functional.h:1018
carto::volumeutil::select_min
Definition: functional.h:1008
cdouble
std::complex< double > cdouble
cfloat
std::complex< float > cfloat