newimage
 
Loading...
Searching...
No Matches
newimage.h
1/* Templated image storage class
2
3 Mark Jenkinson, FMRIB Image Analysis Group
4
5 Copyright (C) 2000 University of Oxford */
6
7/* CCOPYRIGHT */
8
9#if !defined(__newimage_h)
10#define __newimage_h
11
12#define volume4D volume
13
14#include <algorithm>
15#include <cstdlib>
16#include <iostream>
17#include <limits>
18#include <stdexcept>
19#include <string>
20#include <type_traits>
21#include <vector>
22
23#include "armawrap/newmatap.h"
24#include "positerators.h"
25#include "NewNifti/NewNifti.h"
26#include "miscmaths/miscmaths.h"
27#include "miscmaths/kernel.h"
28#include "miscmaths/splinterpolator.h"
29#include "utils/threading.h"
30
31
32namespace NEWIMAGE {
33
34 #define SUBSET -1
35
36 const bool USEASSERT=false; // used inside imthrow to determine type of action
37
38 void imthrow(const std::string& msg, int nierrnum, const bool quiet=false);
39
40 enum extrapolation { zeropad, constpad, extraslice, mirror, periodic,
41 boundsassert, boundsexception, userextrapolation };
42
43 enum interpolation { nearestneighbour, trilinear, sinc, userkernel,
44 userinterpolation, spline };
45
46 enum threshtype { inclusive , exclusive };
47
48 enum constructionMode { CLONE, TEMPLATE, ALIAS };
49
50#define FSL_TYPE_ANALYZE 10
51#define FSL_TYPE_NIFTI 1
52#define FSL_TYPE_NIFTI2 2
53#define FSL_TYPE_NIFTI_PAIR 11
54#define FSL_TYPE_NIFTI2_PAIR 12
55#define FSL_TYPE_ANALYZE_GZ 100
56#define FSL_TYPE_NIFTI_GZ 101
57#define FSL_TYPE_NIFTI2_GZ 102
58#define FSL_TYPE_NIFTI_PAIR_GZ 111
59#define FSL_TYPE_NIFTI2_PAIR_GZ 112
60
61class fileFormats {
62public:
63 static const std::vector<std::pair<std::string,int> > allFormats() {
64 const std::pair<std::string,int> allOutputFormats[10] = {
65 std::pair<std::string,int>("ANALYZE",FSL_TYPE_ANALYZE),
66 std::pair<std::string,int>("ANALYZE_GZ",FSL_TYPE_ANALYZE_GZ),
67 std::pair<std::string,int>("NIFTI",FSL_TYPE_NIFTI),
68 std::pair<std::string,int>("NIFTI_GZ",FSL_TYPE_NIFTI_GZ),
69 std::pair<std::string,int>("NIFTI_STD::PAIR",FSL_TYPE_NIFTI_PAIR),
70 std::pair<std::string,int>("NIFTI_STD::PAIR_GZ",FSL_TYPE_NIFTI_PAIR_GZ),
71 std::pair<std::string,int>("NIFTI2",FSL_TYPE_NIFTI2),
72 std::pair<std::string,int>("NIFTI2_GZ",FSL_TYPE_NIFTI2_GZ),
73 std::pair<std::string,int>("NIFTI2_STD::PAIR",FSL_TYPE_NIFTI2_PAIR),
74 std::pair<std::string,int>("NIFTI2_STD::PAIR_GZ",FSL_TYPE_NIFTI2_PAIR_GZ),
75 };
76 std::vector<std::pair<std::string,int> > formats(10);
77 for(int i=0;i<10;i++)
78 formats[i]=allOutputFormats[i];
79 return formats;
80 }
81};
82
83#define FSL_RADIOLOGICAL -1
84#define FSL_NEUROLOGICAL 1
85#define FSL_INCONSISTENT 0
86#define FSL_ZERODET -101
87
88template<class T> class ShadowVolume;
89template<class T> class volume;
90
91template <class T>
92int readGeneralVolume(volume<T>& target, const std::string& filename,
93 short& dtype, const bool swap2radiological=true,
94 int64_t x0 = -1, int64_t y0 = -1, int64_t z0 = -1, int64_t t0 = -1, int64_t d50 = -1, int64_t d60 = -1, int64_t d70 = -1,
95 int64_t x1 = -1, int64_t y1 = -1, int64_t z1 = -1, int64_t t1 = -1, int64_t d51 = -1, int64_t d61 = -1, int64_t d71 = -1,
96 const bool readAs4D=false);
97
98#pragma interface
99 template <class T>
100 class volume {
101 private:
102 T* Data;
103 T* DataEnd;
104 mutable bool data_owner;
105 mutable double maskDelimiter;
106 int64_t nElements;
107 int64_t nThreads;
108
109 int64_t ColumnsX;
110 int64_t RowsY;
111 int64_t SlicesZ;
112 int64_t dim4;
113 int64_t dim5;
114 int64_t dim6;
115 int64_t dim7;
116 std::vector<int64_t> originalSizes;
117
118 float Xdim;
119 float Ydim;
120 float Zdim;
121 float p_TR;
122 float pxdim5;
123 float pxdim6;
124 float pxdim7;
125
126 // the following matrices map voxel coords to mm coords (nifti conventions)
127 NEWMAT::Matrix StandardSpaceCoordMat;
128 NEWMAT::Matrix RigidBodyCoordMat;
129 int StandardSpaceTypeCode;
130 int RigidBodyTypeCode;
131
132 mutable int IntentCode;
133 mutable float IntentParam1;
134 mutable float IntentParam2;
135 mutable float IntentParam3;
136
137 mutable int SliceOrderingCode;
138
139 int64_t no_voxels;
140
141 mutable SPLINTERPOLATOR::Splinterpolator<T> splint; // Object used to perform spline interpolation
142 mutable std::mutex splinecoef_mutex; // Used to lock recalculation of spline coefficients
143 mutable int splineorder; // Spline-order (typically 3 for cubic splines)
144 mutable bool splineuptodate; // Indicates need to recalculate spline coefficients
145
146 mutable MISCMATHS::kernel interpkernel;
147 mutable extrapolation p_extrapmethod;
148
149 mutable interpolation p_interpmethod;
150 mutable T (*p_userextrap)(const volume<T>& , int64_t, int64_t, int64_t);
151 mutable float (*p_userinterp)(const volume<T>& , float, float, float);
152 mutable T padvalue; // was p_padval in volume4D
153
154 mutable T extrapval; // the reference target for all extrapolations
155 mutable std::vector<bool> ep_valid; // Indicates if extrapolation can produce "valid" values.
156
157 mutable float displayMaximum;
158 mutable float displayMinimum;
159 mutable char auxFile[24];
160
161 // Internal functions
162 inline T* basicptr(int64_t x, int64_t y, int64_t z) {
163 return (nsfbegin() + (z*RowsY + y)*ColumnsX + x);
164 }
165 inline T* basicptr(int64_t x, int64_t y, int64_t z) const
166 { return (Data + (z*RowsY + y)*ColumnsX + x); }
167 int initialize(int64_t xsize, int64_t ysize, int64_t zsize, T *d, bool d_owner, int64_t nthreads); //3D
168 int initialize(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, T *d, bool d_owner, int64_t nthreads); //4D
169 protected:
170 virtual int initialize(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, int64_t d5, int64_t d6, int64_t d7, T *d, bool d_owner, int64_t nthreads); //Master 7D
171 private:
172 void setdefaultproperties();
173 void enforcelimits(std::vector<int>& lims) const;
174 const T& extrapolate(int64_t x, int64_t y, int64_t z) const;
175 float kernelinterpolation(const float x, const float y,
176 const float z) const;
177 float splineinterpolate(float x, float y, float z) const;
178 float spline_interp1partial(float x, float y, float z, int dir, float *deriv) const;
179 float spline_interp3partial(float x, float y, float z, float *dfdx, float *dfdy, float *dfdz) const;
180
181 template <class S, class D> friend
182 void copybasicproperties(const volume<S>& source, volume<D>& dest);
183
184 template <class S> friend
185 int readGeneralVolume(volume<S>& target, const std::string& filename,
186 short& dtype, const bool swap2radiological,
187 int64_t x0, int64_t y0, int64_t z0, int64_t t0, int64_t d50, int64_t d60, int64_t d70,
188 int64_t x1, int64_t y1, int64_t z1, int64_t t1, int64_t d51, int64_t d61, int64_t d71,
189 const bool readAs4D);
190
191 template <class S> friend
192 int read_volume_hdr_only(volume<S>&, const std::string&);
193
194 void basic_swapdimensions(int dim1, int dim2, int dim3, bool keepLRorder, const bool headerOnly=false);
195
196#ifdef EXPOSE_TREACHEROUS
197 public:
198#endif
199 // sampling_mat should now be avoided - use newimagevox2mm_mat instead
200 std::vector<NiftiIO::NiftiExtension> extensions;
201 bool RadiologicalFile;
202 NEWMAT::Matrix sampling_mat() const;
203 void set_sform(int sform_code, const NEWMAT::Matrix& snewmat);
204 void set_qform(int qform_code, const NEWMAT::Matrix& qnewmat);
205 int left_right_order() const; // see also newimagevox2mm_mat()
206 void swapLRorder(const bool headerOnly=false);
207 void setLRorder(int LRorder);
208 void makeradiological(const bool headerOnly=false);
209 void makeneurological();
210 NEWMAT::Matrix sform_mat() const { return StandardSpaceCoordMat; }
211 int sform_code() const { return StandardSpaceTypeCode; }
212 NEWMAT::Matrix qform_mat() const { return RigidBodyCoordMat; }
213 int qform_code() const { return RigidBodyTypeCode; }
214 void set_data_owner(bool d_owner) const { data_owner=d_owner; }
215 public:
216 typedef T* nonsafe_fast_iterator;
217 inline nonsafe_fast_iterator nsfbegin()
218 { this->invalidateSplines(); return Data; }
219 inline nonsafe_fast_iterator nsfend()
220 { return DataEnd; }
221 public:
222 // CONSTRUCTORS AND DESTRUCTORS (including copy, = and reinitialize)
223 volume(Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
224 volume(const volume<T>& source, const bool copyData);
225 volume(const volume<T>& source, const constructionMode mode=CLONE);
226 volume(int64_t xsize, int64_t ysize, int64_t zsize, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
227 volume(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
228 volume(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, int64_t d5, int64_t d6, int64_t d7, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
229 virtual ~volume();
230 void destroy();
231 const volume<T>& operator=(const volume<T>& that) { return this->equals(that); }
232
233 virtual const volume<T>& equals(const volume<T>& );
234
235 void reinitialize(const volume<T>& source, const constructionMode mode=CLONE);
236 int reinitialize(int64_t xsize, int64_t ysize, int64_t zsize, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
237 int reinitialize(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize,T *d=nullptr, bool d_owner=true, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
238 int reinitialize(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, int64_t d5, int64_t d6, int64_t d7, Utilities::NoOfThreads nt = Utilities::NoOfThreads(1));
239 void copyproperties(const volume<T>& source);
240 int copydata(const volume<T>& source, const bool isOwner=true);
241
242 // BASIC PROPERTIES
243 inline int64_t xsize() const { return ColumnsX; }
244 inline int64_t ysize() const { return RowsY; }
245 inline int64_t zsize() const { return SlicesZ; }
246 inline int64_t tsize() const { return dim4; }
247 inline int64_t size1() const { return xsize(); }
248 inline int64_t size2() const { return ysize(); }
249 inline int64_t size3() const { return zsize(); }
250 inline int64_t size4() const { return tsize(); }
251 inline int64_t size5() const { return dim5; }
252 inline int64_t size6() const { return dim6; }
253 inline int64_t size7() const { return dim7; }
254 int64_t size(int n) const;
255
256 int dimensionality() const;
257 void throwsIfNot3D() const { if (this->dimensionality()>3) imthrow("3D only method called by higher-dimensional volume.",75); }
258
259 inline int64_t maxx() const { return xsize()-1; }
260 inline int64_t maxy() const { return ysize()-1; }
261 inline int64_t maxz() const { return zsize()-1; }
262 inline int64_t maxt() const { return tsize()-1; }
263 inline int64_t minx() const { return 0; }
264 inline int64_t miny() const { return 0; }
265 inline int64_t minz() const { return 0; }
266 inline int64_t mint() const { return 0; }
267
268 inline float xdim() const { return Xdim; }
269 inline float ydim() const { return Ydim; }
270 inline float zdim() const { return Zdim; }
271 inline float tdim() const { return p_TR; }
272 inline float TR() const { return p_TR; }
273 inline float pixdim1() const { return Xdim; }
274 inline float pixdim2() const { return Ydim; }
275 inline float pixdim3() const { return Zdim; }
276 inline float pixdim4() const { return p_TR; }
277 inline float pixdim5() const { return pxdim5; }
278 inline float pixdim6() const { return pxdim6; }
279 inline float pixdim7() const { return pxdim7; }
280
281 inline float getDisplayMaximum() const { return displayMaximum; }
282 inline float getDisplayMinimum() const { return displayMinimum; }
283 inline std::string getAuxFile() const { return std::string(auxFile); }
284
285 inline bool hasSplines() const { return(splint.Valid()); }
286
287 inline int64_t nthreads() const { return(nThreads); }
288
289 void setxdim(float x) { Xdim = fabs(x); }
290 void setydim(float y) { Ydim = fabs(y); }
291 void setzdim(float z) { Zdim = fabs(z); }
292 void settdim(float tr) { p_TR = fabs(tr); }
293 void setd5dim(float d5) { pxdim5 = d5; }
294 void setd6dim(float d6) { pxdim6 = d6; }
295 void setd7dim(float d7) { pxdim7 = d7; }
296 void setTR(float tr) { settdim(tr); }
297 void setdims(float x, float y, float z)
298 { setxdim(x); setydim(y); setzdim(z); }
299 void setdims(float x, float y, float z, float tr)
300 { setxdim(x); setydim(y); setzdim(z); settdim(tr); }
301 void setdims(float x, float y, float z, float tr, float d5, float d6, float d7)
302 { setxdim(x); setydim(y); setzdim(z); settdim(tr); setd5dim(d5); setd6dim(d6); setd7dim(d7); }
303 void setDisplayMaximumMinimum(const float maximum, const float minimum) const { displayMaximum=maximum; displayMinimum=minimum; }
304 void setDisplayMaximum(const float maximum) const { setDisplayMaximumMinimum(maximum,displayMinimum); }
305 void setDisplayMinimum(const float minimum) const { setDisplayMaximumMinimum(displayMaximum,minimum); }
306 void setAuxFile(const std::string fileName) { strncpy(auxFile,fileName.c_str(),24); }
307 int64_t nvoxels() const { return no_voxels; }
308 int64_t ntimepoints() const { return tsize(); }
309 int64_t totalElements() const { return nElements; }
310
311 std::vector<int64_t> ptrToCoord(size_t offset) const;
312
313 static const T smallestPositiveT;
314 // MATRIX <-> VOLUME CONVERSIONS
315 NEWMAT::ReturnMatrix vec(const volume<T>& mask) const;
316 NEWMAT::ReturnMatrix vec() const;
317 void insert_vec(const NEWMAT::ColumnVector& pvec, const volume<T>& pmask);
318 void insert_vec(const NEWMAT::ColumnVector& pvec);
319 std::vector<int64_t> labelToCoord(const int64_t label) const;
320 NEWMAT::ReturnMatrix matrix(const volume<T>& mask, bool using_mask=true) const;
321 NEWMAT::ReturnMatrix matrix(const volume<T>& mask, std::vector<int64_t>& voxelLabels) const;
322 NEWMAT::ReturnMatrix matrix() const;
323 void setmatrix(const NEWMAT::Matrix& newmatrix, const volume<T>& mask,
324 const T pad=0, bool using_mask=true);
325 void setmatrix(const NEWMAT::Matrix& newmatrix);
326 volume<int> vol2matrixkey(const volume<T>& mask); //returns a volume with numbers in relating to matrix colnumbers
327 NEWMAT::ReturnMatrix matrix2volkey(volume<T>& mask);
328
329 // SECONDARY PROPERTIES
330 // maps *NEWIMAGE* voxel coordinates to mm (consistent with FSLView mm)
331 // NB: do not try to determine left-right order from this matrix
332 // sampling_mat should now be avoided - use newimagevox2mm_mat instead
333 NEWMAT::Matrix newimagevox2mm_mat() const;
334 NEWMAT::Matrix niftivox2newimagevox_mat() const;
335
336 int intent_code() const { return IntentCode; }
337 float intent_param(int n) const;
338 void set_intent(int intent_code, float p1, float p2, float p3) const;
339
340 inline double maskThreshold() const { return maskDelimiter; }
341 T min() const;
342 T max() const;
343 double sum() const;
344 double sumsquares() const;
345 double mean() const { return sum()/(MISCMATHS::Max(1.0,(double)totalElements()));}
346 double variance() const { double n(totalElements());
347 return (n/(n-1))*(sumsquares()/n - mean()*mean()); }
348 double stddev() const { return sqrt(variance()); }
349 T robustmin() const;
350 T robustmax() const;
351 std::vector<T> robustlimits(const volume<T>& mask) const;
352 NEWMAT::ColumnVector principleaxis(int n) const;
353 NEWMAT::Matrix principleaxes_mat() const;
354 T percentile(float pvalue) const; // argument in range [0.0 , 1.0]
355 NEWMAT::ColumnVector histogram(int nbins) const;
356 NEWMAT::ColumnVector histogram(int nbins, T minval, T maxval) const;
357 NEWMAT::ColumnVector cog(const std::string& coordtype="voxel", bool abscog=false) const;
358
359 // SECONDARY PROPERTIES (using mask)
360 T min(const volume<T>& mask) const;
361 T max(const volume<T>& mask) const;
362 T min(const volume<T>& mask, std::vector<int64_t>& coords) const;
363 T max(const volume<T>& mask, std::vector<int64_t>& coords) const;
364 double sum(const volume<T>& mask) const;
365 double sumsquares(const volume<T>& mask) const;
366 double mean(const volume<T>& mask) const;
367 double variance(const volume<T>& mask) const;
368 double stddev(const volume<T>& mask) const { return sqrt(variance(mask)); }
369 T robustmin(const volume<T>& mask) const;
370 T robustmax(const volume<T>& mask) const;
371 T percentile(float pvalue, const volume<T>& mask) const; // arg in [0,1]
372 NEWMAT::ColumnVector histogram(int nbins, const volume<T>& mask) const;
373 NEWMAT::ColumnVector histogram(int nbins, T minval, T maxval, const volume<T>& mask)
374 const;
375
376 // DATA ACCESS FUNCTIONS (iterators)
377
378 typedef const T* fast_const_iterator;
379
380 inline fast_const_iterator fbegin() const { return Data; }
381 inline fast_const_iterator fend() const { return DataEnd; }
382 inline fast_const_iterator fbegin(const int64_t offset) const { return Data + offset*nvoxels(); }
383 inline fast_const_iterator fend(const int64_t offset) const { return Data + (offset+1)*nvoxels(); }
384
385 // BASIC DATA ACCESS FUNCTIONS
386 template<class A, class B, class C>
387 inline bool in_bounds( A x, B y, C z) const {
388 if ( std::is_integral<A>::value && std::is_integral<B>::value && std::is_integral<C>::value )
389 return ( (x>=0) && (y>=0) && (z>=0) && (x<ColumnsX) && (y<RowsY) && (z<SlicesZ) );
390 return ( (floor(x)>=0) && (floor(y)>=0) && (floor(z)>=0) && (ceil(x)<ColumnsX) && (ceil(y)<RowsY) && (ceil(z)<SlicesZ) );
391 }
392
393 inline bool in_bounds(int64_t t) const
394 { return ( (t>=0) && (t<this->tsize()) ); }
395
396 template<class A, class B, class C, class D>
397 inline bool in_bounds(A x, B y, C z, D t) const
398 { return ( this->in_bounds(t) && this->in_bounds(x,y,z) ); }
399
400 inline bool in_bounds(int64_t x, int64_t y, int64_t z, int64_t t, int64_t d5, int64_t d6=0, int64_t d7=0) const {
401 return (this->in_bounds(t) && this->in_bounds(x,y,z) && (d5>=0) && (d6>=0) && (d7>=0) && (d5<dim5) && (d6<dim6) && (d7<dim7));}
402
403 bool in_extraslice_bounds(float x, float y, float z) const
404 {
405 int64_t ix=((int64_t) floor(x));
406 int64_t iy=((int64_t) floor(y));
407 int64_t iz=((int64_t) floor(z));
408 return((ix>=-1) && (iy>=-1) && (iz>=-1) && (ix<ColumnsX) && (iy<RowsY) && (iz<SlicesZ));
409 }
410 inline bool valid(int64_t x, int64_t y, int64_t z) const
411 {
412 return((ep_valid[0] || (x>=0 && x<ColumnsX)) && (ep_valid[1] || (y>=0 && y<RowsY)) && (ep_valid[2] || (z>=0 && z<SlicesZ)));
413 }
414 bool valid(float x, float y, float z, double tol=1e-8) const
415 {
416 // int ix=((int) floor(x));
417 // int iy=((int) floor(y));
418 // int iz=((int) floor(z));
419 return((ep_valid[0] || (x+tol >= 0.0 && x <= ColumnsX-1+tol)) &&
420 (ep_valid[1] || (y+tol >= 0.0 && y <= RowsY-1+tol)) &&
421 (ep_valid[2] || (z+tol >= 0.0 && z <= SlicesZ-1+tol)));
422 // return((ep_valid[0] || (ix>=0 && (ix+1)<ColumnsX)) && (ep_valid[1] || (iy>=0 && (iy+1)<RowsY)) && (ep_valid[2] || (iz>=0 && (iz+1)<SlicesZ)));
423 }
424
425 inline T& operator()(int64_t x, int64_t y, int64_t z) {
426 this->invalidateSplines();
427 if (in_bounds(x,y,z)) return *(basicptr(x,y,z));
428 else return const_cast<T& > (extrapolate(x,y,z));
429 }
430
431 inline const T& operator()(int64_t x, int64_t y, int64_t z) const {
432 if (in_bounds(x,y,z)) return *(basicptr(x,y,z));
433 else return extrapolate(x,y,z);
434 }
435
436 inline T& operator()(int64_t x, int64_t y, int64_t z, int64_t t) {
437 this->invalidateSplines();
438 if (!in_bounds(t)) imthrow("Out of Bounds (time index)",5);
439 else if (!in_bounds(x,y,z)) return const_cast<T& > (operator[](t).extrapolate(x,y,z));
440 return *(Data + ((zsize()*t + z)*ysize() + y)*xsize() + x);
441 }
442
443 inline const T& operator()(int64_t x, int64_t y, int64_t z, int64_t t) const
444 { if (!in_bounds(t)) imthrow("Out of Bounds (time index)",5);
445 else if (!in_bounds(x,y,z)) return const_cast<T& > (operator[](t).extrapolate(x,y,z));
446 return *(Data + ((zsize()*t + z)*ysize() + y)*xsize() + x); }
447
448 inline T& operator()(int64_t x, int64_t y, int64_t z, int64_t t, int64_t d5, int64_t d6=0, int64_t d7=0) {
449 this->invalidateSplines();
450 if (!in_bounds(x,y,z,t,d5,d6,d7)) imthrow("Out of Bounds (7D index)",5);
451 return *(Data + ((((((size6()*d7+d6)*size5()+d5)*tsize()+t)*zsize()*t+z)*ysize()+y)*xsize()+ x));
452 }
453
454 inline const T& operator()(int64_t x, int64_t y, int64_t z, int64_t t, int64_t d5, int64_t d6=0, int64_t d7=0) const
455 { if (!in_bounds(x,y,z,t,d5,d6,d7)) imthrow("Out of Bounds (7D index)",5);
456 return *(Data + ((((((size6()*d7+d6)*size5()+d5)*tsize()+t)*zsize()*t+z)*ysize()+y)*xsize()+ x)); }
457
458 inline T splineCoef(int64_t x, int64_t y, int64_t z) const
459 {
460 if (!in_bounds(x,y,z)) imthrow("Out of Bounds",5);
461 return(splint.Coef(static_cast<unsigned int>(x),static_cast<unsigned int>(y),static_cast<unsigned int>(z)));
462 }
463
464 float interpolate(float x, float y, float z) const;
465 float interpolate(float x, float y, float z, bool *ep) const;
466 float interp1partial(// Input
467 float x, float y, float z, // Co-ordinates to get value for
468 int dir, // Direction for partial, 0->x, 1->y, 2->z
469 // Output
470 float *pderiv // Derivative returned here
471 ) const;
472 float interp3partial(// Input
473 float x, float y, float z, // Co-ordinate to get value for
474 // Output
475 float *dfdx, float *dfdy, float *dfdz // Partials
476 ) const;
477
478 inline T& value(int64_t x, int64_t y, int64_t z)
479 { return *(basicptr(x,y,z)); }
480 inline const T& value(int64_t x, int64_t y, int64_t z) const
481 { return *(basicptr(x,y,z)); }
482
483 inline T& value(int64_t x, int64_t y, int64_t z, int64_t t)
484 { return *(Data + ((zsize()*t + z)*ysize() + y)*xsize() + x); }
485 inline const T& value(int64_t x, int64_t y, int64_t z, int64_t t) const
486 { return *(Data + ((zsize()*t + z)*ysize() + y)*xsize() + x); }
487 inline T& value(int64_t x, int64_t y, int64_t z, int64_t t, int64_t d5, int64_t d6=0, int64_t d7=0)
488 { return *(Data + ((((((size6()*d7+d6)*size5()+d5)*tsize()+t)*zsize()+z)*ysize()+y)*xsize()+ x)); }
489 inline const T& value(int64_t x, int64_t y, int64_t z, int64_t t, int64_t d5, int64_t d6=0, int64_t d7=0) const
490 { return *(Data + ((((((size6()*d7+d6)*size5()+d5)*tsize()+t)*zsize()+z)*ysize()+y)*xsize()+ x)); }
491
492 float interpolatevalue(float x, float y, float z) const;
493 NEWMAT::ColumnVector ExtractRow(int j, int k) const;
494 NEWMAT::ColumnVector ExtractColumn(int i, int k) const;
495 void SetRow(int j, int k, const NEWMAT::ColumnVector& row);
496 void SetColumn(int j, int k, const NEWMAT::ColumnVector& col);
497
498 // ROI and higher dimensional accesses
499 void addvolume(const volume<T>& source);
500 void deletevolume(int64_t t);
501 void clear(); // deletes all volumes
502 NEWMAT::ReturnMatrix voxelts(int64_t x, int64_t y, int64_t z) const;
503 void setvoxelts(const NEWMAT::ColumnVector& ts, int64_t x, int64_t y, int64_t z);
504 void copySubVolume(const int64_t t, volume<T>& newVolume);
505 ShadowVolume<T> operator[](const int64_t t);
506 const ShadowVolume<T> operator[](const int64_t t) const;
507 const volume<T> constSubVolume(const int64_t t) const;
508 void replaceSubVolume(const int64_t t, const volume<T>& newVolume);
509 void copyROI(const volume<T>& source,const int64_t minx, const int64_t miny, const int64_t minz, const int64_t mint, const int64_t maxx, const int64_t maxy, const int64_t maxz, const int64_t maxt);
510 volume<T> ROI(const int64_t minx, const int64_t miny, const int64_t minz, const int64_t mint, const int64_t maxx, const int64_t maxy, const int64_t maxz, const int64_t maxt) const;
511
512 // SECONDARY FUNCTIONS
513 void setextrapolationmethod(extrapolation extrapmethod) const { p_extrapmethod = extrapmethod; }
514 extrapolation getextrapolationmethod() const { return(p_extrapmethod); }
515
516 void setpadvalue(T padval) const { padvalue = padval; }
517 T getpadvalue() const { return padvalue; }
518 T backgroundval() const;
519
520 void defineuserextrapolation(T (*extrap)(
521 const volume<T>& , int64_t, int64_t, int64_t)) const;
522 void setinterpolationmethod(interpolation interpmethod) const;
523 interpolation getinterpolationmethod() const { return p_interpmethod; }
524 void setsplineorder(int order) const;
525 inline void invalidateSplines() const { splineuptodate = false; }
526 int getsplineorder() const { return(splineorder); }
527 void forcesplinecoefcalculation() const;
528 void setextrapolationvalidity(bool xv, bool yv, bool zv) const { ep_valid[0]=xv; ep_valid[1]=yv; ep_valid[2]=zv; }
529 std::vector<bool> getextrapolationvalidity() const { return(ep_valid); }
530 void defineuserinterpolation(float (*interp)(
531 const volume<T>& , float, float, float)) const;
532 void definekernelinterpolation(const NEWMAT::ColumnVector& kx,
533 const NEWMAT::ColumnVector& ky,
534 const NEWMAT::ColumnVector& kz,
535 int wx, int wy, int wz) const; // full-width
536 void definekernelinterpolation(const volume<T>& vol) const;
537 void definesincinterpolation(const std::string& sincwindowtype,
538 int w, int nstore=1201) const; // full-width
539 void definesincinterpolation(const std::string& sincwindowtype,
540 int wx, int wy, int wz, int nstore=1201) const;
541 // full-width
542
543 inline void getneighbours(int64_t x, int64_t y, int64_t z,
544 T &v000, T &v001, T &v010,
545 T &v011, T &v100, T &v101,
546 T &v110, T &v111) const;
547 inline void getneighbours(int64_t x, int64_t y, int64_t z,
548 T &v000, T &v010,
549 T &v100, T &v110) const;
550
551
552
553 // ARITHMETIC FUNCTIONS
554 T operator=(T val);
555 const volume<T>& operator+=(T val);
556 const volume<T>& operator-=(T val);
557 const volume<T>& operator*=(T val);
558 const volume<T>& operator/=(T val);
559 const volume<T>& operator+=(const volume<T>& source);
560 const volume<T>& operator-=(const volume<T>& source);
561 const volume<T>& operator*=(const volume<T>& source);
562 const volume<T>& operator/=(const volume<T>& source);
563
564 volume<T> operator+(T num) const;
565 volume<T> operator-(T num) const;
566 volume<T> operator*(T num) const;
567 volume<T> operator/(T num) const;
568 volume<T> operator+(const volume<T>& vol2) const;
569 volume<T> operator-(const volume<T>& vol2) const;
570 volume<T> operator*(const volume<T>& vol2) const;
571 volume<T> operator/(const volume<T>& vol2) const;
572
573 template <class S>
574 friend volume<S> operator+(S num, const volume<S>& vol);
575 template <class S>
576 friend volume<S> operator-(S num, const volume<S>& vol);
577 template <class S>
578 friend volume<S> operator*(S num, const volume<S>& vol);
579 template <class S>
580 friend volume<S> operator/(S num, const volume<S>& vol);
581 template <class S>
582 friend volume<S> operator-(const volume<S>& vol);
583
584 // Comparisons. These are used for "spatial" purposes
585 // so that if data is identical and all the "spatial
586 // fields of the header are identical then the volumes
587 // are considered identical.
588 template <class S>
589 friend bool operator==(const volume<S>& v1, const volume<S>& v2);
590 template <class S>
591 friend bool operator!=(const volume<S>& v1, const volume<S>& v2); // { return(!(v1==v2)); }
592
593 // GENERAL MANIPULATION
594
595 //binarise: Values within range are set to 1, values outside to 0
596 //if invertRange is true this behaviour is reversed
597 void binarise(const T lower = smallestPositiveT,
598 const T upper = std::numeric_limits<T>::max(),
599 const threshtype tt = inclusive,
600 const bool invertRange = false);
601
602 //threshold: Values outside of the range will be zeroed
603 //if invertRange is true values inside the range will be zeroed
604 void threshold( const T lower,
605 const T upper = std::numeric_limits<T>::max(),
606 const threshtype tt = inclusive,
607 const bool invertRange = false);
608
609 // valid entries for dims are +/- 1, 2, 3 (and for newx, etc they are x, -x, y, -y, z, -z)
610 void swapdimensions(int dim1, int dim2, int dim3, bool keepLRorder=false);
611 void swapdimensions(const std::string& newx, const std::string& newy, const std::string& newz, const bool keepLRorder=false);
612 NEWMAT::Matrix swapmat(int dim1, int dim2, int dim3) const;
613 NEWMAT::Matrix swapmat(const std::string& newx, const std::string& newy, const std::string& newz) const;
614
615
616 // CONVERSION FUNCTIONS
617 template <class S, class D> friend
618 void copyconvert(const volume<S>& source, volume<D>& dest, const bool copyData);
619 };
620
621template <typename T>
622const T volume<T>::smallestPositiveT = std::is_integral<T>::value ? 1 : std::numeric_limits<T>::min();
623
624template<class T>
625class ShadowVolume : public volume<T> {
626 public:
627 using volume<T>::operator=;
628 private:
629 void setinterpolationmethod(interpolation interp) const { imthrow("Called private shadow method",101); }
630 void setextrapolationmethod(extrapolation extrapmethod) const { imthrow("Called private shadow method",101); }
631 bool assigned;
632 ShadowVolume() : assigned(false) {};
633 virtual int initialize(int64_t xsize, int64_t ysize, int64_t zsize, int64_t tsize, int64_t d5, int64_t d6, int64_t d7, T *d, bool d_owner, int64_t nt);
634 public:
635 const volume<T>& equals(const volume<T>& source);
636 ShadowVolume(const ShadowVolume<T>& source) : assigned(false) { this->reinitialize(source,ALIAS); }
637 ShadowVolume(const volume<T>& source) : assigned(false) { this->reinitialize(source,ALIAS); }
638};
639
640template<class T>
641class baseIterator {
642protected:
643 T* ptr;
644 T* const ptrBegin;
645 T* const ptrEnd;
646 const volume<T>* const source;
647 baseIterator() : ptr(nullptr), ptrBegin(nullptr), ptrEnd(nullptr), source(nullptr) { };
648public:
649 baseIterator(const volume<T>& input) : ptr((T*)input.fbegin()), ptrBegin(ptr), ptrEnd((T*)input.fend()), source(&input) {}
650 baseIterator& operator++() { ++ptr; return *this; }
651 baseIterator operator++(int) {baseIterator temp(*this); operator++(); return temp;}
652 inline T& operator*() { return *ptr;}
653 inline bool isValid() { return ptr != ptrEnd; }
654 inline size_t offset() { return ptr - ptrBegin; }
655};
656
657template<class T>
658 class cyclicIterator : public baseIterator<T> {
659 private:
660 cyclicIterator();
661 public:
662 cyclicIterator(const volume<T>& input) : baseIterator<T>(input) {};
663 cyclicIterator& operator++() { if ( ++baseIterator<T>::ptr == baseIterator<T>::ptrEnd ) baseIterator<T>::ptr=baseIterator<T>::ptrBegin; return *this; }
664 cyclicIterator operator++(int) {cyclicIterator temp(*this); operator++(); return temp;}
665 };
666
667template<class T,class U>
668class maskedIterator : public baseIterator<T> {
669private:
670 const U* mPtr;
671 const U* firstM;
672 const U* const mPtrEnd;
673 T* firstP;
674 const volume<U>* const mask;
675 bool usingMask;
676 maskedIterator() : mPtr(nullptr) {};
677public:
678 maskedIterator(const volume<T>& data,const volume<U>& dataMask,const std::string& exceptionSource="");
679 maskedIterator& operator++();
680 maskedIterator operator++(int) {maskedIterator temp(*this); operator++(); return temp;}
681 inline bool operator!=(const T* oPtr) { return baseIterator<T>::ptr != oPtr;}
682 maskedIterator& begin();
683 const T* const end() { return baseIterator<T>::ptrEnd; }
684 };
685
686template<class T, class U>
687maskedIterator<T,U>& maskedIterator<T,U>::begin() {
688 baseIterator<T>::ptr=firstP;
689 mPtr=firstM;
690 return *this;
691}
692
693template<class T, class U>
694maskedIterator<T,U>::maskedIterator(const volume<T>& data,const volume<U>& dataMask, const std::string& exceptionSource) : baseIterator<T>(data), mPtrEnd(dataMask.fend()), mask(&dataMask)
695{
696 mPtr=mask->fbegin(); //Should be safe even for NULL mask
697 usingMask=( mask->totalElements() != 0 );
698
699 if ( usingMask && !samesize(data,dataMask,SUBSET) )
700 throw std::runtime_error(exceptionSource+"maskedIterator: mask and volume must be the same size");
701 if ( usingMask && *mPtr <= mask->maskThreshold() ) //forward data pointer to first valid mask voxel
702 operator++();
703 firstP=baseIterator<T>::ptr;
704 firstM=mPtr;
705}
706
707
708template<class T, class U>
709maskedIterator<T,U>& maskedIterator<T,U>::operator++() {
710 //Increment the data pointer if currently valid
711 while ( this->isValid() && ++baseIterator<T>::ptr && usingMask ) {
712 //Seek/Cycle to the next valid mask voxel
713 if ( ++mPtr == mPtrEnd )
714 mPtr=mask->fbegin();
715 if ( *mPtr > mask->maskThreshold() ) //Is the candidate masked?
716 break;
717 }
718 return *this;
719}
720
721template<class T>
722 class safeLazyIterator : public baseIterator<T> { //Safe: cannot dereference if invalid, Lazy: invalidate source lazy flag on dereference
723public:
724 bool isValid() { return baseIterator<T>::ptr != nullptr && baseIterator<T>::ptr != baseIterator<T>::source->fend(); }
725 safeLazyIterator(const volume<T>& input) : baseIterator<T>(input) {}
726 safeLazyIterator& operator++() { ++baseIterator<T>::ptr ; return *this; }
727 safeLazyIterator operator++(int) {safeLazyIterator temp(*this); operator++(); return temp;}
728 T& operator*() {
729 this->invalidateSplines();
730 if(isValid()) return *baseIterator<T>::ptr;
731 else throw std::runtime_error("Attempted to dereference an invalid safeLazyIterator");
732 }
733 };
734
735template<class T>
736 class lazyIterator : public baseIterator<T> {
737public:
738 lazyIterator(const volume<T>& input) : baseIterator<T>(input) {}
739 lazyIterator& operator++() { ++baseIterator<T>::ptr ; return *this; }
740 lazyIterator operator++(int) {lazyIterator temp(*this); operator++(); return temp;}
741 T& operator*() {
742 this->invalidateSplines();
743 return *baseIterator<T>::ptr;
744 } //Need to invalidate source
745};
746
747
748
749 // HELPER FUNCTIONS
750 template <class S, class D>
751 void convertbuffer(const S* source, D* dest, size_t len, float slope=1.0, float intercept=0.0);
752
753 template <class S1, class S2>
754 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2);
755 template <class S1, class S2>
756 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2, bool checkdim);
757 template <class S1, class S2>
758 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2, int ndims, bool checkdim=false);
759
760 template <class S1, class S2>
761 bool samedim(const volume<S1>& vol1, const volume<S2>& vol2, int ndims);
762
766
767 template <class T>
768 inline void volume<T>::getneighbours(int64_t x, int64_t y, int64_t z,
769 T &v000, T &v001, T &v010,
770 T &v011, T &v100, T &v101,
771 T &v110, T &v111) const {
772 T *ptr = basicptr(x,y,z);
773 v000 = *ptr;
774 ptr++;
775 v100 = *ptr;
776 ptr+= ColumnsX;
777 v110 = *ptr;
778 ptr--;
779 v010 = *ptr;
780 ptr += RowsY * ColumnsX;
781 v011 = *ptr;
782 ptr++;
783 v111 = *ptr;
784 ptr-= ColumnsX;
785 v101 = *ptr;
786 ptr--;
787 v001 = *ptr;
788 }
789
790
791 template <class T>
792 inline void volume<T>::getneighbours(int64_t x, int64_t y, int64_t z,
793 T &v000, T &v010,
794 T &v100, T &v110) const {
795 T *ptr = basicptr(x,y,z);
796 v000 = *ptr;
797 ptr++;
798 v100 = *ptr;
799 ptr+= ColumnsX;
800 v110 = *ptr;
801 ptr--;
802 v010 = *ptr;
803 }
804
805
806
807
808
812
813 template <class T>
814 int64_t no_mask_voxels(const volume<T>& mask)
815 {
816 int64_t n(0);
817 for(maskedIterator<T,T> it(mask,mask);it.isValid();++it)
818 n++;
819 return n;
820 }
821
822
823 template <class S, class D>
824 void convertbuffer(const S* source, D* dest, size_t len, float slope, float intercept)
825 {
826 if ( slope == 1.0 && intercept == 0.0)
827 for (const S* sptr=source; sptr<(source+len); sptr++)
828 *(dest++) = (D) *sptr;
829 else
830 for (const S* sptr=source; sptr<(source+len); sptr++)
831 *(dest++) = (D) ((*sptr) * slope + intercept);
832 }
833
834 template <class S1, class S2>
835 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2, int ndims, bool checkdim)
836 {
837 int n=ndims;
838
839 if (ndims==SUBSET) { n = MISCMATHS::Min(vol1.dimensionality(),vol2.dimensionality()); }
840 bool same(vol1.xsize()==vol2.xsize());
841 if ( ndims == SUBSET )
842 same &= ( vol1.dimensionality() >= vol2.dimensionality() );
843 if (n>=2) same = same && ( vol1.ysize()==vol2.ysize());
844 if (n>=3) same = same && ( vol1.zsize()==vol2.zsize());
845 if (n>=4) same = same && ( vol1.tsize()==vol2.tsize());
846 if (n>=5) same = same && ( vol1.size5()==vol2.size5());
847 if (n>=6) same = same && ( vol1.size6()==vol2.size6());
848 if (n>=7) same = same && ( vol1.size7()==vol2.size7());
849 if (checkdim)
850 same = same && samedim(vol1,vol2,n);
851 if ( !same ) {
852 std::cerr << "Volumes are not compatible sizes:" << std::endl;
853 std::cerr << "Volume1:\t" << vol1.xsize() << " "
854 << vol1.ysize() << " "
855 << vol1.zsize() << " "
856 << vol1.tsize() << " "
857 << vol1.size5() << " "
858 << vol1.size6() << " "
859 << vol1.size7() << " " << std::endl;
860 std::cerr << "Volume2:\t" << vol2.xsize() << " "
861 << vol2.ysize() << " "
862 << vol2.zsize() << " "
863 << vol2.tsize() << " "
864 << vol2.size5() << " "
865 << vol2.size6() << " "
866 << vol2.size7() << " " << std::endl;
867 }
868 return(same);
869 }
870
871 template <class S1, class S2>
872 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2, bool checkdim)
873 {
874 return samesize(vol1,vol2,7,checkdim);
875 }
876
877 template <class S1, class S2>
878 bool samesize(const volume<S1>& vol1, const volume<S2>& vol2)
879 {
880 return samesize(vol1,vol2,7,false);
881 }
882
883 //The definition of samedim is if the difference in pixdim[i] is no greater than 0.01% of mean pixdim[i]
884 //fabs( dim_image1 - dim_image2 ) / ( ( dim_image1 + dim_image2 ) / 2) <= 0.0001
885 //The test is refactored to cater for the edge case of both pixdims being 0
886 template <class S1, class S2>
887 bool samedim(const volume<S1>& v1, const volume<S2>& v2, int ndims)
888 {
889 bool same=true;
890 same = ( std::fabs( v1.xdim()-v2.xdim() ) <= 5e-5 * ( v1.xdim()+v2.xdim() ) );
891 if (ndims>=2)
892 same = same && ( std::fabs( v1.ydim()-v2.ydim() ) <= 5e-5 * ( v1.ydim()+v2.ydim() ) );
893 if (ndims>=3)
894 same = same && ( std::fabs( v1.zdim()-v2.zdim() ) <= 5e-5 * ( v1.zdim()+v2.zdim() ) );
895 if (ndims>=4)
896 same = same && ( std::fabs( v1.tdim()-v2.tdim() ) <= 5e-5 * ( v1.tdim()+v2.tdim() ) );
897 if (ndims>=5)
898 same = same && ( std::fabs( v1.pixdim5()-v2.pixdim5() ) <= 5e-5 * ( v1.pixdim5()+v2.pixdim5() ) );
899 if (ndims>=6)
900 same = same && ( std::fabs( v1.pixdim6()-v2.pixdim6() ) <= 5e-5 * ( v1.pixdim6()+v2.pixdim6() ) );
901 if (ndims>=7)
902 same = same && ( std::fabs( v1.pixdim7()-v2.pixdim7() ) <= 5e-5 * ( v1.pixdim7()+v2.pixdim7() ) );
903 return same;
904 }
905
906
907
911
912 template <class S, class D>
913 void copybasicproperties(const volume<S>& source, volume<D>& dest)
914 {
915 // set up properties
916 dest.Xdim = source.Xdim;
917 dest.Ydim = source.Ydim;
918 dest.Zdim = source.Zdim;
919 dest.p_TR = source.p_TR;
920 dest.pxdim5 = source.pxdim5;
921 dest.pxdim6 = source.pxdim6;
922 dest.pxdim7 = source.pxdim7;
923
924 dest.originalSizes = source.originalSizes;
925
926 dest.StandardSpaceCoordMat = source.StandardSpaceCoordMat;
927 dest.RigidBodyCoordMat = source.RigidBodyCoordMat;
928 dest.StandardSpaceTypeCode = source.StandardSpaceTypeCode;
929 dest.RigidBodyTypeCode = source.RigidBodyTypeCode;
930
931 dest.RadiologicalFile = source.RadiologicalFile;
932
933 dest.IntentCode = source.IntentCode;
934 dest.IntentParam1 = source.IntentParam1;
935 dest.IntentParam2 = source.IntentParam2;
936 dest.IntentParam3 = source.IntentParam3;
937
938 dest.SliceOrderingCode = source.SliceOrderingCode;
939
940 dest.interpkernel = source.interpkernel;
941 dest.p_interpmethod = source.p_interpmethod;
942 dest.p_extrapmethod = source.p_extrapmethod;
943
944 dest.padvalue = (D) source.padvalue;
945 dest.splineorder = source.splineorder;
946 dest.ep_valid = source.ep_valid;
947
948 dest.displayMaximum=source.displayMaximum;
949 dest.displayMinimum=source.displayMinimum;
950 dest.setAuxFile(source.getAuxFile());
951
952 dest.extensions=source.extensions;
953 }
954
955
956 template <class S, class D>
957 void copyconvert(const volume<S>& source, volume<D>& dest, const bool copyData=true)
958 {
959 // set up basic size and data storage
960 if ( dest.totalElements() != source.totalElements() )
961 dest.reinitialize(source.xsize(),source.ysize(),source.zsize(),source.tsize(),source.size5(),source.size6(),source.size7());
962 // set up properties (except lazy ones)
963 copybasicproperties(source,dest);
964 // now copy across the data
965 if ( copyData )
966 convertbuffer(source.Data, dest.Data, source.totalElements() );
967 }
968
969 template <class S>
970 volume<S> operator+(S num, const volume<S>& vol)
971 { return (vol + num); }
972
973 template <class S>
974 volume<S> operator-(S num, const volume<S>& vol)
975 {
976 volume<S> tmp = vol;
977 tmp=num;
978 tmp-=vol;
979 return tmp;
980 }
981
982 template <class S>
983 volume<S> operator*(S num, const volume<S>& vol)
984 { return (vol * num); }
985
986 template <class S>
987 volume<S> operator/(S num, const volume<S>& vol)
988 {
989 volume<S> tmp = vol;
990 tmp=num;
991 tmp/=vol;
992 return tmp;
993 }
994
995 template <class S>
996 volume<S> operator-(const volume<S>& vol)
997 {
998 return(vol * (static_cast<S>(-1)));
999 }
1000
1001 template <class S>
1002 bool operator==(const volume<S>& v1,
1003 const volume<S>& v2)
1004 {
1005 // Check relevant parts of header
1006 if (!samesize(v1,v2,true)) return(false);
1007 if (v1.sform_code() != v2.sform_code()) return(false);
1008 if (v1.sform_mat() != v2.sform_mat()) return(false);
1009 if (v1.qform_code() != v2.qform_code()) return(false);
1010 if (v1.qform_mat() != v2.qform_mat()) return(false);
1011 // Check data
1012 for (typename volume<S>::fast_const_iterator it1=v1.fbegin(), it_end=v1.fend(), it2=v2.fbegin(); it1 != it_end; ++it1, ++it2) {
1013 if ((*it1) != (*it2)) return(false);
1014 }
1015
1016 return(true);
1017 }
1018 template <class S>
1019 bool operator!=(const volume<S>& v1, const volume<S>& v2) {return(!(v1==v2));}
1020
1021 template < class S, class V >
1022 std::vector<double> calculateSums(const S& inputVolume, const V& mask );
1023 template < class T, class V>
1024 std::vector<T> calculateExtrema(const volume<T>& inputVolume, std::vector<int64_t>& coordinates, const volume<V>& mask );
1025
1026} // end namespace
1027
1028#endif
NEWIMAGE::ShadowVolume
Definition: newimage.h:625
NEWIMAGE::baseIterator
Definition: newimage.h:641
NEWIMAGE::cyclicIterator
Definition: newimage.h:658
NEWIMAGE::fileFormats
Definition: newimage.h:61
NEWIMAGE::lazyIterator
Definition: newimage.h:736
NEWIMAGE::maskedIterator
Definition: newimage.h:668
NEWIMAGE::safeLazyIterator
Definition: newimage.h:722
NEWIMAGE::volume
Definition: newimage.h:100
NEWIMAGE::offset
Definition: newimagefns.h:1940