1 #include "mrilib.h"
2
3 #ifdef USE_OMP /* don't use this yet */
4 #include <omp.h>
5 #endif
6
7 #undef INMASK
8 #define INMASK(i) ( mask == NULL || mask[i] != 0 )
9
10 static int myTHD_extract_nbhd( THD_3dim_dataset *dset, byte *mask,
11 int xx, int yy, int zz, MCW_cluster *nbhd,
12 int *ivar , float *tsar ) ;
13
14 #ifdef USE_OMP
15 #include "cs_pv.c" /* to ensure it's OpenMP-ized as well */
16 #endif
17
18 /*------------------------------------------------------------------------*/
19
vstep_print(void)20 static void vstep_print(void)
21 {
22 static char xx[10] = "0123456789" ; static int vn=0 ;
23 fprintf(stderr , "%c" , xx[vn%10] ) ;
24 if( vn%10 == 9) fprintf(stderr,".") ;
25 vn++ ;
26 }
27
28 /*------------------------------------------------------------------------*/
29 /* Adapted from 3dLocalSVD, to be more efficient when #evals=1 */
30
main(int argc,char * argv[])31 int main( int argc , char *argv[] )
32 {
33 THD_3dim_dataset *inset=NULL , *outset=NULL , *evset=NULL , *outset2=NULL ;
34 MCW_cluster *nbhd=NULL ;
35 byte *mask=NULL ; int mask_nx=0,mask_ny=0,mask_nz=0 , automask=0 ;
36 char *prefix ="./LocalPV" ;
37 char *prefix2=NULL ;
38 char *evprefix=NULL ; int nev ;
39 int iarg=1 , verb=1 , ntype=0 , kk,nx,ny,nz,nxy,nxyz,nt , vstep=0 ;
40 float na,nb,nc , dx,dy,dz ;
41 int do_vnorm=0 , do_vproj=0 , polort=-1 ;
42 int nmask=0 , domean=0 , use_nonmask=0 ;
43 float *evar=NULL,*fvar=NULL ;
44 unsigned int gseed ;
45 int despike=0 ; /* 14 Oct 2010 */
46
47 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
48
49 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
50 printf(
51 "Usage: 3dLocalPV [options] inputdataset\n"
52 "* You may want to use 3dDetrend before running this program,\n"
53 " or at least use the '-polort' option.\n"
54 "* This program is highly experimental. And slowish. Real slowish.\n"
55 "* Computes the SVD of the time series from a neighborhood of each\n"
56 " voxel. An inricate way of 'smoothing' 3D+time datasets, kind of, sort of.\n"
57 "* This is like 3dLocalSVD, except that the '-vproj' option doesn't\n"
58 " allow anything but 1 and 2 dimensional projection. This is because\n"
59 " 3dLocalPV uses a special method to compute JUST the first 1 or 2\n"
60 " principal vectors -- faster than 3dLocalSVD, but less general.\n"
61 "\n"
62 "Options:\n"
63 " -mask mset = restrict operations to this mask\n"
64 " -automask = create a mask from time series dataset\n"
65 " -prefix ppp = save SVD vector result into this new dataset\n"
66 " [default = 'LocalPV']\n"
67 " -prefix2 qqq = save second principal vector into this new dataset\n"
68 " [default = don't save it]\n"
69 " -evprefix ppp = save singular value at each voxel into this dataset\n"
70 " [default = don't save]\n"
71 " -input inputdataset = input time series dataset\n"
72 " -nbhd nnn = e.g., 'SPHERE(5)' 'TOHD(7)' etc.\n"
73 " -despike = remove time series spikes from input dataset\n"
74 " -polort p = detrending\n"
75 " -vnorm = normalize data vectors [strongly recommended]\n"
76 " -vproj [2] = project central data time series onto local SVD vector;\n"
77 " if followed by '2', then the central data time series\n"
78 " will be projected on the 2-dimensional subspace\n"
79 " spanned by the first 2 principal SVD vectors.\n"
80 " [default: just output principal singular vector]\n"
81 " [for 'smoothing' purposes, '-vnorm -vproj' is an idea]\n"
82 #if 0 /* this option not actually implemented */
83 "\n"
84 " -use_nonmask = Allow the computation of the local SVD time series\n"
85 " even from voxels that are NOT in the mask, provided\n"
86 " that there are voxels IN the mask inside the local\n"
87 " neighborhood.\n"
88 " * You could use '-use_nonmask' to compute the principal\n"
89 " SVD vector of local white matter time series, for\n"
90 " example, even at non-WM voxels.\n"
91 " * '-vproj' is not allowed with '-use_nonmask'!\n"
92 #endif
93 "\n"
94 "Notes:\n"
95 "* On my Mac Pro, about 30%% faster than 3dLocalSVD computing the same thing.\n"
96 "* If you're curious, the 'special method' used for the eigensolution is\n"
97 " a variant of matrix power iteration, called 'simultaneous iteration'.\n"
98 "* This method uses pseudo-random numbers to initialize the vector iterations.\n"
99 " If you wish to control that seed, set environment variable\n"
100 " AFNI_RANDOM_SEEDVAL to some nonzero number. Otherwise, a random seed will\n"
101 " be selected from the time, which means otherwise identical runs will give\n"
102 " slightly different results.\n"
103 "* By contrast, 3dLocalSVD uses EISPACK functions for eigensolution-izing.\n"
104 ) ;
105 PRINT_AFNI_OMP_USAGE("3dLocalPV",NULL) ;
106 PRINT_COMPILE_DATE ; exit(0) ;
107 }
108
109 /*---- official startup ---*/
110
111 PRINT_VERSION("3dLocalPV"); mainENTRY("3dLocalPV main"); machdep();
112 AFNI_logger("3dLocalPV",argc,argv); AUTHOR("Emperor Zhark the Iterator");
113
114 gseed = (unsigned int)AFNI_numenv("AFNI_RANDOM_SEEDVAL") ;
115 if( gseed == 0 )
116 gseed = ((unsigned int)time(NULL)) + 17*(unsigned int)getpid() ;
117 INFO_message("random seed set to %u from AFNI_RANDOM_SEEDVAL",gseed) ;
118
119 /*---- loop over options ----*/
120
121 while( iarg < argc && argv[iarg][0] == '-' ){
122
123 if( strcmp(argv[iarg],"-use_nonmask") == 0 ){ /* 17 Jul 2009 */
124 use_nonmask = 1 ; iarg++ ; continue ;
125 }
126 if( strcmp(argv[iarg],"-vnorm") == 0 ){
127 do_vnorm = 1 ; iarg++ ; continue ;
128 }
129 if( strcmp(argv[iarg],"-vproj") == 0 ){
130 if( iarg+1 < argc && isdigit(argv[iarg+1][0]) )
131 do_vproj = (int)strtod(argv[++iarg],NULL) ;
132 if( do_vproj < 1 ){
133 do_vproj = 1 ;
134 } else if( do_vproj > 2 ){
135 do_vproj = 2 ; WARNING_message("-vproj set to 2") ;
136 }
137 iarg++ ; continue ;
138 }
139
140 if( strcmp(argv[iarg],"-despike") == 0 ){ /* 14 Oct 2010 */
141 despike = 1 ; iarg++ ; continue ;
142 }
143
144 if( strcmp(argv[iarg],"-polort") == 0 ){
145 char *cpt ;
146 if( ++iarg >= argc ) ERROR_exit("Need argument after '-polort'") ;
147 polort = (int)strtod(argv[iarg],&cpt) ;
148 if( *cpt != '\0' ) WARNING_message("Illegal non-numeric value after -polort") ;
149 iarg++ ; continue ;
150 }
151
152 if( strcmp(argv[iarg],"-input") == 0 ){
153 if( inset != NULL ) ERROR_exit("Can't have two -input options") ;
154 if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;
155 inset = THD_open_dataset( argv[iarg] ) ;
156 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
157 iarg++ ; continue ;
158 }
159
160 if( strcmp(argv[iarg],"-prefix") == 0 ){
161 if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix'") ;
162 prefix = strdup(argv[iarg]) ;
163 if( ! THD_filename_ok(prefix) )
164 ERROR_exit("-prefix '%s' is illegal",prefix) ;
165 iarg++ ; continue ;
166 }
167
168 if( strcmp(argv[iarg],"-prefix2") == 0 ){
169 if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix2'") ;
170 prefix2 = strdup(argv[iarg]) ;
171 if( ! THD_filename_ok(prefix2) )
172 ERROR_exit("-prefix '%s' is illegal",prefix2) ;
173 iarg++ ; continue ;
174 }
175
176 if( strcmp(argv[iarg],"-evprefix") == 0 ){
177 if( ++iarg >= argc ) ERROR_exit("Need argument after '-evprefix'") ;
178 evprefix = strdup(argv[iarg]) ;
179 iarg++ ; continue ;
180 }
181
182 if( strcmp(argv[iarg],"-mask") == 0 ){
183 THD_3dim_dataset *mset ;
184 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
185 if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;
186 mset = THD_open_dataset( argv[iarg] ) ;
187 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
188 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
189 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
190 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
191 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;
192 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
193 INFO_message("Number of voxels in mask = %d",nmask) ;
194 if( nmask < 2 ) ERROR_exit("Mask is too small to process") ;
195 iarg++ ; continue ;
196 }
197
198 if( strcmp(argv[iarg],"-automask") == 0 ){
199 if( mask != NULL ) ERROR_exit("Can't have -automask and -mask") ;
200 automask = 1 ;
201 iarg++ ; continue ;
202 }
203
204 if( strcmp(argv[iarg],"-nbhd") == 0 ){
205 char *cpt ;
206 if( ntype > 0 ) ERROR_exit("Can't have 2 '-nbhd' options") ;
207 if( ++iarg >= argc ) ERROR_exit("Need argument after '-nbhd'") ;
208
209 cpt = argv[iarg] ;
210 if( strncasecmp(cpt,"SPHERE",6) == 0 ){
211 sscanf( cpt+7 , "%f" , &na ) ;
212 if( na == 0.0f ) ERROR_exit("Can't have a SPHERE of radius 0") ;
213 ntype = NTYPE_SPHERE ;
214 } else if( strncasecmp(cpt,"RECT",4) == 0 ){
215 sscanf( cpt+5 , "%f,%f,%f" , &na,&nb,&nc ) ;
216 if( na == 0.0f && nb == 0.0f && nc == 0.0f )
217 ERROR_exit("'RECT(0,0,0)' is not a legal neighborhood") ;
218 ntype = NTYPE_RECT ;
219 } else if( strncasecmp(cpt,"RHDD",4) == 0 ){
220 sscanf( cpt+5 , "%f" , &na ) ;
221 if( na == 0.0f ) ERROR_exit("Can't have a RHDD of radius 0") ;
222 ntype = NTYPE_RHDD ;
223 } else if( strncasecmp(cpt,"TOHD",4) == 0 ){
224 sscanf( cpt+5 , "%f" , &na ) ;
225 if( na == 0.0f ) ERROR_exit("Can't have a TOHD of radius 0") ;
226 ntype = NTYPE_TOHD ;
227 } else {
228 ERROR_exit("Unknown -nbhd shape: '%s'",cpt) ;
229 }
230 iarg++ ; continue ;
231 }
232
233 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
234
235 } /*--- end of loop over options ---*/
236
237 if( use_nonmask ){
238 if( do_vproj ){
239 do_vproj = 0 ; WARNING_message("-use_nonmask disables -vproj") ;
240 }
241 }
242
243 /*---- deal with input dataset ----*/
244
245 if( inset == NULL ){
246 if( iarg >= argc ) ERROR_exit("No input dataset on command line?") ;
247 inset = THD_open_dataset( argv[iarg] ) ;
248 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
249 }
250 if( !IS_REAL_TYPE(DSET_BRICK_TYPE(inset,0)) )
251 ERROR_exit("Can only process real-valued datasets in this program!") ;
252 nt = DSET_NVALS(inset) ;
253 if( nt < 9 )
254 ERROR_exit("Must have at least 9 values per voxel in time series dataset '%s'",
255 DSET_BRIKNAME(inset) ) ;
256 if( polort+1 >= nt )
257 ERROR_exit("'-polort %d' too big for time series length = %d",polort,nt) ;
258
259 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
260
261 /*-- deal with mask --*/
262
263 nx = DSET_NX(inset) ;
264 ny = DSET_NY(inset) ; nxy = nx*ny ;
265 nz = DSET_NZ(inset) ; nxyz = nxy*nz ;
266
267 if( mask != NULL ){
268 if( mask_nx != DSET_NX(inset) ||
269 mask_ny != DSET_NY(inset) ||
270 mask_nz != DSET_NZ(inset) )
271 ERROR_exit("-mask dataset grid doesn't match input dataset") ;
272
273 } else if( automask ){
274 mask = THD_automask( inset ) ;
275 if( mask == NULL )
276 ERROR_message("Can't create -automask from input dataset?") ;
277 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
278 INFO_message("Number of voxels in automask = %d",nmask) ;
279 if( nmask < 2 ) ERROR_exit("Automask is too small to process") ;
280
281 } else {
282 nmask = nxyz ; /* all voxels */
283 }
284
285 /*---- create neighborhood (as a cluster) -----*/
286
287 if( ntype <= 0 ){ /* default neighborhood */
288 ntype = NTYPE_SPHERE ; na = -1.01f ;
289 INFO_message("Using default neighborhood = self + 6 neighbors") ;
290 }
291
292 switch( ntype ){
293 default:
294 ERROR_exit("WTF? ntype=%d",ntype) ;
295
296 case NTYPE_SPHERE:{
297 if( na < 0.0f ){ dx = dy = dz = 1.0f ; na = -na ; }
298 else { dx = fabsf(DSET_DX(inset)) ;
299 dy = fabsf(DSET_DY(inset)) ;
300 dz = fabsf(DSET_DZ(inset)) ; }
301 nbhd = MCW_spheremask( dx,dy,dz , na ) ;
302 }
303 break ;
304
305 case NTYPE_RECT:{
306 if( na < 0.0f ){ dx = 1.0f; na = -na; } else dx = fabsf(DSET_DX(inset));
307 if( nb < 0.0f ){ dy = 1.0f; nb = -nb; } else dy = fabsf(DSET_DY(inset));
308 if( nc < 0.0f ){ dz = 1.0f; nc = -nc; } else dz = fabsf(DSET_DZ(inset));
309 nbhd = MCW_rectmask( dx,dy,dz , na,nb,nc ) ;
310 }
311 break ;
312
313 case NTYPE_RHDD:{
314 if( na < 0.0f ){ dx = dy = dz = 1.0f ; na = -na ; }
315 else { dx = fabsf(DSET_DX(inset)) ;
316 dy = fabsf(DSET_DY(inset)) ;
317 dz = fabsf(DSET_DZ(inset)) ; }
318 nbhd = MCW_rhddmask( dx,dy,dz , na ) ;
319 }
320 break ;
321
322 case NTYPE_TOHD:{
323 if( na < 0.0f ){ dx = dy = dz = 1.0f ; na = -na ; }
324 else { dx = fabsf(DSET_DX(inset)) ;
325 dy = fabsf(DSET_DY(inset)) ;
326 dz = fabsf(DSET_DZ(inset)) ; }
327 nbhd = MCW_tohdmask( dx,dy,dz , na ) ;
328 }
329 break ;
330 }
331 MCW_radsort_cluster( nbhd , dx,dy,dz ) ; /* 26 Feb 2008 */
332 /* ensures first value is centroid */
333
334 INFO_message("Neighborhood comprises %d voxels",nbhd->num_pt) ;
335 INFO_message("Each time series has %d points",nt) ;
336 #if 0
337 printf("Neighborhood offsets:\n") ;
338 for( kk=0 ; kk < nbhd->num_pt ; kk++ )
339 printf(" [%2d] = %2d %2d %2d\n",kk,nbhd->i[kk],nbhd->j[kk],nbhd->k[kk]) ;
340 #endif
341
342 /** create output dataset **/
343
344 if( prefix != NULL && strcmp(prefix,"NULL") != 0 ){
345 outset = EDIT_empty_copy(inset) ;
346 EDIT_dset_items( outset, ADN_prefix,prefix, ADN_brick_fac,NULL, ADN_none );
347 tross_Copy_History( inset , outset ) ;
348 tross_Make_History( "3dLocalPV" , argc,argv , outset ) ;
349 for( kk=0 ; kk < nt ; kk++ ) /* create bricks */
350 EDIT_substitute_brick( outset , kk , MRI_float , NULL ) ;
351 }
352
353 if( prefix2 != NULL && strcmp(prefix2,"NULL") != 0 ){ /* 27 Apr 2016 */
354 outset2 = EDIT_empty_copy(inset) ;
355 EDIT_dset_items( outset2, ADN_prefix,prefix2, ADN_brick_fac,NULL, ADN_none );
356 tross_Copy_History( inset , outset2 ) ;
357 tross_Make_History( "3dLocalPV" , argc,argv , outset2 ) ;
358 for( kk=0 ; kk < nt ; kk++ ) /* create bricks */
359 EDIT_substitute_brick( outset2 , kk , MRI_float , NULL ) ;
360 }
361
362 if( evprefix != NULL ){
363 evset = EDIT_empty_copy(inset) ;
364 kk = (do_vproj <= 1) ? 1 : 2 ;
365 EDIT_dset_items( evset , ADN_prefix,evprefix, ADN_ntt,0 ,
366 ADN_brick_fac,NULL, ADN_nvals,kk, ADN_none ) ;
367 tross_Copy_History( inset , evset ) ;
368 tross_Make_History( "3dLocalPV" , argc,argv , evset ) ;
369 EDIT_substitute_brick( evset , 0 , MRI_float , NULL ) ;
370 evar = DSET_ARRAY(evset,0) ;
371 if( kk == 2 ){
372 EDIT_substitute_brick( evset , 1 , MRI_float , NULL ) ;
373 fvar = DSET_ARRAY(evset,1) ;
374 }
375 }
376
377 if( despike ){ /* 14 Oct 2010 */
378 THD_3dim_dataset *qset ;
379 INFO_message("Despiking input dataset") ;
380 qset = THD_despike9_dataset( inset , mask ) ;
381 if( qset == NULL ) ERROR_exit("Despiking fails!?") ;
382 DSET_delete(inset) ; inset = qset ;
383 }
384
385 if( polort >= 0 ){
386 float **polyref ; THD_3dim_dataset *qset ;
387 INFO_message("Detrending input dataset with %d basis functions",polort+1) ;
388 polyref = THD_build_polyref( polort+1 , nt ) ;
389 qset = THD_detrend_dataset( inset , polort+1,polyref , 2,0,mask,NULL ) ;
390 if( qset == NULL ) ERROR_exit("Detrending fails!?") ;
391 free(polyref) ; DSET_delete(inset) ; inset = qset ;
392 }
393
394 #ifndef USE_OMP
395 vstep = (verb && nxyz > 999) ? nxyz/50 : 0 ;
396 if( vstep ) fprintf(stderr,"++ voxel loop: ") ;
397 #endif
398
399 /*** the real work now begins ***/
400
401 AFNI_OMP_START ;
402 #pragma omp parallel if( nmask > 666 )
403 { int kk , xx,yy,zz , vv,ii , mm ;
404 float *zar=NULL , *nbar ; int *ivar ;
405 float *tsar , *uvec , *vvec=NULL , *ws=NULL ;
406 unsigned short xran[3] ; int ithr=0 ; float sval,tval=0.0 ;
407
408 #ifdef USE_OMP
409 ithr = omp_get_thread_num() ;
410 if( ithr == 0 )
411 INFO_message("Starting voxel loop: OpenMP threads=%d",omp_get_num_threads());
412 #endif
413
414 #pragma omp critical (MALLOC)
415 { uvec = (float *)malloc(sizeof(float)*nt) ;
416 ws = pv_get_workspace(nt,nbhd->num_pt) ;
417 nbar = (float *)malloc(sizeof(float)*nt*nbhd->num_pt) ;
418 ivar = (int *)malloc(sizeof(int)*nbhd->num_pt) ;
419 if( do_vproj == 2 || outset2 != NULL )
420 vvec = (float *)malloc(sizeof(float)*nt) ;
421 if( do_vproj )
422 zar = (float *)malloc(sizeof(float)*nt) ;
423 }
424
425 xran[2] = ( gseed & 0xffff) + (unsigned short)ithr ;
426 xran[1] = ((gseed >> 16) & 0xffff) - (unsigned short)ithr ;
427 xran[0] = 0x330e + (unsigned short)ithr ;
428
429 #pragma omp for
430 for( kk=0 ; kk < nxyz ; kk++ ){
431
432 #if 0
433 #pragma omp critical
434 { if( vstep && kk%vstep==vstep-1 ) vstep_print() ; }
435 #endif
436
437 if( !INMASK(kk) ) continue ;
438 IJK_TO_THREE( kk , xx,yy,zz , nx,nxy ) ;
439
440 mm = myTHD_extract_nbhd( inset , mask , xx,yy,zz , nbhd , ivar,nbar ) ;
441 if( mm <= 0 ) continue ; /* no data? */
442 if( do_vproj ){
443 for( ii=0 ; ii < nt ; ii++ ) zar[ii] = nbar[ii] ;
444 }
445 if( do_vnorm ){
446 register float sum ;
447 for( vv=0 ; vv < mm ; vv++ ){
448 tsar = nbar + nt*vv ;
449 for( sum=0.0f,ii=0 ; ii < nt ; ii++ ) sum += tsar[ii]*tsar[ii] ;
450 if( sum > 0.0f ){
451 for( sum=1.0f/sqrtf(sum),ii=0 ; ii < nt ; ii++ ) tsar[ii] *= sum ;
452 }
453 }
454 }
455 tsar = nbar ;
456 if( do_vproj <= 1 && outset2 == NULL ){
457 sval = principal_vector( nt, mm, 0,nbar, uvec, (do_vproj) ? NULL : tsar, ws,xran ) ;
458 } else {
459 float_pair spair ;
460 spair = principal_vector_pair( nt,mm, 0,nbar, uvec,vvec, NULL, ws,xran ) ;
461 sval = spair.a ; tval = spair.b ;
462 }
463 switch( do_vproj ){
464 case 1:{
465 register float sum ;
466 for( sum=0.0f,ii=0 ; ii < nt ; ii++ ) sum += zar[ii]*uvec[ii] ;
467 for( ii=0 ; ii < nt ; ii++ ) uvec[ii] *= sum ;
468 }
469 break ;
470
471 case 2:{
472 register float sum ;
473 for( sum=0.0f,ii=0 ; ii < nt ; ii++ ) sum += zar[ii]*uvec[ii] ;
474 for( ii=0 ; ii < nt ; ii++ ) uvec[ii] *= sum ;
475 for( sum=0.0f,ii=0 ; ii < nt ; ii++ ) sum += zar[ii]*vvec[ii] ;
476 for( ii=0 ; ii < nt ; ii++ ) uvec[ii] += sum*vvec[ii] ;
477 }
478 }
479 if( outset != NULL )
480 THD_insert_series( kk, outset , nt, MRI_float, uvec, 0 ) ;
481 if( outset2 != NULL )
482 THD_insert_series( kk, outset2, nt, MRI_float, vvec, 0 ) ;
483 if( evar != NULL ) evar[kk] = sval ;
484 if( fvar != NULL ) fvar[kk] = tval ;
485 }
486
487 #pragma omp critical (MALLOC)
488 { free(uvec) ; free(nbar) ; free(ws) ; free(ivar) ;
489 if( vvec != NULL ) free(vvec) ;
490 if( zar != NULL ) free(zar) ; }
491 } /* end OpenMP */
492 AFNI_OMP_END ;
493
494 #ifndef USE_OMP
495 if( vstep ) fprintf(stderr,"\n") ;
496 #else
497 ININFO_message("Voxel loop finished!") ;
498 #endif
499
500 /*** cleanup and exit ***/
501
502 DSET_delete(inset) ;
503
504 if( outset != NULL ){
505 DSET_write(outset) ; WROTE_DSET(outset) ; DSET_delete(outset) ;
506 }
507 if( outset2 != NULL ){
508 DSET_write(outset2) ; WROTE_DSET(outset2) ; DSET_delete(outset2) ;
509 }
510
511 if( evset != NULL ){
512 DSET_write(evset) ; WROTE_DSET(evset) ; DSET_delete(evset) ;
513 }
514
515 exit(0) ;
516 }
517
518 /*------------------------------------------------------------------------*/
519
520 #undef TS
521 #define TS(i,j) tsar[(i)+(j)*nv]
522
myTHD_extract_nbhd(THD_3dim_dataset * dset,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,int * ivar,float * tsar)523 static int myTHD_extract_nbhd( THD_3dim_dataset *dset, byte *mask,
524 int xx, int yy, int zz, MCW_cluster *nbhd,
525 int *ivar , float *tsar )
526 {
527 int nv , ival , kk ;
528 int nvox , nx,ny,nz,nxy,nxyz , npt , aa,bb,cc , ii,vv ;
529
530 nx = DSET_NX(dset) ;
531 ny = DSET_NY(dset) ; nxy = nx*ny ;
532 nz = DSET_NZ(dset) ; nxyz = nxy*nz ; npt = nbhd->num_pt ;
533 nv = dset->dblk->nvals ;
534
535 for( nvox=ii=0 ; ii < npt ; ii++ ){
536 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
537 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
538 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
539 vv = aa + bb*nx + cc*nxy ;
540 if( mask == NULL || mask[vv] ) ivar[nvox++] = vv ;
541 }
542 if( nvox == 0 ) return 0 ; /* nothing to extract */
543
544 /* fill the output */
545
546 switch( DSET_BRICK_TYPE(dset,0) ){
547
548 default: /* don't know what to do --> return nada */
549 return 0 ;
550
551 case MRI_byte:{
552 byte *bar ;
553 for( ival=0 ; ival < nv ; ival++ ){
554 bar = (byte *)DSET_ARRAY(dset,ival) ;
555 for( kk=0 ; kk < nvox ; kk++ ) TS(ival,kk) = (float)bar[ivar[kk]] ;
556 }
557 }
558 break ;
559
560 case MRI_short:{
561 short *bar ;
562 for( ival=0 ; ival < nv ; ival++ ){
563 bar = (short *)DSET_ARRAY(dset,ival) ;
564 for( kk=0 ; kk < nvox ; kk++ ) TS(ival,kk) = (float)bar[ivar[kk]] ;
565 }
566 }
567 break ;
568
569 case MRI_float:{
570 float *bar ;
571 for( ival=0 ; ival < nv ; ival++ ){
572 bar = (float *) DSET_ARRAY(dset,ival) ;
573 for( kk=0 ; kk < nvox ; kk++ ) TS(ival,kk) = bar[ivar[kk]] ;
574 }
575 }
576 break ;
577 }
578
579 /* scale outputs, if needed */
580
581 if( THD_need_brick_factor(dset) ){
582 float fac ;
583 for( ival=0 ; ival < nv ; ival++ ){
584 fac = DSET_BRICK_FACTOR(dset,ival) ;
585 for( kk=0 ; kk < nvox ; kk++ ) TS(ival,kk) *= fac ;
586 }
587 }
588
589 return nvox ;
590 }
591