cartodata  4.5.0
volumeutil_d.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_D_H
35 #define CARTODATA_VOLUME_VOLUMEUTIL_D_H
36 
46 #include <limits>
47 
48 /* bug in gcc 4.2, accumulate doesn't accept std::max<T> or std::min<T> as
49  argument, but accepts a function without a namespace */
50 /* in C++11 (g++ --std=c++11) std::min and std::max<T> cannot be passed to a
51  function also - because there are other variants, like
52  max( std::initializer_list<T> ilist ). */
53 
54 namespace
55 {
56  template <typename T>
57  inline
58  T internal_max( const T x, const T y )
59  { return std::max<T>( x, y ); }
60 
61  template <typename T>
62  inline
63  T internal_min( const T x, const T y )
64  { return std::min<T>( x, y ); }
65 }
66 
67 
68 namespace carto
69 {
70 
71  template <typename T> template <class UnaryFunction>
72  VolumeRef<T> VolumeUtil<T>::apply( UnaryFunction f, const VolumeRef<T> & o )
73  {
74  VolumeRef<T> res
75  ( new Volume<T>( o->getSizeX(), o->getSizeY(), o->getSizeZ(),
76  o->getSizeT() ) );
77  res->header() = o->header();
78 
79  unsigned x, nx = o->getSizeX(), y, ny = o->getSizeY(),
80  z, nz = o->getSizeZ(), t, nt = o->getSizeT();
81  T *op, *rp;
82  for( t=0; t<nt; ++t )
83  for( z=0; z<nz; ++z )
84  for( y=0; y<ny; ++y )
85  {
86  op = &o->at( 0, y, z, t );
87  rp = &res->at( 0, y, z, t );
88  for( x=0; x<nx; ++x )
89  *rp++ = f( *op++ );
90  }
91  return res;
92  }
93 
94 
95  namespace internal
96  {
97 
98  template <typename T> inline T _neutral()
99  {
100  return T();
101  }
102 
103  template<> inline int8_t _neutral<int8_t>()
104  {
105  return 0;
106  }
107 
108  template<> inline uint8_t _neutral<uint8_t>()
109  {
110  return 0;
111  }
112 
113  template<> inline int16_t _neutral<int16_t>()
114  {
115  return 0;
116  }
117 
118  template<> inline uint16_t _neutral<uint16_t>()
119  {
120  return 0;
121  }
122 
123  template<> inline int32_t _neutral<int32_t>()
124  {
125  return 0;
126  }
127 
128  template<> inline uint32_t _neutral<uint32_t>()
129  {
130  return 0;
131  }
132 
133  template<> inline float _neutral<float>()
134  {
135  return 0;
136  }
137 
138  template<> inline double _neutral<double>()
139  {
140  return 0;
141  }
142 
143  }
144 
145 
146  template <typename T> template <class BinaryFunction>
148  const VolumeRef<T> & o1,
149  const VolumeRef<T> & o2 )
150  {
151  int x, nx = o1->getSizeX(), y, ny = o1->getSizeY(),
152  z, nz = o1->getSizeZ(), t, nt = o1->getSizeT(), nx2, ny2, nz2, nt2, nx3;
153  if( o2->getSizeX() < nx )
154  {
155  nx2 = nx;
156  nx = o2->getSizeX();
157  }
158  else
159  nx2 = o2->getSizeX();
160  if( o2->getSizeY() < ny )
161  {
162  ny2 = ny;
163  ny = o2->getSizeY();
164  }
165  else
166  ny2 = o2->getSizeY();
167  if( o2->getSizeZ() < nz )
168  {
169  nz2 = nz;
170  nz = o2->getSizeZ();
171  }
172  else
173  nz2 = o2->getSizeZ();
174  if( o2->getSizeT() < nt )
175  {
176  nt2 = nt;
177  nt = o2->getSizeT();
178  }
179  else
180  nt2 = o2->getSizeT();
181 
182  VolumeRef<T> res( new Volume<T>( nx2, ny2, nz2, nt2 ) );
183  res->header() = o1->header();
184 
185  bool x1, y1, z1, x2, y2, z2;
186  T *o1p, *o2p, *rp;
187  for( t=0; t<nt2; ++t )
188  {
189  if( t >= nt )
190  if( nt == o1->getSizeT() )
191  {
192  z1 = false;
193  z2 = true;
194  }
195  else
196  {
197  z1 = true;
198  z2 = false;
199  }
200  else
201  {
202  z1 = true;
203  z2 = true;
204  }
205 
206  for( z=0; z<nz2; ++z )
207  {
208  y1 = z1;
209  y2 = z2;
210  if( z >= nz )
211  {
212  if( nz == o1->getSizeZ() )
213  y1 = false;
214  else
215  y2 = false;
216  }
217 
218  for( y=0; y<ny2; ++y )
219  {
220  x1 = y1;
221  x2 = y2;
222  if( y >= ny )
223  {
224  if( ny == o1->getSizeY() )
225  x1 = false;
226  else
227  x2 = false;
228  }
229  if( x1 )
230  o1p = &o1->at( 0, y, z, t );
231  else
232  o1p = 0;
233  if( x2 )
234  o2p = &o2->at( 0, y, z, t );
235  else
236  o2p = 0;
237  rp = &res->at( 0, y, z, t );
238 
239  /*
240  std::cout << "pos: " << t << ", " << z << ", " << y
241  << ", x1: " << x1 << ", " << y1 << ", " << z1
242  << ", x2: " << x2 << ", " << y2 << ", " << z2
243  << std::endl;
244  */
245 
246  x = 0;
247  if( x1 && x2 )
248  {
249  for( ; x<nx; ++x )
250  *rp++ = f( *o1p++, *o2p++ );
251  nx3 = nx2;
252  if( nx == o1->getSizeX() )
253  x1 = false;
254  else
255  x2 = false;
256  }
257  else if( x1 )
258  nx3 = o1->getSizeX();
259  else if( x2 )
260  nx3 = o2->getSizeX();
261  else
262  nx3 = 0;
263  if( !x1 && x2 )
264  for( ; x<nx3; ++x )
265  *rp++ = f( internal::_neutral<T>(), *o2p++ );
266  else if( x1 && !x2 )
267  for( ; x<nx3; ++x )
268  *rp++ = f( *o1p++, internal::_neutral<T>() );
269  for( ; x<nx2; ++x )
270  *rp++ = f( internal::_neutral<T>(),
271  internal::_neutral<T>() );
272  }
273  }
274  }
275  return res;
276  }
277 
278 
279  template <typename T> template <class UnaryFunction>
280  void VolumeUtil<T>::selfApply( UnaryFunction f, VolumeRef<T> & o )
281  {
282  int x, nx = o->getSizeX(), y, ny = o->getSizeY(),
283  z, nz = o->getSizeZ(), t, nt = o->getSizeT();
284  T *op;
285  for( t=0; t<nt; ++t )
286  for( z=0; z<nz; ++z )
287  for( y=0; y<ny; ++y )
288  {
289  op = &o->at( 0, y, z, t );
290  for( x=0; x<nx; ++x, ++op )
291  *op = f( *op );
292  }
293  }
294 
295 
296  template <typename T> template <class BinaryFunction>
297  void VolumeUtil<T>::selfApply( BinaryFunction f, VolumeRef<T> & o1,
298  const VolumeRef<T> & o2 )
299  {
300  int x, nx = o1->getSizeX(), y, ny = o1->getSizeY(),
301  z, nz = o1->getSizeZ(), t, nt = o1->getSizeT();
302 
303  if( o2->getSizeX() < nx )
304  nx = o2->getSizeX();
305  if( o2->getSizeY() < ny )
306  ny = o2->getSizeY();
307  if( o2->getSizeZ() < nz )
308  nz = o2->getSizeZ();
309  if( o2->getSizeT() < nt )
310  nt = o2->getSizeT();
311 
312  T *o1p, *o2p;
313  for( t=0; t<nt; ++t )
314  for( z=0; z<nz; ++z )
315  for( y=0; y<ny; ++y )
316  {
317  o1p = &o1->at( 0, y, z, t );
318  o2p = &o2->at( 0, y, z, t );
319  for( x=0; x<nx; ++x, ++o1p )
320  *o1p = f( *o1p, *o2p++ );
321  }
322  }
323 
324 
325  template <typename T> template <class BinaryFunction>
326  T VolumeUtil<T>::accumulate( BinaryFunction f,
327  const Volume<T> & o, T initial )
328  {
329  T res = initial;
330 
331  unsigned x, nx = o.getSizeX(), y, ny = o.getSizeY(),
332  z, nz = o.getSizeZ(), t, nt = o.getSizeT();
333  const T *op;
334  for( t=0; t<nt; ++t )
335  for( z=0; z<nz; ++z )
336  for( y=0; y<ny; ++y )
337  {
338  op = &o.at( 0, y, z, t );
339  for( x=0; x<nx; ++x )
340  res = f( *op++, res );
341  }
342  return res;
343  }
344 
345  template <typename T>
347  {
349  internal_min<T>, o, std::numeric_limits<T>::max() );
350  }
351 
352  template <typename T>
354  {
356  internal_max<T>, o, -std::numeric_limits<T>::max() );
357  }
358 
359 
360 }
361 
362 #endif
363 
4D Volume main class
static VolumeRef< T > apply(UnaryFunction f, const VolumeRef< T > &o)
applies a unary function to each voxel of a volume
Definition: volumeutil_d.h:72
int getSizeT() const
Definition: volumeproxy.h:118
const T & at(long x, long y=0, long z=0, long t=0) const
Convenient handle for a Volume - this is normally the entry point for all volumes handling...
int getSizeZ() const
Definition: volumeproxy.h:108
int getSizeX() const
Definition: volumeproxy.h:88
const T & at(long x, long y=0, long z=0, long t=0) const
const PropertySet & header() const
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:1267
static void selfApply(UnaryFunction f, VolumeRef< T > &o)
same as apply() except that the input volume is used to store the result
Definition: volumeutil_d.h:280
int getSizeY() const
Definition: volumeproxy.h:98