newimage
 
Loading...
Searching...
No Matches
newimagefns.h
1/* newimagefns.h
2
3 Mark Jenkinson et al, FMRIB Image Analysis Group
4
5 Copyright (C) 2000-2006 University of Oxford */
6
7/* CCOPYRIGHT */
8
9// General image processing functions
10
11#if !defined(__newimagefns_h)
12#define __newimagefns_h
13
14#include <string>
15#include <cstdlib>
16#include <iostream>
17#include <fstream>
18#include <sstream>
19#include <queue>
20#include "armawrap/newmat.h"
21#include "miscmaths/miscmaths.h"
22#include "newimage.h"
23#include "complexvolume.h"
24#include "imfft.h"
25
26
27#ifndef MAX
28#define MAX(a,b) (((a)>(b))?(a):(b))
29#endif
30
31#ifndef MIN
32#define MIN(a,b) (((a)<(b))?(a):(b))
33#endif
34
35namespace NEWIMAGE {
36
37
38 // The following lines are ignored by the current SGI compiler
39 // (version egcs-2.91.57)
40 // A temporary fix of including the std:: in front of all abs() etc
41 // has been done for now
42 using std::abs;
43 using std::sqrt;
44
46
47 // BASIC IMAGE SUPPORT FUNCTIONS
48
49 // Complex volume support
50 volume<float> abs(const volume<float>& realvol,
51 const volume<float>& imagvol);
52
53 volume<float> phase(const volume<float>& realvol,
54 const volume<float>& imagvol);
55
56 volume<float> real(const volume<float>& absvol,
57 const volume<float>& phasevol);
58
59 volume<float> imag(const volume<float>& absvol,
60 const volume<float>& phasevol);
61
62 // Basic Arithmetic Operations
63 template <class T>
64 volume<T> abs(const volume<T>& vol);
65
66 template <class T, class S>
67 volume<T> divide(const volume<T>& numervol, const volume<T>& denomvol,
68 const volume<S>& mask);
69
70 template <class T, class S>
71 volume<T> mask_volume( const volume<T>& invol,const volume<S>& mask );
72
73 // General Mathematical Operations
74
75 template <class T>
76 void clamp(volume<T>& vol, T minval, T maxval);
77
78
79 template <class T>
80 volume<T> binarise(const volume<T>& vol, T lowerth, T upperth, threshtype tt=inclusive, bool invert=false);
81 template <class T>
82 volume<T> binarise(const volume<T>& vol, T thresh, bool invert=false);
83
84
85 template <class T>
86 volume<T> threshold(const volume<T>& vol, T lowerth, T upperth, threshtype tt=inclusive);
87 template <class T>
88 volume<T> threshold(const volume<T>& vol, T thresh);
89
90
91 template <class T>
92 volume<T> min(const volume<T>& v1, const volume<T>& v2);
93 template <class T>
94 volume<T> max(const volume<T>& v1, const volume<T>& v2);
95
96 volume<float> sqrt(const volume<char>& vol);
97 volume<float> sqrt(const volume<short>& vol);
98 volume<float> sqrt(const volume<int>& vol);
99 volume<float> sqrt(const volume<float>& vol);
100 volume<double> sqrt(const volume<double>& vol);
101 template <class T>
102 volume<float> sqrt_float(const volume<T>& vol);
103
104 template <class T>
105 volume<float> meanvol(const volume4D<T>& vol4);
106 template <class T>
107 volume<float> stddevvol(const volume4D<T>& vol4);
108 template <class T>
109 volume<float> variancevol(const volume4D<T>& vol4);
110 template <class T>
111 volume<float> sumvol(const volume4D<T>& vol4);
112 template <class T>
113 volume<float> sumsquaresvol(const volume4D<T>& vol4);
114 template <class T>
115 volume<float> dotproductvol(const volume4D<T>& vol4,
116 const NEWMAT::ColumnVector& vec);
117
118 template <class T>
119 void pad(const volume<T>& vol, volume<T>& paddedvol);
120 template <class T>
121 void pad(const volume<T>& vol, volume<T>& paddedvol,
122 int offsetx, int offsety, int offsetz);
123
124 // Considers each volume as a vector and returns the dotproduct
125 template <class T>
126 double dotproduct(const volume<T>& vol1,
127 const volume<T>& vol2);
128 template <class T, class S>
129 double dotproduct(const volume<T>& vol1,
130 const volume<T>& vol2,
131 const volume<S>& mask);
132 template <class T, class S>
133 double dotproduct(const volume<T>& vol1,
134 const volume<T>& vol2,
135 const volume<S> *mask);
136
137 // This is a bit of a special-needs funtion. If we consider vol1 and vol2
138 // as column vectors v1 and v2 then powerdotproduct returns (v1.^n1)' * (v2.^n2)
139 template <class T>
140 double powerdotproduct(const volume<T>& vol1,
141 unsigned int n1,
142 const volume<T>& vol2,
143 unsigned int n2);
144 template <class T, class S>
145 double powerdotproduct(const volume<T>& vol1,
146 unsigned int n1,
147 const volume<T>& vol2,
148 unsigned int n2,
149 const volume<S>& mask);
150 template <class T, class S>
151 double powerdotproduct(const volume<T>& vol1,
152 unsigned int n1,
153 const volume<T>& vol2,
154 unsigned int n2,
155 const volume<S> *mask);
156
157 // These are global functions that duplicate some of the member functions.
158 // The reason is that we want to be able to double template them in a
159 // convenient manner (e.g. use a <char> mask for <float> data).
160 template <class T>
161 double mean(const volume<T>& vol);
162 template <class T, class S>
163 double mean(const volume<T>& vol,
164 const volume<S>& mask);
165 template <class T, class S>
166 double mean(const volume<T>& vol,
167 const volume<S> *mask);
168
169 template <class T>
170 double sum(const volume<T>& vol);
171 template <class T, class S>
172 double sum(const volume<T>& vol,
173 const volume<S>& mask);
174 template <class T, class S>
175 double sum(const volume<T>& vol,
176 const volume<S> *mask);
177
178
180
181 // IMAGE PROCESSING ROUTINES
182
183 template <class T>
184 void affine_transform(const volume<T>& vin, volume<T>& vout,
185 const NEWMAT::Matrix& aff, float paddingsize=0.0);
186 template <class T>
187 void affine_transform(const volume<T>& vin, volume<T>& vout,
188 const NEWMAT::Matrix& aff, float paddingsize, bool set_backgnd);
189
190 template <class T>
191 void affine_transform(const volume<T>& vin, volume<T>& vout,
192 const NEWMAT::Matrix& aff, interpolation interptype,
193 float paddingsize=0.0);
194 template <class T>
195 volume<T> affine_transform_mask(const volume<T>& vin, const volume<T>& vout,
196 const NEWMAT::Matrix& aff, float padding=0.0);
197
198 template <class T>
199 void get_axis_orientations(const volume<T>& inp1,
200 int& icode, int& jcode, int& kcode);
201
202
203 // the following convolve function do not attempt to normalise the kernel
204 template <class T, class S>
205 volume<T> convolve(const volume<T>& source, const volume<S>& kernel);
206 template <class T, class S>
207 volume<T> convolve_separable(const volume<T>& source,
208 const NEWMAT::ColumnVector& kernelx,
209 const NEWMAT::ColumnVector& kernely,
210 const NEWMAT::ColumnVector& kernelz);
211
212 // the following convolve functions take in a mask and also renormalise
213 // the result according to the overlap of kernel and mask at each point
214 // NB: these functions should NOT be used with zero-sum kernels (eg.Laplacian)
215 template <class T, class S, class M>
216 volume<T> convolve(const volume<T>& source, const volume<S>& kernel,
217 const volume<M>& mask, bool ignoremask=false, bool renormalise=true);
218 template <class T, class M>
219 volume<T> convolve_separable(const volume<T>& source,
220 const NEWMAT::ColumnVector& kernelx,
221 const NEWMAT::ColumnVector& kernely,
222 const NEWMAT::ColumnVector& kernelz,
223 const volume<M>& mask, bool ignoremask=false, bool renormalise=true);
224
225 //This implements Professor Smith's SUSAN convolve algorithm, note the number of optional parameters
226 template <class T, class S>
227 volume<T> susan_convolve(volume<T> source, const volume<S>& kernel, const float sigmabsq, const bool use_median, int num_usan,volume<T>* usan_area = new volume<T>(1,1,1),volume<T> usan_vol1=volume<T>(1,1,1),const float sigmab1sq=0,volume<T> usan_vol2 = volume<T>(1,1,1),const float sigmab2sq=0);
228
229 template <class T, class S>
230 volume4D<T> generic_convolve(const volume4D<T>& source, const volume<S>& kernel, bool seperable=false, bool renormalise=true);
231 template <class T, class S>
232 volume<T> efficient_convolve(const volume<T>& vin, const volume<S>& vker);
233 template <class T, class S>
234 int insertpart(volume<T>& v1, const volume<S>& v2);
235 template <class T, class S, class U>
236 volume<S> extractpart(const volume<T>& v1, const volume<S>& v2, const volume<U>& kernel) ;
237 float fsllog2(float x);
238
239
240 template <class T>
241 volume4D<T> bandpass_temporal_filter(volume4D<T>& source,double hp_sigma, double lp_sigma);
242
243
244 template <class T, class S>
245 volume<T> morphfilter(const volume<T>& source, const volume<S>& kernel,
246 const std::string& filtertype);
247
248 template <class T>
249 volume<T> dilall(const volume<T>& im, volume<T>& mask);
250 template <class T>
251 volume<T> fill_holes(const volume<T>& im, int connectivity=6);
252
253 template <class T>
254 volume<T> isotropic_resample(const volume<T>& aniso, float scale);
255
256 template <class T>
257 volume<T> subsample_by_2(const volume<T>& refvol, bool centred=true);
258 template <class T>
259 int upsample_by_2(volume<T>& highresvol, const volume<T>& lowresvol,
260 bool centred=true);
261 template <class T>
262 volume<T> upsample_by_2(const volume<T>& lowresvol, bool centred=true);
263
264 // for all blur functions the size of blurring is in mm
265 void make_blur_mask(NEWMAT::ColumnVector& bmask, const float final_vox_dim,
266 const float init_vox_dim);
267 template <class T>
268 volume<T> blur(const volume<T>& source, const NEWMAT::ColumnVector& resel_size);
269 template <class T>
270 volume<T> blur(const volume<T>& source, float iso_resel_size);
271
272 /*template <class T>
273 volume<T> smooth(const volume<T>& source, float sigma_mm);
274 template <class T>
275 volume<T> smooth2D(const volume<T>& source, float sigma_mm, int nulldir=3);*/
276
277
278 NEWMAT::ColumnVector gaussian_kernel1D(float sigma, int radius);
279 volume<float> gaussian_kernel2D(float sigma, int radius);
280 volume<float> gaussian_kernel3D(float sigma, int radius);
281 volume<float> gaussian_kernel3D(float sigma,float xdim,float ydim,float zdim,float cutoff=4.0);
282 volume<float> spherical_kernel(float radius, float xdim, float ydim, float zdim);
283 volume<float> box_kernel(float length,float xdim,float ydim,float zdim); //mm dimensions
284 volume<float> box_kernel(int x,int y, int z); //voxel dimensions
285
286 void make_grad_masks(volume<float>& maskx, volume<float>& masky,
287 volume<float>& maskz);
288
289 template <class T>
290 volume<float> gradient(const volume<T>& source); // in voxel coords
291
292 template <class T>
293 void gradient(const volume<T>& source,volume<float>& grad); // in voxel coords
294
295 // separate left and right gradients (changes at voxel mid-point)
296 template <class T>
297 volume4D<float> lrxgrad(const volume<float>& im, const volume<T>& mask);
298 template <class T>
299 volume4D<float> lrygrad(const volume<float>& im, const volume<T>& mask);
300 template <class T>
301 volume4D<float> lrzgrad(const volume<float>& im, const volume<T>& mask);
302
303
304 template <class T>
305 volume<T> log_edge_detect(const volume<T>& source,
306 float sigma1, float sigma2,
307 int mode=0); //mode==0: 3D, mode==1: 2D x, mode==2: 2D y, mode==3: 2D z
308 template <class T>
309 volume<T> fixed_edge_detect(const volume<T>& source, float threshold,
310 bool twodimensional=false);
311 template <class T>
312 volume4D<T> edge_strengthen(const volume4D<T>& source);
313
314
315
316 template <class T>
317 volume<int> connected_components(const volume<T>& vol, NEWMAT::ColumnVector& clustersize, int numconnected=26);
318
319 template <class T>
320 volume<int> connected_components(const volume<T>& vol,
321 const volume<T>& mask,
322 bool (*binaryrelation)(T , T), NEWMAT::ColumnVector& clustersize);
323
324 template <class T>
325 volume<int> connected_components(const volume<T>& vol,
326 int numconnected=26);
327 template <class T>
328 volume<int> connected_components(const volume<T>& vol,
329 const volume<T>& mask,
330 bool (*binaryrelation)(T , T));
331
332 template <class T>
333 volume<float> distancemap(const volume<T>& binaryvol);
334 template <class T>
335 volume<float> distancemap(const volume<T>& binaryvol, const volume<T>& mask);
336
337 template <class T>
338 volume4D<float> sparseinterpolate(const volume4D<T>& sparsesamps,
339 const volume<T>& mask,
340 const std::string& interpmethod="general");
341 // can have "general" or "nearestneighbour" (or "nn") for interpmethod
342
343
344 template <class T>
345 NEWMAT::Matrix NewimageVox2NewimageVoxMatrix(const NEWMAT::Matrix& flirt_in2ref,
346 const volume<T>& invol, const volume<T>& refvol);
347
348
349
352 // DEBUG
353 template <class T>
354 void print_volume_info(const volume<T>& source, const std::string& name, std::ostream& output=std::cout)
355 {
356 output << name << ":: Size = (" << source.xsize() << ","
357 << source.ysize() << "," << source.zsize() << "," << source.tsize() << ")" << std::endl;
358 output << name << ":: Dims = (" << source.xdim() << ","
359 << source.ydim() << "," << source.zdim() << "," << source.tdim() << ")" << std::endl;
360 output << name << ":: ROI Size = (" << source.maxx() - source.minx() + 1
361 << "," << source.maxy() - source.miny() + 1
362 << "," << source.maxz() - source.minz() + 1
363 << "," << source.maxt() - source.mint() + 1 << ")" << std::endl;
364 output << name << ":: Minimum and maximum intensities are: "
365 << source.min() << " and " << source.max() << std::endl;
366 }
367
369 // BASIC IMAGE SUPPORT FUNCTIONS
371
372 template <class T>
373 void clamp(volume<T>& vol, T minval, T maxval)
374 {
375 if (maxval < minval) return;
376 for ( typename volume<T>::nonsafe_fast_iterator it=vol.nsfbegin(), itEnd=vol.nsfend(); it < itEnd; it++ ) {
377 if ( *it > maxval )
378 *it = maxval;
379 else if ( *it < minval )
380 *it = minval;
381 }
382 }
383
384 template <class T>
385 volume<T> abs(const volume<T>& vol)
386 {
387 volume<T> newvol(vol);
388 for ( typename volume<T>::nonsafe_fast_iterator it=newvol.nsfbegin(), itEnd=newvol.nsfend(); it < itEnd; it++ ) {
389 *it = (T)fabs( (double)*it );
390 }
391 return newvol;
392 }
393
394 template <class T>
395 volume<T> binarise(const volume<T>& vol, T lowerth, T upperth, threshtype tt, bool invert)
396 {
397 volume<T> newvol(vol);
398 newvol.binarise(lowerth,upperth,tt,invert);
399 return newvol;
400 }
401
402 template <class T>
403 volume<T> binarise(const volume<T>& vol, T lowthresh, bool invert)
404 {
405 return binarise(vol,lowthresh,vol.max(),inclusive, invert);
406 }
407
408 template <class T>
409 volume<T> threshold(const volume<T>& vol, T lowerth, T upperth, threshtype tt)
410 {
411 volume<T> newvol(vol);
412 newvol.threshold(lowerth,upperth,tt);
413 return newvol;
414 }
415
416 template <class T>
417 volume<T> threshold(const volume<T>& vol, T thresh)
418 {
419 return threshold(vol,thresh,vol.max(),inclusive);
420 }
421
422 template <class T>
423 volume<T> min(const volume<T>& v1, const volume<T>& v2)
424 {
425 if (!samesize(v1,v2)) {
426 imthrow("Must use volumes of same size in min(v1,v2)",3);
427 }
428 volume<T> newvol(v1, false);
429 typename volume<T>::nonsafe_fast_iterator it=newvol.nsfbegin();
430 for ( typename volume<T>::fast_const_iterator v1it=v1.fbegin(), v1itEnd=v1.fend(), v2it=v2.fbegin(); v1it < v1itEnd; it++, v1it++, v2it++ )
431 *it=MISCMATHS::Min( *v1it, *v2it );
432 return newvol;
433 }
434
435 template <class T>
436 volume<T> max(const volume<T>& v1, const volume<T>& v2)
437 {
438 if (!samesize(v1,v2)) {
439 imthrow("Must use volumes of same size in max(v1,v2)",3);
440 }
441 volume<T> newvol;
442 copyconvert(v1, newvol, false);
443 typename volume<T>::nonsafe_fast_iterator it=newvol.nsfbegin();
444 for ( typename volume<T>::fast_const_iterator v1it=v1.fbegin(), v1itEnd=v1.fend(), v2it=v2.fbegin(); v1it < v1itEnd; it++, v1it++, v2it++ )
445 *it=MISCMATHS::Max( *v1it, *v2it );
446 return newvol;
447 }
448
449 template <class T>
450 volume<float> sqrt_float(const volume<T>& vol)
451 {
452 volume<float> retvol;
453 copyconvert(vol, retvol, false);
454 typename volume<float>::nonsafe_fast_iterator it=retvol.nsfbegin();
455 for ( typename volume<T>::fast_const_iterator vit=vol.fbegin(), vend=vol.fend(); vit < vend; it++, vit++ ) {
456 if ( *vit > 0 )
457 *it=sqrt( (double) *vit);
458 else
459 *it=0;
460 }
461 return retvol;
462 }
463
464 template <class T>
465 volume<float> sumvol(const volume<T>& vol4) //4D only
466 {
467 if (vol4.tsize()<1) { volume<float> newvol; return newvol; }
468 volume<float> SumVol, dummy;
469 copyconvert(vol4.constSubVolume(0),SumVol);
470 for (int ctr= 1; ctr < vol4.tsize(); ctr++) {
471 copyconvert(vol4.constSubVolume(ctr),dummy);
472 SumVol += dummy;
473 }
474 return SumVol;
475 }
476
477
478 template <class T>
479 volume<float> meanvol(const volume<T>& vol4)//4D only
480 {
481 if (vol4.tsize()<1) { volume<float> newvol; return newvol; }
482 volume<float> MeanVol = sumvol(vol4) / (float)vol4.tsize();
483 return MeanVol;
484 }
485
486
487 template <class T>
488 volume<float> sumsquaresvol(const volume<T>& vol4) //4D only
489 {
490 if (vol4.tsize()<1) { volume<float> newvol; return newvol; }
491 volume<float> SumSq, dummy;
492 copyconvert( vol4.constSubVolume(0) * vol4.constSubVolum(0) , SumSq);
493 for (int ctr=1; ctr < vol4.tsize(); ctr++) {
494 copyconvert( vol4.constSubVolume(ctr) * vol4.constSubVolume(ctr), dummy);
495 SumSq += dummy;
496 }
497 return SumSq;
498 }
499
500 //This rewrite of variancevol gives the same output as doing the sumsquaresvol etc in double precision internally
501 template <class T>
502 volume<float> variancevol(const volume<T>& vol4) //4D only
503 {
504 volume<float> variance;
505 if (vol4.tsize()<1)
506 return variance;
507 volume<float> Mean = meanvol(vol4);
508 variance = Mean*0.0f;
509
510 for (int z=0; z<vol4.zsize(); z++)
511 for (int y=0; y<vol4.ysize(); y++)
512 for (int x=0; x<vol4.xsize(); x++) {
513 double total(0);
514 for (int t=0; t<vol4.tsize(); t++)
515 total+=pow(vol4(x,y,z,t)-Mean(x,y,z),2.0);
516 variance(x,y,z)=(float)total;
517 }
518 variance /= (float) (vol4.tsize()-1.0);
519 return variance;
520 }
521
522 template <class T>
523 volume<float> stddevvol(const volume<T>& vol4) //4D only
524 {
525 if (vol4.tsize()<1) { volume<float> newvol; return newvol; }
526 volume<float> StdVol = variancevol(vol4);
527 for (int z=0; z<StdVol.zsize(); z++)
528 for (int y=0; y<StdVol.ysize(); y++)
529 for (int x=0; x<StdVol.xsize(); x++)
530 StdVol(x,y,z) = sqrt(StdVol(x,y,z));
531
532 return StdVol;
533 }
534
535 template <class T>
536 volume<float> dotproductvol(const volume4D<T>& vol4, //4D only
537 const NEWMAT::ColumnVector& vec)
538 {
539 if (vol4.mint()<0) { volume<float> newvol; return newvol; }
540 if ( vol4.tsize() != vec.Nrows() )
541 {
542 std::cerr << "ERROR::Time series length differs from vector length in"
543 << " dotproductvol()" << std::endl;
544 volume<float> newvol; return newvol;
545 }
546 volume<float> vol4copy;
547 copyconvert(vol4[vol4.mint()],vol4copy);
548 volume<float> DotVol(vol4copy);
549 DotVol *= (float) vec(1);
550 for (int n=vol4.mint() + 1; n <= vol4.maxt(); n++) {
551 copyconvert(vol4[n],vol4copy);
552 DotVol += (vol4copy * (float) vec(1 + n - vol4.mint()));
553 }
554 return DotVol;
555 }
556
557
558 template <class T>
559 void pad(const volume<T>& vol, volume<T>& paddedvol)
560 {
561 // The default type of padding is central padding
562 if ( ( paddedvol.xsize() < vol.xsize() ) || ( paddedvol.ysize() < vol.ysize() ) || ( paddedvol.zsize() < vol.zsize() ) ) {
563 imthrow("Cannot pad when target volume is smaller than original",7);
564 }
565 int64_t offx = ( paddedvol.xsize() - vol.xsize() ) / 2;
566 int64_t offy = ( paddedvol.ysize() - vol.ysize() ) / 2;
567 int64_t offz = ( paddedvol.zsize() - vol.zsize() ) / 2;
568 pad(vol,paddedvol,offx,offy,offz);
569 }
570
571 template <class T>
572 double dotproduct(const volume<T>& vol1,
573 const volume<T>& vol2)
574 {
575 if (!samesize(vol1,vol2,true)) imthrow("dotproduct: Image dimension mismatch",99);
576
577 double rval = 0.0;
578 for (typename volume<T>::fast_const_iterator it1=vol1.fbegin(), it_end=vol1.fend(), it2=vol2.fbegin(); it1 != it_end; ++it1, ++it2) {
579 rval += static_cast<double>((*it1)*(*it2));
580 }
581 return(rval);
582 }
583
584 template <class T, class S>
585 double dotproduct(const volume<T>& vol1,
586 const volume<T>& vol2,
587 const volume<S>& mask)
588 {
589 if (!samesize(vol1,vol2,true) || !samesize(vol1,mask,true)) imthrow("dotproduct: Image dimension mismatch",99);
590
591 double rval = 0.0;
592 typename volume<S>::fast_const_iterator itm = mask.fbegin();
593 for (typename volume<T>::fast_const_iterator it1=vol1.fbegin(), it_end=vol1.fend(), it2=vol2.fbegin(); it1 != it_end; ++it1, ++it2, ++itm) {
594 if (*itm > 0.5) {
595 rval += static_cast<double>((*it1)*(*it2));
596 }
597 }
598 return(rval);
599 }
600
601 template <class T, class S>
602 double dotproduct(const volume<T>& vol1,
603 const volume<T>& vol2,
604 const volume<S> *mask)
605 {
606 if (mask) return(dotproduct(vol1,vol2,*mask));
607 else return(dotproduct(vol1,vol2));
608 }
609
610 template <class T>
611 double powerdotproduct(const volume<T>& vol1,
612 unsigned int n1,
613 const volume<T>& vol2,
614 unsigned int n2)
615 {
616 if (!samesize(vol1,vol2,true)) imthrow("powerdotproduct: Image dimension mismatch",99);
617
618 double rval = 0.0;
619 for (typename volume<T>::fast_const_iterator it1=vol1.fbegin(), it_end=vol1.fend(), it2=vol2.fbegin(); it1 != it_end; ++it1, ++it2) {
620 double val1 = 1.0;
621 for (unsigned int i=0; i<n1; i++) val1 *= *it1;
622 double val2 = 1.0;
623 for (unsigned int i=0; i<n2; i++) val2 *= *it2;
624 rval += static_cast<double>(val1*val2);
625 }
626 return(rval);
627 }
628
629 template <class T, class S>
630 double powerdotproduct(const volume<T>& vol1,
631 unsigned int n1,
632 const volume<T>& vol2,
633 unsigned int n2,
634 const volume<S>& mask)
635 {
636 if (!samesize(vol1,vol2,true) || !samesize(vol1,mask,true)) imthrow("powerdotproduct: Image dimension mismatch",99);
637
638 double rval = 0.0;
639 typename volume<S>::fast_const_iterator itm = mask.fbegin();
640 for (typename volume<T>::fast_const_iterator it1=vol1.fbegin(), it_end=vol1.fend(), it2=vol2.fbegin(); it1 != it_end; ++it1, ++it2, ++itm) {
641 if (*itm > 0.5) {
642 double val1 = 1.0;
643 for (unsigned int i=0; i<n1; i++) val1 *= *it1;
644 double val2 = 1.0;
645 for (unsigned int i=0; i<n2; i++) val2 *= *it2;
646 rval += static_cast<double>(val1*val2);
647 }
648 }
649 return(rval);
650 }
651
652 template <class T, class S>
653 double powerdotproduct(const volume<T>& vol1,
654 unsigned int n1,
655 const volume<T>& vol2,
656 unsigned int n2,
657 const volume<S> *mask)
658 {
659 if (mask) return(powerdotproduct(vol1,n1,vol2,n2,*mask));
660 else return(powerdotproduct(vol1,n1,vol2,n2));
661 }
662
663 template <class T>
664 double mean(const volume<T>& vol)
665 {
666 double rval = sum(vol);
667 rval /= static_cast<double>( vol.totalElements() );
668 return(rval);
669 }
670
671 template <class T, class S>
672 double mean(const volume<T>& vol,
673 const volume<S>& mask)
674 {
675 if (!samesize(vol,mask,true)) imthrow("mean: Image-Mask dimension mismatch",99);
676
677 double rval(0);
678 int64_t n(0);
679 typename volume<S>::fast_const_iterator itm = mask.fbegin();
680 for (typename volume<T>::fast_const_iterator it=vol.fbegin(), it_end=vol.fend(); it != it_end; ++it, ++itm) {
681 if (*itm > 0.5) {
682 n++;
683 rval += static_cast<double>(*it);
684 }
685 }
686 rval /= static_cast<double>(n);
687 return(rval);
688 }
689
690 template <class T, class S>
691 double mean(const volume<T>& vol,
692 const volume<S> *mask)
693 {
694 if (mask) return(mean(vol,*mask));
695 else return(mean(vol));
696 }
697
698 template <class T>
699 double sum(const volume<T>& vol)
700 {
701 double rval(0);
702 for (typename volume<T>::fast_const_iterator it=vol.fbegin(), it_end=vol.fend(); it != it_end; ++it) {
703 rval += static_cast<double>(*it);
704 }
705 return(rval);
706 }
707
708 template <class T, class S>
709 double sum(const volume<T>& vol,
710 const volume<S>& mask)
711 {
712 if (!samesize(vol,mask,true)) imthrow("sum: Image-Mask dimension mismatch",99);
713
714 double rval(0);
715 typename volume<S>::fast_const_iterator itm = mask.fbegin();
716 for (typename volume<T>::fast_const_iterator it=vol.fbegin(), it_end=vol.fend(); it != it_end; ++it, ++itm) {
717 if (*itm > 0.5) {
718 rval += static_cast<double>(*it);
719 }
720 }
721 return(rval);
722 }
723
724 template <class T, class S>
725 double sum(const volume<T>& vol,
726 const volume<S> *mask)
727 {
728 if (mask) return(sum(vol,*mask));
729 else return(sum(vol));
730 }
731
732 template <class T, class S>
733 volume<T> divide(const volume<T>& numervol, const volume<T>& denomvol,
734 const volume<S>& mask)
735 {
736 if ((!samesize(numervol,denomvol)) || (!samesize(mask,denomvol,SUBSET))) {
737 imthrow("Attempted to divide images of different sizes",3);
738 }
739 volume<T> resvol(numervol);
740
741 typename volume<S>::fast_const_iterator maskIt = mask.fbegin(), maskEnd=mask.fend();
742 typename volume<T>::nonsafe_fast_iterator it=resvol.nsfbegin();
743 for (typename volume<T>::fast_const_iterator itd=denomvol.fbegin(), itd_end=denomvol.fend(); itd != itd_end; ++it, ++itd, maskIt++) {
744 if ( *maskIt != 0 ) {
745 *it /= *itd;
746 } else {
747 *it = 0;
748 }
749 if ( maskIt == maskEnd )
750 maskIt=mask.fbegin();
751 }
752 return resvol;
753 }
754
755
756
757template <class T, class S>
758 volume<T> mask_volume( const volume<T>& invol, const volume<S>& mask)
759 {
760 if ( !samesize(invol,mask,SUBSET)) {
761 imthrow("Attempted to mask_volume with wrong sized mask",3);
762 }
763 volume<T> resvol(invol, false);
764
765 typename volume<S>::fast_const_iterator maskIt = mask.fbegin(), maskEnd=mask.fend();
766 typename volume<T>::nonsafe_fast_iterator itr=resvol.nsfbegin();
767 for (typename volume<T>::fast_const_iterator it=invol.fbegin(), it_end=invol.fend(); it != it_end; ++it, ++itr, maskIt++) {
768 if ( *maskIt != 0 ) {
769 *itr = *it;
770 } else {
771 *itr = 0;
772 }
773 if ( maskIt == maskEnd )
774 maskIt=mask.fbegin();
775 }
776 return resvol;
777 }
778
779
780 // AFFINE TRANSFORM
781 template <class T>
782 void raw_affine_transform(const volume<T>& vin, volume<T>& vout,
783 const NEWMAT::Matrix& aff);
784
785 template <class T>
786 void affine_transform_mask(const volume<T>& vin, volume<T>& vout,
787 const NEWMAT::Matrix& aff, float padding, const T padval);
788
789 template <class T>
790 volume<T> affine_transform_mask(const volume<T>& vin, const volume<T>& vout,
791 const NEWMAT::Matrix& aff, float padding)
792 {
793 volume<T> affmask;
794 affmask = vout;
795 affmask = (T) 1;
796 affine_transform_mask(vin,affmask,aff,padding,(T) 0);
797 return affmask;
798 }
799
800
801 template <class T>
802 void affine_transform(const volume<T>& vin, volume<T>& vout,
803 const NEWMAT::Matrix& aff, float paddingsize, bool set_backgnd)
804 {
805 T padval=0;
806 extrapolation oldex;
807 if (set_backgnd) {
808 padval = vin.getpadvalue();
809 oldex = vin.getextrapolationmethod();
810 vin.setpadvalue(vin.backgroundval());
811 vin.setextrapolationmethod(extraslice);
812 }
813
814
815 raw_affine_transform(vin,vout,aff);
816 // now mask the output to eliminate streaks formed by the sinc interp...
817 affine_transform_mask(vin,vout,aff,paddingsize,vin.getpadvalue());
818 if (set_backgnd) {
819 vin.setpadvalue(padval);
820 vin.setextrapolationmethod(oldex);
821 }
822 }
823
824 template <class T>
825 void affine_transform(const volume<T>& vin, volume<T>& vout,
826 const NEWMAT::Matrix& aff, float paddingsize)
827 {
828 affine_transform(vin,vout,aff,paddingsize,true);
829 }
830
831 template <class T>
832 void affine_transform(const volume<T>& vin, volume<T>& vout,
833 const NEWMAT::Matrix& aff, interpolation interptype,
834 float paddingsize)
835 {
836 interpolation oldinterp;
837 oldinterp = vin.getinterpolationmethod();
838 vin.setinterpolationmethod(interptype);
839 affine_transform(vin,vout,aff,paddingsize);
840 vin.setinterpolationmethod(oldinterp);
841 }
842
843
844
846 // CONVOLVE
847 template <class T, class S>
848 volume<T> susan_convolve(const volume<T> source, const volume<S>& kernel, const float sigmabsq, const bool use_median, int num_usan,volume<T>* usan_area,volume<T> usan_vol1,const float sigmab1sq,volume<T> usan_vol2,const float sigmab2sq)
849 //template <class T, class S, class U, class V, class W>
850 //volume<T> susan_convolve(const volume<T>& source, const volume<S>& kernel, const float sigmabsq, const bool use_median, int num_usan,volume<U>* usan_area = new volume<T>(1,1,1),const volume<V>& usan_vol1=volume<T>(1,1,1),const float sigmab1sq=0,const volume<W>& usan_vol2 = volume<T>(1,1,1),const float sigmab2sq=0)
851 //Note that the commented out declaration won't work with the optional arguements (since U,V,W need to be defined in call...). Code is provided for possible
852 //future improvments, as it is all usans are templated as input
853{
854//need to use a pointer for usan_area as creating a default parameter for a pass-by-reference gives
855//a "assignment to tempory memory" warning in gcc
856//default values for lut1 etc
857 if ( (( (kernel.maxz() - kernel.minz()) % 2)==1) ||
858 (( (kernel.maxy() - kernel.miny()) % 2)==1) ||
859 (( (kernel.maxx() - kernel.minx()) % 2)==1) )
860 std::cerr << "WARNING:: Off-centre convolution being performed as kernel has even dimensions" << std::endl;
861 if ((num_usan>=1 && !samesize(source,usan_vol1)) || (num_usan>=2 && !samesize(source,usan_vol2)) || (num_usan>=1 && !samesize(source,*usan_area)) )
862 {
863 std::cerr << "Warning: an external usan or output usan is not the same size as the source image, reverting to num_usans=0 mode" << std::endl;
864 num_usan=0;
865 }
866 int midx, midy, midz,lz,uz,lx,ux,ly,uy,lutsize=16384;
867 lz=source.minz();
868 uz=source.maxz();
869 lx=source.minx();
870 ux=source.maxx();
871 ly=source.miny();
872 uy=source.maxy();
873 volume<T> result(source);
874 midz=(kernel.maxz() - kernel.minz())/2;
875 midy=(kernel.maxy() - kernel.miny())/2;
876 midx=(kernel.maxx() - kernel.minx())/2;
877 //generate look up table
878 float range1=1,range2=1,range = (source.max() - source.min())/(float)lutsize;
879 float **lut=new float *[3];
880 for(int i=0;i<=num_usan;i++) lut[i]=new float[2*lutsize+1];
881 for(int i=0;i<=num_usan;i++) lut[i]+=lutsize;
882 for (int i=0;i<=lutsize;i++) lut[0][-i]=lut[0][i]= exp(-pow(i*range,2.0)/sigmabsq);
883 if (num_usan>=1)
884 {
885 range1= (usan_vol1.max() - usan_vol1.min())/(float)lutsize;
886 for (int i=0;i<=lutsize;i++) lut[1][-i]=lut[1][i]= exp(-pow(i*range1,2.0)/sigmab1sq);
887 }
888 if (num_usan>=2)
889 {
890 range2= (usan_vol2.max() - usan_vol2.min())/(float)lutsize;
891 for (int i=0;i<=lutsize;i++) lut[2][-i]=lut[2][i]= exp(-pow(i*range2,2.0)/sigmab2sq);
892 }
893 NEWMAT::ColumnVector mediankernel((kernel.zsize()>1)?26:8); // cube or square, minus central voxel
894 int medoffst=((kernel.zsize()>1)?1:0);
895 for (int z=lz; z<=uz; z++)
896 for (int y=ly; y<=uy; y++)
897 for (int x=lx; x<=ux; x++)
898 {
899 int xmin=x-midx,ymin=y-midy,zmin=z-midz;
900 int xmax=MIN(x+midx,ux),ymax=MIN(y+midy,uy),zmax=MIN(z+midz,uz);
901 float num=0, denom=0,center_val1=0,center_val2=0,factor;
902 float center_val=source.value(x,y,z);
903 if (num_usan>=1) center_val1=usan_vol1.value(x,y,z);
904 if (num_usan>=2) center_val2=usan_vol2.value(x,y,z);
905 for(int mz=MAX(zmin,lz); mz<=zmax; mz++)
906 for(int my=MAX(ymin,ly); my<=ymax; my++)
907 for(int mx=MAX(xmin,lx); mx<=xmax; mx++)
908 if ((factor=(float)kernel.value(mx-xmin,my-ymin,mz-zmin)))
909 {
910 if (num_usan==0) factor*= lut[0][(int)((source.value(mx,my,mz)-center_val)/range)];
911 else factor*= lut[1][(int)((usan_vol1.value(mx,my,mz)-center_val1)/range1)];
912 if (num_usan>=2) factor*=lut[2][(int)((usan_vol2.value(mx,my,mz)-center_val2)/range2)];
913 num+=source.value(mx,my,mz) * factor;
914 denom+=factor;
915 }
916 if (num_usan>=1) usan_area->value(x,y,z)=(T) denom;
917 if (use_median && denom<1.5)
918 {
919 int count=1;
920 for(int x2=MAX(x-1,lx);x2<=MIN(x+1,ux);x2++)
921 for(int y2=MAX(y-1,ly);y2<=MIN(y+1,uy);y2++)
922 for(int z2=MAX(z-medoffst,lz);z2<=MIN(z+medoffst,uz);z2++)
923 if ( (x2-x) || (y2-y) || (z2-z) ) mediankernel(count++)=source.value(x2,y2,z2);
924 NEWMAT::ColumnVector subkernel = mediankernel.SubMatrix(1,count-1,1,1);
925 SortAscending(subkernel);
926 result(x,y,z) = (T)((subkernel(count/2)+subkernel((count+1)/2))/2.0);
927 }
928 else result.value(x,y,z)=(T) (num/denom);
929 }
930 for(int i=0;i<=num_usan;i++) delete[] (lut[i]-lutsize);
931 delete[] lut;
932 return result;
933}
934
935 template <class T, class S>
936 volume<T> convolve(const volume<T>& source, const volume<S>& kernel)
937 {
938 extrapolation oldex = source.getextrapolationmethod();
939 int offset=0;
940 if ((oldex==boundsassert) || (oldex==boundsexception))
941 { source.setextrapolationmethod(constpad); }
942 volume<T> result(source);
943 if ( (kernel.zsize()%2==0) || (kernel.ysize()%2==0) || (kernel.xsize()%2==0) )
944 {
945 std::cerr << "WARNING:: Off-centre convolution being performed as kernel"
946 << " has even dimensions" << std::endl;
947 //offset=2;
948 //offset gives correction to convolve to match results with fft for even kernel
949 //not even kernel with normalise (e.g. -fmean in fslmaths) still has edge problems
950 }
951 int midx, midy, midz;
952 midz=(kernel.zsize()-1)/2 + offset;
953 midy=(kernel.ysize()-1)/2 + offset;
954 midx=(kernel.xsize()-1)/2 + offset;
955
956 float val;
957 for (int z=0; z<result.zsize(); z++)
958 for (int y=0; y<result.ysize(); y++)
959 for (int x=0; x<result.xsize(); x++) {
960 val=0.0;
961 for (int mz=0; mz<kernel.zsize(); mz++)
962 for (int my=0; my<kernel.ysize(); my++)
963 for (int mx=0; mx<kernel.xsize(); mx++) {
964 val+=source(x+mx-midx,y+my-midy,z+mz-midz) * kernel(mx,my,mz);
965 }
966 result(x,y,z)=(T) val;
967 }
968
969
970 source.setextrapolationmethod(oldex);
971 return result;
972 }
973
974 template <class T, class S, class M>
975 volume<T> convolve(const volume<T>& source, const volume<S>& kernel,
976 const volume<M>& mask, bool ignoremask, bool renormalise)
977 {
978 if (!ignoremask && !samesize(mask, source))
979 imthrow("convolve: mask and source are not the same size",10);
980
981 if ( (( (kernel.maxz() - kernel.minz()) % 2)==1) ||
982 (( (kernel.maxy() - kernel.miny()) % 2)==1) ||
983 (( (kernel.maxx() - kernel.minx()) % 2)==1) )
984 {
985 std::cerr << "WARNING:: Off-centre convolution being performed as kernel"
986 << " has even dimensions" << std::endl;
987 }
988 volume<T> result(source);
989 int midx, midy, midz;
990 int lz,uz,lx,ux,ly,uy;
991 lz=result.minz();
992 uz=result.maxz();
993 lx=result.minx();
994 ux=result.maxx();
995 ly=result.miny();
996 uy=result.maxy();
997 midz=(kernel.maxz() - kernel.minz())/2;
998 midy=(kernel.maxy() - kernel.miny())/2;
999 midx=(kernel.maxx() - kernel.minx())/2;
1000 for (int z=lz; z<=uz; z++)
1001 for (int y=ly; y<=uy; y++)
1002 for (int x=lx; x<=ux; x++)
1003 if (ignoremask || mask(x,y,z)>0.5)
1004 {
1005 float val(0),norm(0);
1006 int x3,y3,z3;
1007 x3=x-midx;
1008 y3=y-midy;
1009 z3=z-midz;
1010 for (int mz=kernel.minz(); mz<=kernel.maxz(); mz++)
1011 for (int my=kernel.miny(); my<=kernel.maxy(); my++)
1012 for (int mx=kernel.minx(); mx<=kernel.maxx(); mx++)
1013 {
1014 int x2,y2,z2;
1015 x2=x3+mx;
1016 y2=y3+my;
1017 z2=z3+mz;
1018 if ((ignoremask && (x2<=ux && x2>=lx && y2<=uy && y2>=ly && z2<=uz && z2>=lz)) || (!ignoremask && mask(x2,y2,z2)>0.5))
1019 {
1020 val+=source.value(x2,y2,z2) * kernel.value(mx,my,mz);
1021 norm+=kernel.value(mx,my,mz);
1022 }
1023 }
1024 if (renormalise && fabs(norm)>1e-12) result.value(x,y,z)=(T) (val/norm);
1025 else result.value(x,y,z)=(T) val;
1026 }
1027 return result;
1028 }
1029
1030 template <class T>
1031 volume<T> convolve_separable(const volume<T>& source,
1032 const NEWMAT::ColumnVector& kernelx,
1033 const NEWMAT::ColumnVector& kernely,
1034 const NEWMAT::ColumnVector& kernelz)
1035 {
1036 volume<T> result(source);
1037 volume<double> kerx(kernelx.Nrows(),1,1);
1038 volume<double> kery(1,kernely.Nrows(),1);
1039 volume<double> kerz(1,1,kernelz.Nrows());
1040 for (int n=1; n<=kernelx.Nrows(); n++) kerx.value(n-1,0,0) = kernelx(n);
1041 for (int n=1; n<=kernely.Nrows(); n++) kery.value(0,n-1,0) = kernely(n);
1042 for (int n=1; n<=kernelz.Nrows(); n++) kerz.value(0,0,n-1) = kernelz(n);
1043 result = convolve(result,kerx);
1044 result = convolve(result,kery);
1045 result = convolve(result,kerz);
1046 return result;
1047 }
1048
1049 template <class T, class M>
1050 volume<T> convolve_separable(const volume<T>& source,
1051 const NEWMAT::ColumnVector& kernelx,
1052 const NEWMAT::ColumnVector& kernely,
1053 const NEWMAT::ColumnVector& kernelz,
1054 const volume<M>& mask, bool ignoremask, bool renormalise)
1055 {
1056 volume<T> result(source);
1057 volume<double> kerx(kernelx.Nrows(),1,1);
1058 volume<double> kery(1,kernely.Nrows(),1);
1059 volume<double> kerz(1,1,kernelz.Nrows());
1060 for (int n=1; n<=kernelx.Nrows(); n++) kerx.value(n-1,0,0) = kernelx(n);
1061 for (int n=1; n<=kernely.Nrows(); n++) kery.value(0,n-1,0) = kernely(n);
1062 for (int n=1; n<=kernelz.Nrows(); n++) kerz.value(0,0,n-1) = kernelz(n);
1063 result = convolve(result,kerx,mask,ignoremask,renormalise);
1064 result = convolve(result,kery,mask,ignoremask,renormalise);
1065 result = convolve(result,kerz,mask,ignoremask,renormalise);
1066 return result;
1067 }
1068
1070 // GENERAL DILATION INCLUDING MEDIAN FILTERING
1071
1072 template <class T, class S>
1073 volume<T> morphfilter(const volume<T>& source, const volume<S>& kernel,
1074 const std::string& filtertype)
1075 {
1076 // implements a whole range of filtering, set by the string filtertype:
1077 // dilateM (mean), dilateD (mode), median, max or dilate, min or erode,
1078 // erodeS (set to zero if any neighbour is zero)
1079 extrapolation oldex = source.getextrapolationmethod();
1080 if ((oldex==boundsassert) || (oldex==boundsexception))
1081 { source.setextrapolationmethod(constpad); }
1082 volume<T> result(source);
1083 result = 0;
1084
1085 int nker;
1086 {
1087 volume<S> dummy(kernel);
1088 dummy.binarise((S)0.5); //new cast to avoid warnings for int-type templates when compiling
1089 nker = (int) dummy.sum();
1090 }
1091
1092 int midx, midy, midz;
1093 if ( (( (kernel.maxz() - kernel.minz()) % 2)==1) ||
1094 (( (kernel.maxy() - kernel.miny()) % 2)==1) ||
1095 (( (kernel.maxx() - kernel.minx()) % 2)==1) )
1096 {
1097 std::cerr << "WARNING:: Off-centre morphfilter being performed as kernel"
1098 << " has even dimensions" << std::endl;
1099 }
1100 midz=(kernel.maxz() - kernel.minz())/2;
1101 midy=(kernel.maxy() - kernel.miny())/2;
1102 midx=(kernel.maxx() - kernel.minx())/2;
1103 int count=1;
1104
1105 for (int z=result.minz(); z<=result.maxz(); z++) {
1106 for (int y=result.miny(); y<=result.maxy(); y++) {
1107 for (int x=result.minx(); x<=result.maxx(); x++) {
1108 NEWMAT::ColumnVector vals(nker);
1109 count=1;
1110 for (int mz=MISCMATHS::Max(kernel.minz(),result.minz()-z+midz);
1111 mz<=MISCMATHS::Min(kernel.maxz(),result.maxz()-z+midz); mz++) {
1112 for (int my=MISCMATHS::Max(kernel.miny(),result.miny()-y+midy);
1113 my<=MISCMATHS::Min(kernel.maxy(),result.maxy()-y+midy); my++) {
1114 for (int mx=MISCMATHS::Max(kernel.minx(),result.minx()-x+midx);
1115 mx<=MISCMATHS::Min(kernel.maxx(),result.maxx()-x+midx); mx++) {
1116 if (kernel(mx,my,mz)>0.5) {
1117 if ((filtertype!="dilateM" && filtertype!="dilateD") || source(x+mx-midx,y+my-midy,z+mz-midz)) vals(count++) = source(x+mx-midx,y+my-midy,z+mz-midz);
1118 }
1119 }
1120 }
1121 }
1122 if (count>1) {
1123 NEWMAT::ColumnVector littlevals;
1124 littlevals = vals.SubMatrix(1,count-1,1,1);
1125 if (filtertype=="median") {
1126 SortAscending(littlevals); //count/2 works for odd kernel (even count) count+1 gives edge compatibility
1127 result(x,y,z) = (T)littlevals(MISCMATHS::Max(1,(count+1)/2)); //with steves IP, otherwise gives the IP median-1 element
1128 } else if ((filtertype=="max") || (filtertype=="dilate")) result(x,y,z) = (T)littlevals.Maximum();
1129 else if ((filtertype=="min") || (filtertype=="erode")) result(x,y,z) = (T)littlevals.Minimum();
1130 else if (filtertype=="erodeS") {if (source(x,y,z)!=0 && littlevals.Minimum()==0) result(x,y,z) = 0; else result(x,y,z) = source(x,y,z);}
1131 else if (filtertype=="dilateM") {if (source(x,y,z)==0) result(x,y,z) = (T)(littlevals.Sum()/--count); else result(x,y,z) = source(x,y,z);}
1132 else if (filtertype=="dilateD") {if (source(x,y,z)==0){
1133 SortDescending(littlevals);
1134 double max=littlevals(1);
1135 int maxn=1;
1136 double current=littlevals(1);
1137 int currentn=1;
1138 for(int i=2;i<count;i++)
1139 {
1140 if (littlevals(i)==current) currentn++;
1141 else
1142 {
1143 current=littlevals(i);
1144 if (currentn>maxn)
1145 {
1146 max=littlevals(i-1);
1147 maxn=currentn;
1148 }
1149 currentn=1;
1150 }
1151 }
1152 result(x,y,z) = (T)max;} else result(x,y,z) = source(x,y,z);
1153 }
1154 else imthrow("morphfilter:: Filter type " + filtertype + "unsupported",7);
1155 } else result(x,y,z) = source(x,y,z); // THE DEFAULT CASE (leave alone)
1156 }
1157 }
1158 }
1159 source.setextrapolationmethod(oldex);
1160 return result;
1161 }
1162
1163 template <class T>
1164 double dilateval(const volume<T>& im, const volume<T>& mask, int x, int y, int z)
1165 {
1166 double sum=0;
1167 int n=0;
1168 for (int zz=MISCMATHS::Max(0,z-1); zz<=MISCMATHS::Min(z+1,im.maxz()); zz++) {
1169 for (int yy=MISCMATHS::Max(0,y-1); yy<=MISCMATHS::Min(y+1,im.maxy()); yy++) {
1170 for (int xx=MISCMATHS::Max(0,x-1); xx<=MISCMATHS::Min(x+1,im.maxx()); xx++) {
1171 if ((mask(xx,yy,zz)>(T) 0.5) && !((xx==x) && (yy==y) && (zz==z))) {
1172 sum += im(xx,yy,zz);
1173 n++;
1174 }
1175 }
1176 }
1177 }
1178 return sum/(MISCMATHS::Max(n,1));
1179 }
1180
1181
1182 template <class T>
1183 bool ext_edge(const volume<T>& mask, int x, int y, int z)
1184 {
1185 if (mask(x,y,z)>(T) 0.5) { return false; }
1186 for (int zz=MISCMATHS::Max(0,z-1); zz<=MISCMATHS::Min(z+1,mask.maxz()); zz++) {
1187 for (int yy=MISCMATHS::Max(0,y-1); yy<=MISCMATHS::Min(y+1,mask.maxy()); yy++) {
1188 for (int xx=MISCMATHS::Max(0,x-1); xx<=MISCMATHS::Min(x+1,mask.maxx()); xx++) {
1189 if ((mask(xx,yy,zz)>(T) 0.5) && !((xx==x) && (yy==y) && (zz==z))) {
1190 return true;
1191 }
1192 }
1193 }
1194 }
1195 return false;
1196 }
1197
1198 struct vec3_temp { int x; int y; int z; };
1199
1200
1201 template <class T>
1202 volume<T> dilall(const volume<T>& input, volume<T>& mask)
1203 {
1204 volume<T> im(input);
1205 if (!samesize(input,mask)) { std::cerr << "ERROR::dilall::image are not the same size" << std::endl; }
1206 std::deque<vec3_temp> ptlist;
1207 vec3_temp v, newv;
1208 // initial pass
1209 for (int z=0; z<=im.maxz(); z++) {
1210 for (int y=0; y<=im.maxy(); y++) {
1211 for (int x=0; x<=im.maxx(); x++) {
1212 if (ext_edge(mask,x,y,z)) {
1213 v.x=x; v.y=y; v.z=z;
1214 ptlist.push_front(v);
1215 }
1216 }
1217 }
1218 }
1219 while (!ptlist.empty()) {
1220 v = ptlist.back();
1221 ptlist.pop_back();
1222 if (mask(v.x,v.y,v.z)<=(T) 0.5) { // only do it if the voxel is still unset
1223 im(v.x,v.y,v.z)=dilateval(im,mask,v.x,v.y,v.z);
1224 mask(v.x,v.y,v.z)=(T)1;
1225 // check neighbours and add them to the list if necessary
1226 for (int zz=MISCMATHS::Max(0,v.z-1); zz<=MISCMATHS::Min(v.z+1,im.maxz()); zz++) {
1227 for (int yy=MISCMATHS::Max(0,v.y-1); yy<=MISCMATHS::Min(v.y+1,im.maxy()); yy++) {
1228 for (int xx=MISCMATHS::Max(0,v.x-1); xx<=MISCMATHS::Min(v.x+1,im.maxx()); xx++) {
1229 if (ext_edge(mask,xx,yy,zz)) { newv.x=xx; newv.y=yy; newv.z=zz; ptlist.push_front(newv); }
1230 }
1231 }
1232 }
1233 }
1234 }
1235 return im;
1236 }
1237
1238
1239 template <class T>
1240 volume<T> fill_holes(const volume<T>& im, int connectivity)
1241 {
1242 volume<int> mask;
1243 std::set<int> edgelabs;
1244 mask = connected_components((T)1-im,connectivity);
1245 // go over edge of FOV and for each each non-zero value there store this
1246 // The loops below go over the yz, xz and xy faces of the FOV and double count edges and corners, but this doesn't matter here
1247 for (int z=0; z<=im.maxz(); z++) {
1248 for (int y=0; y<=im.maxy(); y++) {
1249 if (mask(0,y,z)!=0) edgelabs.insert(mask(0,y,z));
1250 if (mask(im.maxx(),y,z)!=0) edgelabs.insert(mask(im.maxx(),y,z));
1251 }
1252 }
1253 for (int z=0; z<=im.maxz(); z++) {
1254 for (int x=0; x<=im.maxx(); x++) {
1255 if (mask(x,0,z)!=0) edgelabs.insert(mask(x,0,z));
1256 if (mask(x,im.maxy(),z)!=0) edgelabs.insert(mask(x,im.maxy(),z));
1257 }
1258 }
1259 for (int y=0; y<=im.maxy(); y++) {
1260 for (int x=0; x<=im.maxx(); x++) {
1261 if (mask(x,y,0)!=0) edgelabs.insert(mask(x,y,0));
1262 if (mask(x,y,im.maxz())!=0) edgelabs.insert(mask(x,y,im.maxz()));
1263 }
1264 }
1265
1266 // zero any voxel containing a detected edge label value
1267 while (!edgelabs.empty()) {
1268 int val=*(edgelabs.begin());
1269 for (int z=0; z<=im.maxz(); z++) {
1270 for (int y=0; y<=im.maxy(); y++) {
1271 for (int x=0; x<=im.maxx(); x++) {
1272 if (mask(x,y,z)==val) mask(x,y,z)=0;
1273 }
1274 }
1275 }
1276 edgelabs.erase(val);
1277 }
1278 volume<T> imcopy(im);
1279 for (int z=0; z<=im.maxz(); z++)
1280 for (int y=0; y<=im.maxy(); y++)
1281 for (int x=0; x<=im.maxx(); x++)
1282 if (mask(x,y,z)>(T)0.5) imcopy(x,y,z)=1;
1283 return imcopy;
1284 }
1285
1286
1288 // RESAMPLE
1289
1290 template <class T>
1291 volume<T> upsample_by_2(const volume<T>& lowresvol, bool centred)
1292 {
1293 volume<T> res;
1294 upsample_by_2(res,lowresvol,centred);
1295 return res;
1296 }
1297
1298
1299
1301 // BLURRING
1302 template <class T>
1303 volume<T> blur(const volume<T>& source, const NEWMAT::ColumnVector& resel_size)
1304 {
1305 NEWMAT::ColumnVector bmaskx, bmasky, bmaskz;
1306 make_blur_mask(bmaskx,resel_size(1),source.xdim());
1307 make_blur_mask(bmasky,resel_size(2),source.ydim());
1308 make_blur_mask(bmaskz,resel_size(3),source.zdim());
1309 return convolve_separable(source,bmaskx,bmasky,bmaskz);
1310 }
1311
1312
1313 template <class T>
1314 volume<T> blur(const volume<T>& source, float iso_resel_size)
1315 {
1316 NEWMAT::ColumnVector resel_size(3);
1317 resel_size = iso_resel_size;
1318 return blur(source,resel_size);
1319 }
1320
1321 template <class T>
1322 volume<T> smooth(const volume<T>& source, float sigma_mm)
1323 {
1324 float sigmax, sigmay, sigmaz;
1325 sigmax = sigma_mm/source.xdim();
1326 sigmay = sigma_mm/source.ydim();
1327 sigmaz = sigma_mm/source.zdim();
1328 int nx=((int) (sigmax-0.001))*2 + 3;
1329 int ny=((int) (sigmay-0.001))*2 + 3;
1330 int nz=((int) (sigmaz-0.001))*2 + 3;
1331 NEWMAT::ColumnVector kernelx, kernely, kernelz;
1332 kernelx = gaussian_kernel1D(sigmax,nx);
1333 kernely = gaussian_kernel1D(sigmay,ny);
1334 kernelz = gaussian_kernel1D(sigmaz,nz);
1335 return convolve_separable(source,kernelx,kernely,kernelz);
1336 }
1337
1338
1339
1340
1341 template <class T>
1342 volume<T> smooth2D(const volume<T>& source, float sigma_mm, int nulldir=3)
1343 {
1344 float sigmax, sigmay, sigmaz;
1345 sigmax = sigma_mm/source.xdim();
1346 sigmay = sigma_mm/source.ydim();
1347 sigmaz = sigma_mm/source.zdim();
1348 int nx=((int) (sigmax-0.001))*2 + 3;
1349 int ny=((int) (sigmay-0.001))*2 + 3;
1350 int nz=((int) (sigmaz-0.001))*2 + 3;
1351 NEWMAT::ColumnVector kernelx, kernely, kernelz, nullker(1);
1352 kernelx = gaussian_kernel1D(sigmax,nx);
1353 kernely = gaussian_kernel1D(sigmay,ny);
1354 kernelz = gaussian_kernel1D(sigmaz,nz);
1355 nullker = 1;
1356 if (nulldir==1) {
1357 return convolve_separable(source,nullker,kernely,kernelz);
1358 } else if (nulldir==2) {
1359 return convolve_separable(source,kernelx,nullker,kernelz);
1360 } else {
1361 // Smoothing in the x-y plane is the default!
1362 return convolve_separable(source,kernelx,kernely,nullker);
1363 }
1364 }
1365
1366template <class T, class S>
1367 int insertpart(volume<T>& v1, const volume<S>& v2) //N.B. This has superficial similarities
1368{ //to copyconvert, but is in fact quite
1369 for (int z=v2.minz(); z<=v2.maxz(); z++) { //different...
1370 for (int y=v2.miny(); y<=v2.maxy(); y++) {
1371 for (int x=v2.minx(); x<=v2.maxx(); x++) {
1372 v1(x,y,z)=(T) v2(x,y,z);
1373 }
1374 }
1375 }
1376 return 0;
1377}
1378
1379
1380 template <class T, class S,class U>
1381volume<S> extractpart(const volume<T>& v1, const volume<S>& v2, const volume<U>& kernel)
1382{
1383 volume<S> vout=v2;
1384 vout = (S) 0.0;
1385 int kxoff = (kernel.xsize()-1)/2;
1386 int kyoff = (kernel.ysize()-1)/2;
1387 int kzoff = (kernel.zsize()-1)/2;
1388 for (int z=v2.minz(); z<=v2.maxz(); z++) {
1389 for (int y=v2.miny(); y<=v2.maxy(); y++) {
1390 for (int x=v2.minx(); x<=v2.maxx(); x++) {
1391 vout(x,y,z)=(S) v1(x+kxoff,y+kyoff,z+kzoff);
1392 }
1393 }
1394 }
1395 return vout;
1396}
1397
1399// Efficient FFT-based convolve
1400template <class T, class S>
1401volume<T> efficient_convolve(const volume<T>& vin, const volume<S>& vker)
1402{
1403 bool usefft=true;
1404 // estimate calculation time for the two methods and pick the best
1405 float offt = 2 * vin.nvoxels() * fsllog2(2 * vin.nvoxels());
1406 float osum = (float)vin.nvoxels() * (float)vker.nvoxels();
1407 // float cast to avoud overflow for large int multiplication
1408 float scalefactor = 44; // relative unit operation cost for fft vs sum
1409 usefft = (osum > offt * scalefactor);
1410 //cout << usefft << endl;
1411 if (usefft) {
1412 int sx = MISCMATHS::Max(vin.xsize(),vker.xsize())*2;
1413 int sy = MISCMATHS::Max(vin.ysize(),vker.ysize())*2;
1414 int sz = MISCMATHS::Max(vin.zsize(),vker.zsize())*2;
1415 complexvolume vif, vkf;
1416 vif.re().reinitialize(sx,sy,sz);
1417 //vif.re().copyproperties(vin); Is this needed check with MJ...
1418 vif.re() = 0.0;
1419 vif.im() = vif.re();
1420 vkf = vif;
1421 insertpart(vif.re(),vin);
1422 insertpart(vkf.re(),vker);
1423 fft3(vif);
1424 fft3(vkf);
1425 vif *= vkf;
1426 ifft3(vif);
1427 return extractpart(vif.re(),vin,vker);
1428 } else return convolve(vin,vker);
1429}
1430
1431 template <class T, class S>
1432 volume4D<T> generic_convolve(const volume4D<T>& source, const volume<S>& kernel, bool seperable, bool renormalise)
1433 {
1434 volume4D<T> result(source);
1435 if (seperable)
1436 {
1437 volume<double> kerx(kernel.xsize(),1,1);
1438 volume<double> kery(1,kernel.ysize(),1);
1439 volume<double> kerz(1,1,kernel.zsize());
1440 for (int n=0; n<kernel.xsize(); n++) kerx.value(n,0,0) = kernel(n,kernel.ysize()/2,kernel.zsize()/2);
1441 for (int n=0; n<kernel.ysize(); n++) kery.value(0,n,0) = kernel(kernel.xsize()/2,n,kernel.zsize()/2);
1442 for (int n=0; n<kernel.zsize(); n++) kerz.value(0,0,n) = kernel(kernel.xsize()/2,kernel.ysize()/2,n);
1443 volume<T> mask(1,1,1);
1444 for (int t=source.mint(); t<=source.maxt(); t++) result.replaceSubVolume(t,convolve(result.constSubVolume(t),kerx,mask,true,renormalise));
1445 for (int t=source.mint(); t<=source.maxt(); t++) result.replaceSubVolume(t,convolve(result.constSubVolume(t),kery,mask,true,renormalise));
1446 for (int t=source.mint(); t<=source.maxt(); t++) result.replaceSubVolume(t,convolve(result.constSubVolume(t),kerz,mask,true,renormalise));
1447 }
1448 else
1449 {
1450 volume<S> norm_kernel(kernel);
1451 if (kernel.sum()) norm_kernel/=kernel.sum();
1452 for (int t=source.mint(); t<=source.maxt(); t++) result.replaceSubVolume(t,efficient_convolve(source.constSubVolume(t),norm_kernel));
1453 result.copyproperties(source);
1454 if(renormalise)
1455 {
1456 volume4D<T> unitary_mask(source);
1457 unitary_mask=1;
1458 for (int t=source.mint(); t<=source.maxt(); t++) unitary_mask.replaceSubVolume(t,efficient_convolve(unitary_mask.constSubVolume(t),norm_kernel));
1459 result/=unitary_mask;
1460 }
1461 }
1462 return result;
1463 }
1464
1465
1467 // GRADIENT
1468
1469
1470 template <class T>
1471 volume<float> gradient(const volume<T>& source)
1472 {
1473 // in voxel coordinates (not mm)
1474 volume<float> maskx,masky,maskz;
1475 make_grad_masks(maskx,masky,maskz);
1476 volume<float> grad(source);
1477 float valx, valy, valz;
1478 int midx, midy, midz;
1479 midz=maskx.xsize()/2;
1480 midy=maskx.ysize()/2;
1481 midx=maskx.zsize()/2;
1482 for (int z=0; z<grad.zsize(); z++) {
1483 for (int y=0; y<grad.ysize(); y++) {
1484 for (int x=0; x<grad.xsize(); x++) {
1485 valx=0.0; valy=0.0; valz=0.0;
1486 for (int mz=-midz; mz<=midz; mz++) {
1487 for (int my=-midy; my<=midy; my++) {
1488 for (int mx=-midx; mx<=midx; mx++) {
1489 valx+=source(x+mx,y+my,z+mz) * maskx(mx+midx,my+midy,mz+midz);
1490 valy+=source(x+mx,y+my,z+mz) * masky(mx+midx,my+midy,mz+midz);
1491 valz+=source(x+mx,y+my,z+mz) * maskz(mx+midx,my+midy,mz+midz);
1492 }
1493 }
1494 }
1495 grad(x,y,z)=sqrt(MISCMATHS::Sqr(valx) + MISCMATHS::Sqr(valy) + MISCMATHS::Sqr(valz));
1496 }
1497 }
1498 }
1499 return grad;
1500 }
1501
1502
1503
1504 template <class T>
1505 void gradient(const volume<T>& source,volume<float>& grad)
1506 {
1507 volume<float> maskx,masky,maskz;
1508 make_grad_masks(maskx,masky,maskz);
1509 grad.reinitialize(source.xsize(),source.ysize(),source.zsize(),3);
1510 copybasicproperties(source,grad);
1511 float valx, valy, valz;
1512 int midx, midy, midz;
1513 midz=maskx.xsize()/2;
1514 midy=maskx.ysize()/2;
1515 midx=maskx.zsize()/2;
1516 for (int z=0; z<grad.zsize(); z++) {
1517 for (int y=0; y<grad.ysize(); y++) {
1518 for (int x=0; x<grad.xsize(); x++) {
1519 valx=0.0; valy=0.0; valz=0.0;
1520 for (int mz=-midz; mz<=midz; mz++) {
1521 for (int my=-midy; my<=midy; my++) {
1522 for (int mx=-midx; mx<=midx; mx++) {
1523 valx+=source(x+mx,y+my,z+mz) * maskx(mx+midx,my+midy,mz+midz);
1524 valy+=source(x+mx,y+my,z+mz) * masky(mx+midx,my+midy,mz+midz);
1525 valz+=source(x+mx,y+my,z+mz) * maskz(mx+midx,my+midy,mz+midz);
1526 }
1527 }
1528 }
1529 grad(x,y,z,0)=valx;
1530 grad(x,y,z,1)=valy;
1531 grad(x,y,z,2)=valz;
1532 }
1533 }
1534 }
1535
1536 }
1537
1538
1539
1540 template <class T>
1541 volume<float> lrxgrad(const volume<float>& im, const volume<T>& mask)
1542 {
1543 // calculates separate left and right gradients (with copying when at
1544 // borders of the mask) or zero for points outside of the mask
1545 // returns the two as part of a volume4D (vol[0] = left grad, vol[1] = right grad)
1546 volume4D<float> grad(im);
1547 grad.addvolume(im);
1548 if (!samesize(im,mask)) imthrow("Mask and image not the same size",20);
1549
1550 for (int z=0; z<im.zsize(); z++) {
1551 for (int y=0; y<im.ysize(); y++) {
1552 for (int x=0; x<im.xsize(); x++) {
1553 // default 0, only overridden if inside mask and can calculate the gradients
1554 if (mask(x,y,z)>0.5) {
1555 // left gradient
1556 if (x>0) {
1557 if (mask(x-1,y,z)>0.5) {
1558 grad(x,y,z,0) = im(x,y,z) - im(x-1,y,z);
1559 }
1560 }
1561 // right gradient
1562 if (x<im.xsize()-1) {
1563 if (mask(x+1,y,z)>0.5) {
1564 grad(x,y,z,1) = im(x+1,y,z) - im(x,y,z);
1565 } else {
1566 // if couldn't calculate right grad, then copy left
1567 grad(x,y,z,1) = grad(x,y,z,0);
1568 }
1569 }
1570 // if couldn't calculate left grad, then copy right
1571 if ( (x>0) && (mask(x-1,y,z)<=0.5) ) {
1572 grad(x,y,z,0) = grad(x,y,z,1);
1573 }
1574 // NB: if couldn't calculate either (but mask>0.5) then it is still 0
1575 }
1576 }
1577 }
1578 }
1579 return grad;
1580 }
1581
1582
1583 template <class T>
1584 volume<float> lrygrad(const volume<float>& im, const volume<T>& mask)
1585 {
1586 // calculates separate left and right gradients (with copying when at
1587 // borders of the mask) or zero for points outside of the mask
1588 // returns the two as part of a volume4D (vol[0] = left grad, vol[1] = right grad)
1589 volume4D<float> grad(im);
1590 grad.addvolume(im);
1591 if (!samesize(im,mask)) imthrow("Mask and image not the same size",20);
1592
1593 for (int z=0; z<im.zsize(); z++) {
1594 for (int y=0; y<im.ysize(); y++) {
1595 for (int x=0; x<im.xsize(); x++) {
1596 // default 0, only overridden if inside mask and can calculate the gradients
1597 if (mask(x,y,z)>0.5) {
1598 // left gradient
1599 if (y>0) {
1600 if (mask(x,y-1,z)>0.5) {
1601 grad(x,y,z,0) = im(x,y,z) - im(x,y-1,z);
1602 }
1603 }
1604 // right gradient
1605 if (y<im.ysize()-1) {
1606 if (mask(x,y+1,z)>0.5) {
1607 grad(x,y,z,1) = im(x,y+1,z) - im(x,y,z);
1608 } else {
1609 // if couldn't calculate right grad, then copy left
1610 grad(x,y,z,1) = grad(x,y,z,0);
1611 }
1612 }
1613 // if couldn't calculate left grad, then copy right
1614 if ( (y>0) && (mask(x,y-1,z)<=0.5) ) {
1615 grad(x,y,z,0) = grad(x,y,z,1);
1616 }
1617 // NB: if couldn't calculate either (but mask>0.5) then it is still 0
1618 }
1619 }
1620 }
1621 }
1622 return grad;
1623 }
1624
1625
1626 template <class T>
1627 volume<float> lrzgrad(const volume<float>& im, const volume<T>& mask)
1628 {
1629 // calculates separate left and right gradients (with copying when at
1630 // borders of the mask) or zero for points outside of the mask
1631 // returns the two as part of a volume4D (vol[0] = left grad, vol[1] = right grad)
1632 volume4D<float> grad(im);
1633 grad.addvolume(im);
1634 if (!samesize(im,mask)) imthrow("Mask and image not the same size",20);
1635 for (int z=0; z<im.zsize(); z++) {
1636 for (int y=0; y<im.ysize(); y++) {
1637 for (int x=0; x<im.xsize(); x++) {
1638 // default 0, only overridden if inside mask and can calculate the gradients
1639 if (mask(x,y,z)>0.5) {
1640 // left gradient
1641 if (z>0) {
1642 if (mask(x,y,z-1)>0.5) {
1643 grad(x,y,z,0) = im(x,y,z) - im(x,y,z-1);
1644 }
1645 }
1646 // right gradient
1647 if (z<im.zsize()-1) {
1648 if (mask(x,y,z+1)>0.5) {
1649 grad(x,y,z,1) = im(x,y,z+1) - im(x,y,z);
1650 } else {
1651 // if couldn't calculate right grad, then copy left
1652 grad(x,y,z,1) = grad(x,y,z,0);
1653 }
1654 }
1655 // if couldn't calculate left grad, then copy right
1656 if ( (z>0) && (mask(x,y,z-1)<=0.5) ) {
1657 grad(x,y,z,0) = grad(x,y,z,1);
1658 }
1659 // NB: if couldn't calculate either (but mask>0.5) then it is still 0
1660 }
1661 }
1662 }
1663 }
1664 return grad;
1665 }
1666
1667
1668
1669
1670 template <class T>
1671 volume<T> bandpass_temporal_filter(volume<T>& source,double hp_sigma, double lp_sigma)
1672 {
1673 int hp_mask_size_PLUS, lp_mask_size_PLUS, hp_mask_size_MINUS, lp_mask_size_MINUS;
1674 double *hp_exp=NULL, *lp_exp=NULL, *array, *array2;
1675 volume<T> result(source);
1676
1677 if (hp_sigma<=0) hp_mask_size_MINUS=0;
1678 else hp_mask_size_MINUS=(int)(hp_sigma*3); /* this isn't a linear filter, so small hard cutoffs at ends don't matter */
1679 hp_mask_size_PLUS=hp_mask_size_MINUS;
1680 if (lp_sigma<=0) lp_mask_size_MINUS=0;
1681 else lp_mask_size_MINUS=(int)(lp_sigma*20)+2; /* this will be small, so we might as well be careful */
1682 lp_mask_size_PLUS=lp_mask_size_MINUS;
1683
1684 array=new double[source.tsize()];
1685 array2=new double[source.tsize()];
1686
1687 if (hp_sigma>0)
1688 {
1689 hp_exp=new double[hp_mask_size_MINUS+hp_mask_size_PLUS+1];
1690 hp_exp+=hp_mask_size_MINUS;
1691 for(int t=-hp_mask_size_MINUS; t<=hp_mask_size_PLUS; t++)
1692 hp_exp[t] = exp( -0.5 * ((double)(t*t)) / (hp_sigma*hp_sigma) );
1693 }
1694
1695 if (lp_sigma>0)
1696 {
1697 double total=0;
1698 lp_exp=new double[lp_mask_size_MINUS+lp_mask_size_PLUS+1];
1699 lp_exp+=lp_mask_size_MINUS;
1700 for(int t=-lp_mask_size_MINUS; t<=lp_mask_size_PLUS; t++)
1701 {
1702 lp_exp[t] = exp( -0.5 * ((double)(t*t)) / (lp_sigma*lp_sigma) );
1703 total += lp_exp[t];
1704 }
1705
1706 for(int t=-lp_mask_size_MINUS; t<=lp_mask_size_PLUS; t++)
1707 lp_exp[t] /= total;
1708 }
1709 for(int z=0;z<source.zsize();z++)
1710 for(int y=0;y<source.ysize();y++)
1711 for(int x=0;x<source.xsize();x++)
1712 {
1713 for(int t=0; t<source.tsize(); t++)
1714 array[t] = (double)source.value(x,y,z,t);
1715 if (hp_sigma>0)
1716 {
1717 int done_c0=0;
1718 double c0=0;
1719 double mean(0);
1720 for(int t=0; t<source.tsize(); t++)
1721 {
1722 int tt;
1723 double c, w, A=0, B=0, C=0, D=0, N=0, tmpdenom;
1724 for(tt=MAX(t-hp_mask_size_MINUS,0); tt<=MIN(t+hp_mask_size_PLUS,source.tsize()-1); tt++)
1725 {
1726 int dt=tt-t;
1727 w = hp_exp[dt];
1728 A += w * dt;
1729 B += w * array[tt];
1730 C += w * dt * dt;
1731 D += w * dt * array[tt];
1732 N += w;
1733 }
1734 tmpdenom=C*N-A*A;
1735 if (tmpdenom!=0)
1736 {
1737 c = (B*C-A*D) / tmpdenom;
1738 if (!done_c0)
1739 {
1740 c0=c;
1741 done_c0=1;
1742 }
1743 array2[t] = c0 + array[t] - c;
1744 }
1745 else array2[t] = array[t];
1746 mean+=array2[t];
1747 }
1748 //Demean timeseries
1749 mean/=source.tsize();
1750 for(int t=0; t<source.tsize(); t++)
1751 array2[t]-=mean;
1752 memcpy(array,array2,sizeof(double)*source.tsize());
1753 }
1754 /* {{{ apply lowpass filter to 1D array */
1755 if (lp_sigma>0)
1756 {
1757 for(int t=0; t<source.tsize(); t++)
1758 {
1759 double total=0;
1760 double sum(0);
1761 int tt;
1762 for(tt=MAX(t-lp_mask_size_MINUS,0); tt<=MIN(t+lp_mask_size_PLUS,source.tsize()-1); tt++) {
1763 total += array[tt] * lp_exp[tt-t];
1764 sum+=lp_exp[tt-t];
1765 }
1766 if (sum>0)
1767 array2[t] = total/sum;
1768 else
1769 array2[t] = total;
1770 }
1771 memcpy(array,array2,sizeof(double)*source.tsize());
1772 }
1773 /* {{{ write 1D array back to input 4D data */
1774 for(int t=0; t<source.tsize(); t++)
1775 result.value(x,y,z,t)= (T)array[t] ;
1776
1777 }
1778 return result;
1779 }
1780
1781
1783 // EDGE DETECTION
1784
1785 /* detects closed contour edges in a volume as the zero crossings of the
1786 Laplacian of a Gaussian (LoG). This is implemented by convolving the
1787 data with a kernel formed by subtracting a Gaussian of radius sigma1
1788 from a second Gaussian kernel of radius sigma2 (sigma1 < sigma2) */
1789
1790 template <class T>
1791 volume<T> log_edge_detect(const volume<T>& source,
1792 float sigma1, float sigma2, int mode)
1793 {
1794
1795 /*
1796 * mode:
1797 * 0 : 3D edge detection
1798 * 1 : 2D edge detection in X plane
1799 * 2 : 2D edge detection in Y plane
1800 * 3 : 2D edge detection in Z plane
1801 */
1802
1803 int radius1;
1804 volume<float> log_kern, temp_kern;
1805 volume<T> log_result(source);
1806 volume<T> zero_crossing_result(source);
1807 zero_crossing_result = 0;
1808
1809 radius1 = (int)(4*sigma2);
1810 if (mode==0) {
1811 log_kern = gaussian_kernel2D(sigma2, radius1);
1812 temp_kern = gaussian_kernel2D(sigma1, radius1);
1813 } else {
1814 log_kern = gaussian_kernel3D(sigma2, radius1);
1815 temp_kern = gaussian_kernel3D(sigma1, radius1);
1816 }
1817
1818 log_kern -= temp_kern;
1819
1820 log_result = convolve(source, log_kern);
1821 for(int t=0;t<log_result.tsize();t++)
1822 for(int z=1;z<log_result.zsize()-1;z++)
1823 for(int y=1;y<log_result.ysize()-1;y++)
1824 for(int x=1;x<log_result.xsize()-1;x++) {
1825 float val = log_result.value(x,y,z,t);
1826 float ival = -val;
1827
1828 if (
1829 (val > 0.0) && (
1830 // 6 connectivity
1831 (log_result.value(x-1,y,z,t) <= ival) || // 1
1832 (log_result.value(x+1,y,z,t) <= ival) || // 2
1833 (log_result.value(x,y-1,z,t) <= ival) || // 3
1834 (log_result.value(x,y+1,z,t) <= ival) || // 4
1835 (log_result.value(x,y,z-1,t) <= ival) || // 5
1836 (log_result.value(x,y,z+1,t) <= ival) // 6
1837 )
1838 ){
1839
1840 zero_crossing_result(x,y,z,t) = 1.0;
1841 }
1842
1843 if (
1844 (val < 0.0) && (
1845 // 6 connectivity
1846 (log_result.value(x-1,y,z,t) > ival) || // 1
1847 (log_result.value(x+1,y,z,t) > ival) || // 2
1848 (log_result.value(x,y-1,z,t) > ival) || // 3
1849 (log_result.value(x,y+1,z,t) > ival) || // 4
1850 (log_result.value(x,y,z-1,t) > ival) || // 5
1851 (log_result.value(x,y,z+1,t) > ival) // 6
1852 )
1853 ){
1854
1855 zero_crossing_result(x,y,z,t) = 1.0;
1856 }
1857 }
1858
1859 return zero_crossing_result;
1860 }
1861
1862 template <class T>
1863 volume<T> fixed_edge_detect(const volume<T>& source, float threshold,
1864 bool twodimensional)
1865 {
1866 volume<T> result = source;
1867 int zsize = 3;
1868 if (twodimensional) zsize=1;
1869
1870 volume<float> log_kern(3,3,zsize);
1871 log_kern = -1;
1872 log_kern(1,1,(zsize-1)/2) = 8;
1873
1874 extrapolation oldex = source.getextrapolationmethod();
1875 source.setextrapolationmethod(mirror);
1876 result = convolve(source, log_kern);
1877 source.setextrapolationmethod(oldex);
1878 result.binarise(threshold);
1879
1880 return result;
1881 }
1882
1883
1884
1886 // EDGE STRENGTHEN (from avwmaths)
1887 template <class T>
1888 volume<T> edge_strengthen(const volume<T>& source)
1889 {
1890 float tmpf=2*sqrt(1/(pow((double)source.xdim(),2.0)) + 1/(pow((double)source.ydim(),2.0)) + 1/(pow((double)source.zdim(),2.0)));
1891 volume<T> result;
1892 result=source;
1893 if (source.zsize()>2)
1894 {
1895 for(int t=0;t<source.tsize();t++)
1896 for(int z=1;z<source.zsize()-1;z++)
1897 for(int y=1;y<source.ysize()-1;y++)
1898 for(int x=1;x<source.xsize()-1;x++)
1899 {
1900 double temp1 = pow(double(source.value(x,y,z+1,t)-source.value(x,y,z-1,t)),2.0)/ pow((double)source.zdim(),2.0);
1901 double temp2 = pow(double(source.value(x,y+1,z,t)-source.value(x,y-1,z,t)),2.0)/ pow((double)source.ydim(),2.0);
1902 double temp3 = pow(double(source.value(x+1,y,z,t)-source.value(x-1,y,z,t)),2.0)/ pow((double)source.xdim(),2.0);
1903 result(x,y,z,t)=(T)(sqrt(temp1+temp2+temp3)/tmpf);
1904 }
1905 }
1906 else
1907 {
1908 for(int t=0;t<source.tsize();t++)
1909 for(int z=0;z<source.zsize();z++)
1910 for(int y=1;y<source.ysize()-1;y++)
1911 for(int x=1;x<source.xsize()-1;x++)
1912 {
1913 double temp1=pow(double(source.value(x,y+1,z,t)-source.value(x,y-1,z,t)),2.0)/ pow((double)source.ydim(),2.0);
1914 double temp2=pow(double(source.value(x+1,y,z,t)-source.value(x-1,y,z,t)),2.0)/ pow((double)source.xdim(),2.0);
1915 result(x,y,z,t) = (T)(sqrt (temp1+temp2) / tmpf);
1916 }
1917 }
1918 return result;
1919 }
1920
1921
1923 // CONNECTED COMPONENTS
1924
1925 // support functions for connected components
1926
1927 int find_first_nonzero(const NEWMAT::Matrix& mat);
1928
1929 void addpair2set(int x, int y, std::vector<int>& sx, std::vector<int>& sy);
1930
1931 void relabel_components_uniquely(volume<int>& labelvol,
1932 const std::vector<int>& equivlista,
1933 const std::vector<int>& equivlistb, NEWMAT::ColumnVector& clustersizes);
1934
1935 void relabel_components_uniquely(volume<int>& labelvol,
1936 const std::vector<int>& equivlista,
1937 const std::vector<int>& equivlistb);
1938
1940 struct offset {
1941 int64_t x,y,z;
1942 offset(const int64_t ix,const int64_t iy,const int64_t iz) : x(ix),y(iy),z(iz) {}
1943 };
1944
1945 std::vector<offset> backConnectivity(int nDirections);
1946
1947 template <class T>
1948 void nonunique_component_labels(const volume<T>& vol,
1949 volume<int>& labelvol,
1950 std::vector<int>& equivlista,
1951 std::vector<int>& equivlistb,
1952 int numconnected)
1953 {
1954 copyconvert(vol,labelvol);
1955 labelvol = 0;
1956
1957 int labelnum(0);
1958 equivlista.erase(equivlista.begin(),equivlista.end());
1959 equivlistb.erase(equivlistb.begin(),equivlistb.end());
1960 std::vector<offset> neighbours(backConnectivity(numconnected));
1961 for (int z=vol.minz(); z<=vol.maxz(); z++)
1962 for (int y=vol.miny(); y<=vol.maxy(); y++)
1963 for (int x=vol.minx(); x<=vol.maxx(); x++) {
1964 T val(vol(x,y,z));
1965 if (val>0.5) { // The eligibility test
1966 int lval(labelvol(x,y,z));
1967 for (std::vector<offset>::iterator it=neighbours.begin();it!=neighbours.end();it++) {
1968 int xnew(x+it->x),ynew(y+it->y),znew(z+it->z);
1969 if ( (xnew>=vol.minx()) && (ynew>=vol.miny()) && (znew>=vol.minz())
1970 && (MISCMATHS::round(vol(xnew,ynew,znew)))==val) {
1971 // Binary relation
1972 int lval2 = labelvol(xnew,ynew,znew);
1973 if (lval != lval2) {
1974 if (lval!=0)
1975 addpair2set(lval2,lval,equivlista,equivlistb);
1976 labelvol(x,y,z) = lval2;
1977 lval = lval2;
1978 }
1979 }
1980 }
1981 if (lval==0)
1982 labelvol(x,y,z) = ++labelnum;
1983 }
1984 }
1985 }
1986
1987 template <class T>
1988 void nonunique_component_labels(const volume<T>& vol,
1989 volume<int>& labelvol,
1990 std::vector<int>& equivlista,
1991 std::vector<int>& equivlistb,
1992 bool (*binaryrelation)(T , T),
1993 int numconnected)
1994 {
1995 copyconvert(vol,labelvol);
1996 labelvol = 0;
1997
1998 int labelnum(0);
1999 equivlista.erase(equivlista.begin(),equivlista.end());
2000 equivlistb.erase(equivlistb.begin(),equivlistb.end());
2001 std::vector<offset> neighbours(backConnectivity(numconnected));
2002 for (int z=vol.minz(); z<=vol.maxz(); z++)
2003 for (int y=vol.miny(); y<=vol.maxy(); y++)
2004 for (int x=vol.minx(); x<=vol.maxx(); x++) {
2005 T val(vol(x,y,z));
2006 if (val>0.5) { // The eligibility test
2007 int lval(labelvol(x,y,z));
2008 for (std::vector<offset>::iterator it=neighbours.begin();it!=neighbours.end();it++) {
2009 int xnew(x+it->x),ynew(y+it->y),znew(z+it->z);
2010 if ((xnew>=vol.minx()) && (ynew>=vol.miny()) && (znew>=vol.minz())
2011 && ((*binaryrelation)(vol(xnew,ynew,znew),val))) {
2012 // Binary relation
2013 int lval2 = labelvol(xnew,ynew,znew);
2014 if (lval != lval2) {
2015 if (lval!=0)
2016 addpair2set(lval2,lval,equivlista,equivlistb);
2017 labelvol(x,y,z) = lval2;
2018 lval = lval2;
2019 }
2020 }
2021 }
2022 if (lval==0)
2023 labelvol(x,y,z) = ++labelnum;
2024 }
2025 }
2026 }
2027
2028 template <class T>
2029 volume<int> connected_components(const volume<T>& vol, NEWMAT::ColumnVector& clustersize, int numconnected)
2030 {
2031 volume<int> labelvol;
2032 copyconvert(vol,labelvol);
2033 std::vector<int> equivlista, equivlistb;
2034 nonunique_component_labels(vol,labelvol,
2035 equivlista,equivlistb,numconnected);
2036 relabel_components_uniquely(labelvol,equivlista,equivlistb,clustersize);
2037 return labelvol;
2038 }
2039
2040
2041 template <class T>
2042 volume<int> connected_components(const volume<T>& vol,
2043 const volume<T>& mask,
2044 bool (*binaryrelation)(T , T), NEWMAT::ColumnVector& clustersize)
2045 {
2046 volume<int> labelvol;
2047 copyconvert(vol,labelvol);
2048 std::vector<int> equivlista, equivlistb;
2049 nonunique_component_labels(vol,mask,labelvol,equivlista,equivlistb,
2050 binaryrelation);
2051 relabel_components_uniquely(labelvol,equivlista,equivlistb,clustersize);
2052 return labelvol;
2053 }
2054
2055
2056 template <class T>
2057 volume<int> connected_components(const volume<T>& vol, int numconnected){
2058 NEWMAT::ColumnVector clustersize;
2059 return connected_components(vol,clustersize,numconnected);
2060 }
2061
2062
2063 template <class T>
2064 volume<int> connected_components(const volume<T>& vol,const volume<T>& mask,
2065 bool (*binaryrelation)(T , T)){
2066 NEWMAT::ColumnVector clustersize;
2067 return connected_components(vol,mask,binaryrelation,clustersize);
2068 }
2069
2070
2072
2073
2074class rowentry { public: int x; int y; int z; float d; } ;
2075
2076bool rowentry_lessthan(const rowentry& r1, const rowentry& r2);
2077
2078template <class T>
2079class distancemapper {
2080private:
2081 const volume<T> &bvol;
2082 const volume<T> &mask;
2083 const volume<T> &bvolneg;
2084 std::vector<rowentry> schedule;
2085 std::vector<offset> octantsign;
2086 bool dualmasks;
2087public:
2088 // basic constructor takes binaryvol (mask of valid values)
2089 // + maskvol (non-zero at desired calculated points only)
2090 // Dual binaryvol version (maskpos and maskneg) for signed distance to nearest mask
2091 distancemapper(const volume<T>& binarypos, const volume<T>& binaryneg, const volume<T>& maskvol);
2092 distancemapper(const volume<T>& binaryvol, const volume<T>& maskvol);
2093 ~distancemapper();
2094 const volume<float> distancemap();
2095 volume4D<float> sparseinterpolate(const volume4D<float>& values,
2096 const std::string& interpmethod="general");
2097private:
2098 int setup_globals();
2099 int find_nearest(int x, int y, int z, int& x1, int& y1, int& z1,
2100 bool findav, NEWMAT::ColumnVector& localav, const volume4D<float>& vals);
2101 int find_nearest(int x, int y, int z, int& x1, int& y1, int& z1);
2102 int create_distancemap(volume4D<float>& vout, const volume4D<float>& valim,
2103 const std::string& interpmethod="none");
2104 int basic_create_distancemap(volume4D<float>& vout,
2105 const volume4D<float>& valim,
2106 const std::string& interpmethod);
2107};
2108
2109template <class T>
2110distancemapper<T>::distancemapper(const volume<T>& binarypos, const volume<T>& binaryneg, const volume<T>& maskvol) :
2111 bvol(binarypos), mask(maskvol), bvolneg(binaryneg)
2112{
2113 if (!samesize(bvol,mask)) imthrow("Mask and image not the same size",20);
2114 if (!samesize(bvol,bvolneg)) imthrow("Two binary images are not the same size",20);
2115 dualmasks=true;
2116 setup_globals();
2117}
2118
2119
2120template <class T>
2121distancemapper<T>::distancemapper(const volume<T>& binaryvol, const volume<T>& maskvol) :
2122 bvol(binaryvol), mask(maskvol), bvolneg(binaryvol)
2123{
2124 if (!samesize(bvol,mask)) imthrow("Mask and image not the same size",20);
2125 dualmasks=false;
2126 setup_globals();
2127}
2128
2129template <class T>
2130distancemapper<T>::~distancemapper()
2131{
2132 // destructors for the private members should do the job
2133}
2134
2135template <class T>
2136int distancemapper<T>::setup_globals()
2137{
2138 // octantsign gives the 8 different octant sign combinations for coord
2139 // offsets
2140 for (int p=-1; p<=1; p+=2)
2141 for (int q=-1; q<=1; q+=2)
2142 for (int r=-1; r<=1; r+=2)
2143 octantsign.push_back(offset(p,q,r));
2144
2145 // construct list of displacements (in one octant) in ascending
2146 // order of distance
2147 for (int z=bvol.minz(); z<=bvol.maxz(); z++) {
2148 for (int y=bvol.miny(); y<=bvol.maxy(); y++) {
2149 for (int x=bvol.minx(); x<=bvol.maxx(); x++) {
2150 rowentry newrow;
2151 newrow.x=x;
2152 newrow.y=y;
2153 newrow.z=z;
2154 float d2 = MISCMATHS::norm2sq(x*bvol.xdim(),y*bvol.ydim(),z*bvol.zdim());
2155 newrow.d=d2;
2156 schedule.push_back(newrow);
2157 }
2158 }
2159 }
2160
2161 // sort schedule to get ascending d2
2162 sort(schedule.begin(),schedule.end(),NEWIMAGE::rowentry_lessthan);
2163 return 0;
2164}
2165
2166// findav determines whether to do interpolation calculations or just
2167// return the location only
2168template <class T>
2169int distancemapper<T>::find_nearest(int x, int y, int z, int& x1, int& y1, int& z1,
2170 bool findav, NEWMAT::ColumnVector& localav,
2171 const volume4D<float>& vals)
2172{
2173 float sumw=0.0, mindist=0, maxdist=0.0, weight;
2174 float minVoxSize = std::min(bvol.xdim(),std::min(bvol.ydim(),bvol.zdim()));
2175 NEWMAT::ColumnVector sumvw;
2176 if (findav) {
2177 localav.ReSize(vals.tsize());
2178 localav=0.0;
2179 sumvw.ReSize(vals.tsize());
2180 sumvw=0.0;
2181 }
2182 for (std::vector<rowentry>::iterator it=schedule.begin(); it!=schedule.end();it++) {
2183 for (std::vector<offset>::iterator oit=octantsign.begin(); oit!=octantsign.end();oit++) {
2184 x1=x+it->x*oit->x;
2185 y1=y+it->y*oit->y;
2186 z1=z+it->z*oit->z;
2187 if (bvol.in_bounds(x1,y1,z1)) {
2188 bool primaryValid( bvol.value(x1,y1,z1) > 0.5);
2189 if ( primaryValid || ( dualmasks && (bvolneg.value(x1,y1,z1)>0.5))) {
2190 if ((!findav) || (dualmasks)) {
2191 return primaryValid ? 1 : -1; //Either primary or secondary is valid, so if primary is false...
2192 } else if (mindist==0.0) { // first time a point is encountered
2193 mindist=it->d;
2194 // select distance band to average over (farther -> more)
2195 maxdist=MISCMATHS::Max(mindist+minVoxSize,mindist*1.5);
2196 } else if (it->d>maxdist) { // stop after maxdist reached
2197 localav=sumvw/sumw;
2198 return 1;
2199 }
2200 weight = it->d ? mindist/it->d : 1;
2201 sumw += weight;
2202 for (int t=0; t<vals.tsize(); t++)
2203 sumvw(t+1) += weight * vals.value(x1,y1,z1,t);
2204 }
2205 }
2206 }
2207 }
2208 // return most distant point (as last resort)
2209 if (findav) { if (sumw>0) { localav=sumvw/sumw; return 1; }}//should be (findav) and not not?
2210 return 0; // not found any binary voxel (true for a zero image)
2211}
2212
2213
2214template <class T>
2215int distancemapper<T>::find_nearest(int x, int y, int z, int& x1, int& y1, int& z1)
2216{
2217 NEWMAT::ColumnVector dummy;
2218 volume4D<float> dummyvol;
2219 return this->find_nearest(x,y,z,x1,y1,z1,false,dummy,dummyvol);
2220}
2221
2222
2223// create the distance map as vout
2224// if return_distance is false then interpolate the value of the input volume
2225// at the output location rather than store the distance to it
2226template <class T>
2227int distancemapper<T>::create_distancemap(volume4D<float>& vout,
2228 const volume4D<float>& valim,
2229 const std::string& interpmethod)
2230{
2231 if (interpmethod!="general")
2232 return this->basic_create_distancemap(vout,valim,interpmethod);
2233 // only get to here for sparseinterpolation
2234 float meanvoxsize = pow(valim.xdim()*valim.ydim()*valim.zdim(),1.0/3.0);
2235 int nsubsamp=0;
2236 if (meanvoxsize<4.0) { nsubsamp=1; }
2237 if (meanvoxsize<2.0) { nsubsamp=2; }
2238 if (meanvoxsize<1.0) { nsubsamp=3; }
2239 // for the straightforward case (>4mm resolution)
2240 if (nsubsamp==0) {
2241 return basic_create_distancemap(vout,valim,interpmethod);
2242 }
2243 // otherwise run subsamplings, sparseinterp and upsamplings
2244 // this is to fill in the large gaps in the image, since this
2245 // takes a *huge* amount of time otherwise
2246 // NB: mask in distancemapper is 1 where the new values are to go
2247 // but in this function we use mask=1 where valid samples are
2248 std::vector<volume4D<float> > sampledData;
2249 std::vector<volume<float> > sampledMask;
2250 sampledData.resize(nsubsamp+1);
2251 sampledMask.resize(nsubsamp+1);
2252 sampledMask[0]=1.0f - mask; // sampledMask is now 1 for all valid sample points
2253 sampledData[0]=valim*sampledMask[0];
2254 for (int n=1; n<=nsubsamp; n++) {
2255 sampledMask[n]=subsample_by_2(sampledMask[n-1],false);
2256 sampledData[n].reinitialize(sampledMask[n].xsize(),sampledMask[n].ysize(),sampledMask[n].zsize(),valim.tsize());
2257 sampledData[n].copyproperties(sampledMask[n]);
2258 for (int t=0; t<valim.tsize(); t++)
2259 (sampledData[n])[t]=divide(subsample_by_2((sampledData[n-1])[t],false),sampledMask[n],sampledMask[n]);
2260
2261 sampledMask[n].binarise(1e-4); // include any voxel with any partial overlap
2262 }
2263 // run the sparseinterpolate function (at 8mm-ish resolution)
2264 // - not the method in this object, but a new one
2265 sampledData[nsubsamp]=NEWIMAGE::sparseinterpolate(sampledData[nsubsamp],sampledMask[nsubsamp]);
2266 // dilate and invert mask
2267 volume<float> kernel(3,3,3);
2268 kernel=1.0f;
2269 sampledMask[nsubsamp]=morphfilter(sampledMask[nsubsamp],kernel,"dilate");
2270 sampledMask[nsubsamp]= 1.0f - sampledMask[nsubsamp]; // now the mask is 1 for all the more distant "gaps"
2271 // upsample mask and spareinterpolated result
2272 volume<float> tmp;
2273 for (int n=nsubsamp; n>0; n--) {
2274 for (int t=0; t<valim.tsize(); t++) {
2275 ShadowVolume<float> shadow(sampledData[n-1][t]);
2276 upsample_by_2(shadow,sampledData[n][t],false);
2277 }
2278 ShadowVolume<float> shadow(sampledMask[n-1]);
2279 upsample_by_2(shadow,sampledMask[n],false);
2280 }
2281 sampledMask[0].binarise(0.5f);
2282 // only keep results in zero part of mask as well as the
2283 // original points
2284 volume4D<float> vres(valim*(1.0f-mask) + sampledData[0]*sampledMask[0]);
2285 sampledMask[0] = mask - sampledMask[0]; // only calculate new points in the "gap"
2286 // now run basic_create_distancemap at full resolution, but with the
2287 // new mask extra filled-in areas from above
2288 volume<float> invMask(1.0f - sampledMask[0]);
2289 distancemapper<float> newdmapper(invMask,sampledMask[0]);
2290 return newdmapper.basic_create_distancemap(vout,vres,interpmethod);
2291}
2292
2293// The following function creates the distancemap or interpolated image
2294// from valim (previously the mask of valid voxels = binaryvol and the
2295// mask of voxels to calculate = maskvol, must have been specified)
2296template <class T>
2297int distancemapper<T>::basic_create_distancemap(volume4D<float>& vout,
2298 const volume4D<float>& valim,
2299 const std::string& interpmethod)
2300{
2301 int x1, y1, z1;
2302 NEWMAT::ColumnVector localav;
2303 int interp=0, distsign=1;
2304 if ((interpmethod=="nn") || (interpmethod=="nearestneighbour")) interp=1;
2305 if (interpmethod=="general") interp=2;
2306 if (interp>0) { vout = valim; } else { vout = bvol; vout *= 0.0f; }
2307 if ((interp>0) && (!samesize(bvol,valim.constSubVolume(0))))
2308 { print_volume_info(bvol,std::string("bvol"),std::cerr); print_volume_info(mask,std::string("mask"),std::cerr); print_volume_info(valim.constSubVolume(0),std::string("valim"),std::cerr);
2309 imthrow("Binary image and interpolant not the same size",21); }
2310 for (int z=0; z<vout.zsize(); z++) {
2311 for (int y=0; y<vout.ysize(); y++) {
2312 for (int x=0; x<vout.xsize(); x++) {
2313 if (mask(x,y,z)>((T) 0.5)) {
2314 distsign=find_nearest(x,y,z,x1,y1,z1,interp>=2,localav,valim);
2315 switch (interp) {
2316 case 2:
2317 for (int t=0;t<valim.tsize();t++) { vout(x,y,z,t)=localav(t+1); }
2318 break;
2319 case 1:
2320 for (int t=0;t<valim.tsize();t++) { vout(x,y,z,t)=valim(x1,y1,z1,t); }
2321 break;
2322 case 0:
2323 default:
2324 vout(x,y,z,0)=distsign*sqrt(MISCMATHS::norm2sq((x1-x)*bvol.xdim(),
2325 (y1-y)*bvol.ydim(),(z1-z)*bvol.zdim()));
2326 }
2327 }
2328 }
2329 }
2330 }
2331 return 0;
2332}
2333
2334
2335template <class T>
2336const volume<float> distancemapper<T>::distancemap()
2337{
2338 volume4D<float> dmap;
2339 create_distancemap(dmap,dmap,"none");
2340 return dmap[0];
2341}
2342
2343template <class T>
2344volume4D<float> distancemapper<T>::sparseinterpolate(const volume4D<float>& values,
2345 const std::string& interpmethod)
2346{
2347 volume4D<float> vout;
2348 create_distancemap(vout,values,interpmethod);
2349 return vout;
2350}
2351
2352
2353template <class T>
2354volume<float> distancemap(const volume<T>& binaryvol)
2355{
2356 volume<T> mask;
2357 mask = ((T) 1) - binarise(binaryvol,((T) 0.5));
2358 return distancemap(binaryvol,mask);
2359}
2360
2361template <class T>
2362volume<float> distancemap(const volume<T>& binaryvol, const volume<T>& mask)
2363{
2364 distancemapper<T> dmapper(binaryvol,mask);
2365 return dmapper.distancemap();
2366}
2367
2368template <class T>
2369volume<float> distancemap(const volume<T>& binarypos, const volume<T>& binaryneg, const volume<T>& maskvol)
2370{
2371 distancemapper<T> dmapper(binarypos,binaryneg,maskvol);
2372 return dmapper.distancemap();
2373}
2374
2375template <class T>
2376volume4D<float> sparseinterpolate(const volume4D<T>& sparsesamps,
2377 const volume<T>& mask,
2378 const std::string& interpmethod)
2379{
2380 // can have "general" or "nearestneighbour" (or "nn") for interpmethod
2381 volume<T> invmask;
2382 invmask=((T) 1) - mask;
2383 distancemapper<T> dmapper(mask,invmask);
2384 return dmapper.sparseinterpolate(sparsesamps,interpmethod);
2385}
2386
2387
2389
2390template <class T>
2391void tfce_orig_slow(volume<T>& VolIntn, float H, float E, int NumConn, float minT, float deltaT)
2392{
2393 float maxval=VolIntn.max();
2394 volume<float> clusterenhance;
2395 copyconvert(VolIntn,clusterenhance);
2396 clusterenhance=0;
2397
2398 if (deltaT==0)
2399 deltaT = (maxval - minT)/100.0; // this needs fixing!!
2400
2401 for (float thresh=minT+deltaT; thresh<=maxval; thresh+=deltaT)
2402 {
2403 volume<float> clusters;
2404 copyconvert(VolIntn,clusters);
2405 clusters.binarise(thresh);
2406
2407 NEWMAT::ColumnVector clustersizes;
2408 volume<int>tmpvol=connected_components(clusters,clustersizes,NumConn);
2409 clustersizes = pow(clustersizes,E) * pow(thresh,H);
2410 for(int z=0;z<VolIntn.zsize();z++)
2411 for(int y=0;y<VolIntn.ysize();y++)
2412 for(int x=0;x<VolIntn.xsize();x++)
2413 if (tmpvol.value(x,y,z)>0)
2414 clusterenhance.value(x,y,z) += clustersizes(tmpvol.value(x,y,z));
2415 }
2416 copyconvert(clusterenhance,VolIntn);
2417 return;
2418}
2419
2420class VecSort{
2421 public:
2422 int Sx, Sy, Sz, Sl;
2423 double Sv;
2424 bool operator<(const VecSort& other) const{
2425 return Sv < other.Sv;
2426 }
2427};
2428
2429template <class T>
2430void tfce(volume<T>& data, float H, float E, int NumConn, float minT, float deltaT)
2431{
2432 volume<int> VolLabl; copyconvert(data, VolLabl);
2433 volume<float> VolEnhn; copyconvert(data, VolEnhn); VolEnhn=0;
2434 bool doIT=false;
2435 const int INIT=-1, MASK=-2;
2436 int curlab=0;
2437 int pX, pY, pZ, qX, qY, qZ, rX, rY, rZ;
2438 int FldCntr=0, FldCntri=0, tfceCntr=0, xFC = 0;
2439 int minX=1, minY=1, minZ=1, maxX=data.maxx()-1, maxY=data.maxy()-1, maxZ=data.maxz()-1;
2440 int sizeC=maxX*maxY*maxZ;
2441 int counter=0, edsta[27];
2442 float maxT=data.max();
2443 if(deltaT==0)
2444 deltaT=maxT/100.0;
2445 if(deltaT<=0)
2446 throw std::runtime_error("Error: tfce requires a positive deltaT input.");
2447 if ( data.xsize() < 3 || data.ysize() < 3 || data.zsize() < 3 )
2448 throw std::runtime_error("Error: tfce currently requires an input with at least 3 voxels extent into each dimension.");
2449 if( data.max()/deltaT > 10000 )
2450 std::cout << "Warning: tfce has detected a large number of integral steps. This operation may require a great deal of time to complete." << std::endl;
2451 std::vector<VecSort> VecSortI(sizeC);
2452 std::queue<int> Qx, Qy, Qz;
2453 for(int z0=-1; z0<=1; z0++)
2454 for(int y0=-1; y0<=1; y0++)
2455 for(int x0=-1; x0<=1; x0++){
2456 edsta[counter++] = (x0*x0+y0*y0+z0*z0);
2457 if (edsta[counter-1]<2 && edsta[counter-1]!=0) edsta[counter-1]=6;
2458 if (edsta[counter-1]<3 && edsta[counter-1]!=0) edsta[counter-1]=18;
2459 if (edsta[counter-1]<4 && edsta[counter-1]!=0) edsta[counter-1]=26;
2460 if (edsta[counter-1]==0) edsta[counter-1]=100;
2461 }
2462 FldCntr=0;
2463 for(int z=minZ; z<=maxZ; z++)
2464 for(int y=minY; y<=maxY; y++)
2465 for(int x=minX; x<=maxX; x++) {
2466 float iVal=data.value(x,y,z);
2467 if( iVal > minT ) {
2468 VecSortI[FldCntr].Sx=x; VecSortI[FldCntr].Sy=y; VecSortI[FldCntr].Sz=z;
2469 VecSortI[FldCntr++].Sv=iVal;
2470 }
2471 }
2472 sizeC=FldCntr;
2473 VecSortI.resize(sizeC);
2474 sort(VecSortI.begin(), VecSortI.end());
2475 for(float curThr=minT; curThr<(maxT+deltaT);curThr+=deltaT){
2476 VolLabl = INIT; FldCntr = xFC; curlab = 0;
2477 while( (VecSortI[FldCntr].Sv<=curThr) && (FldCntr<sizeC) ){
2478 VecSortI[FldCntr].Sl=0; VecSortI[FldCntr++].Sv=0;
2479 }
2480 xFC=FldCntr;
2481 while( VecSortI[FldCntr].Sv>curThr && (FldCntr<sizeC) ){
2482 pX=VecSortI[FldCntr].Sx; pY=VecSortI[FldCntr].Sy; pZ=VecSortI[FldCntr++].Sz;
2483 VolLabl.value(pX,pY,pZ)=MASK;
2484 }
2485 for(FldCntri=xFC; FldCntri<FldCntr; FldCntri++){
2486 pX=VecSortI[FldCntri].Sx; pY=VecSortI[FldCntri].Sy; pZ=VecSortI[FldCntri].Sz;
2487 if(VolLabl.value(pX, pY, pZ)==MASK){//sI
2488 curlab+=1;
2489 Qx.push(pX); Qy.push(pY); Qz.push(pZ);
2490 VolLabl.value(pX, pY, pZ)=curlab;
2491 while(!Qx.empty()){//sW
2492 qX=Qx.front(); qY=Qy.front(); qZ=Qz.front();
2493 Qx.pop(); Qy.pop(); Qz.pop();
2494 counter=0;
2495 for(int z0=-1; z0<=1; z0++)
2496 for(int y0=-1; y0<=1; y0++)
2497 for(int x0=-1; x0<=1; x0++){
2498 rX=qX+x0; rY=qY+y0; rZ=qZ+z0;
2499 doIT =(NumConn>=edsta[counter++]);
2500 if(doIT && (VolLabl.value(rX, rY, rZ)==MASK)){
2501 Qx.push(rX); Qy.push(rY);Qz.push(rZ);
2502 VolLabl.value(rX, rY, rZ)=curlab;
2503 }
2504 }
2505 }//eW
2506 }//eI
2507 }
2508 NEWMAT::ColumnVector ClusterSizes(curlab), ClusterSizesI(curlab);
2509 ClusterSizes=0;
2510 for(tfceCntr=xFC; tfceCntr<sizeC; tfceCntr++){
2511 VecSortI[tfceCntr].Sl=VolLabl.value(VecSortI[tfceCntr].Sx, VecSortI[tfceCntr].Sy, VecSortI[tfceCntr].Sz);
2512 if ( VecSortI[tfceCntr].Sl>0 )
2513 ClusterSizes(VecSortI[tfceCntr].Sl)+=1;
2514 }
2515 float HH=pow(curThr, H);
2516 ClusterSizesI=pow(ClusterSizes, E)*HH;
2517 for(tfceCntr=xFC; tfceCntr<sizeC; tfceCntr++){
2518 if ( VecSortI[tfceCntr].Sl>0 ){
2519 VolEnhn.value(VecSortI[tfceCntr].Sx, VecSortI[tfceCntr].Sy, VecSortI[tfceCntr].Sz) += ClusterSizesI(VecSortI[tfceCntr].Sl);
2520 }
2521 }
2522 }//end curThr
2523 copyconvert(VolEnhn,data);
2524 return;
2525}
2526
2527template <class T>
2528void tfce_support(volume<T>& VolIntn, float H, float E, int NumConn, float minT, float deltaT, int Xoi, int Yoi, int Zoi, float threshTFCE)
2529{
2530 volume<float> VolEnhn; copyconvert(VolIntn, VolEnhn);
2531 volume<float> VolTemp; copyconvert(VolIntn, VolTemp);
2532 int minX=0, minY=0, minZ=0, maxX=VolIntn.maxx(), maxY=VolIntn.maxy(), maxZ=VolIntn.maxz();
2533 float maxT = VolIntn.value(Xoi,Yoi,Zoi);
2534 int thrCntr = 0; int thrNum = 0;
2535 if(deltaT==0){
2536 deltaT = maxT/100;
2537 thrNum = 101;
2538 }
2539 else
2540 thrNum = int(maxT/deltaT)+1;
2541 NEWMAT::ColumnVector Clusters(thrNum), Thresholds(thrNum), ClusterSizes;
2542 for(float thresh=minT; thresh<maxT; thresh+=deltaT){
2543 copyconvert(VolEnhn, VolTemp);
2544 VolTemp.binarise(thresh);
2545 volume<int> VolLabl = connected_components(VolTemp, ClusterSizes, NumConn);
2546 Clusters(thrCntr+1) = ClusterSizes(VolLabl.value(Xoi, Yoi, Zoi));
2547 Thresholds(thrCntr+1) = thresh;
2548 for(int z=minZ; z<=maxZ; z++)
2549 for(int y=minY; y<=maxY; y++)
2550 for(int x=minX; x<=maxX; x++)
2551 if( VolLabl.value(x,y,z) != VolLabl.value(Xoi,Yoi,Zoi) )
2552 VolEnhn.value(x,y,z) = 0;
2553 thrCntr++;
2554 }
2555 float summ=0, thre=0;
2556 for (int i=0; i<thrCntr; i++){
2557 summ+=std::pow((float)Clusters(thrCntr-i), E)*std::pow((float)Thresholds(thrCntr-i), H);
2558 thre = Thresholds(thrCntr-i);
2559 if(summ>=threshTFCE)
2560 break;
2561 }
2562 if(summ<threshTFCE)
2563 std::cout<<"it doesn't reach to specified threshold"<<std::endl;
2564 copyconvert(VolIntn, VolEnhn);
2565 copyconvert(VolIntn, VolTemp); VolTemp.binarise(thre);
2566 volume<int> VolLabl = connected_components(VolTemp, ClusterSizes, NumConn);
2567 for(int z=minZ; z<=maxZ; z++)
2568 for(int y=minY; y<=maxY; y++)
2569 for(int x=minX; x<=maxX; x++)
2570 if( VolLabl.value(x,y,z)!=VolLabl(Xoi, Yoi, Zoi) )
2571 VolEnhn.value(x,y,z) = 0;
2572 copyconvert(VolEnhn,VolIntn);
2573 return;
2574}
2575
2576
2577
2578
2580
2581}
2582
2583
2584#endif
NEWIMAGE::VecSort
Definition: newimagefns.h:2420
NEWIMAGE::distancemapper
Definition: newimagefns.h:2079
NEWIMAGE::rowentry
Definition: newimagefns.h:2074
NEWIMAGE::volume
Definition: newimage.h:100
NEWIMAGE::offset
Definition: newimagefns.h:1940
NEWIMAGE::vec3_temp
Definition: newimagefns.h:1198