1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*-- 14 Oct 1999: modified to allow for inverse warping --*/
10
11 /*-- 18 Oct 1999: modified to implement -preserve option --*/
12
13 THD_fvec3 AFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv ) ;
14 /*== Note any changes here should be reflected in both ================
15 plug_drawdset.c:DRAW_ttatlas_CB()
16 3dfractionize.c:main()
17 or updated to use the same function instead in the future....
18 Unless voting, this function doesn't really need to work this way;
19 it can instead use the 3dresample or future warp-on-demand way.
20 Resulting ROIs are different between this method and "Show Atlas Colors"
21 ==========================================================================
22 */
main(int argc,char * argv[])23 int main( int argc , char * argv[] )
24 {
25 char *prefix = "fractionize" ;
26 THD_3dim_dataset *tset=NULL , *iset=NULL , *dset=NULL , *wset=NULL ;
27 int iarg=1 ;
28 int nxin,nyin,nzin , nxyin ;
29 float dxin,dyin,dzin , xorgin,yorgin,zorgin , clip=0.0 ;
30 int nxout,nyout,nzout , nxyout , nvoxout ;
31 float dxout,dyout,dzout , xorgout,yorgout,zorgout ;
32 float f1,f2,f , g1,g2,g , h1,h2,h , sx,sy,sz , tx,ty,tz ;
33 float x1,x2 , y1,y2 , z1,z2 , xx1,xx2 , yy1,yy2 , zz1,zz2 ;
34 int i,ip , j,jp , k,kp , iv,jv,kv , vtype , ijk ;
35 float *voxout ;
36 byte *bin = NULL ;
37 short *sin = NULL , sclip=0 ;
38 float *fin = NULL , fclip=0.00001 ;
39 void *voxin = NULL ;
40 THD_fvec3 vv ;
41
42 THD_warp *warp=NULL ; /* 14 Oct 1999 */
43
44 int do_vote=0 ; /* 18 Oct 1999 */
45 int *vote_val = NULL ;
46 int nvote_val=0 , ivote , voter=0 , vote_print=0 ;
47 byte *vote_bout = NULL ;
48 short *vote_sout = NULL ;
49 float *vote_best = NULL ;
50
51 /*-- help? --*/
52
53 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
54 printf("Usage: 3dfractionize [options]\n"
55 "\n"
56 "* For each voxel in the output dataset, computes the fraction\n"
57 " of it that is occupied by nonzero voxels from the input.\n"
58 "* The fraction is stored as a short in the range 0..10000,\n"
59 " indicating fractions running from 0..1.\n"
60 "* The template dataset is used only to define the output grid;\n"
61 " its brick(s) will not be read into memory. (The same is\n"
62 " true of the warp dataset, if it is used.)\n"
63 "* The actual values stored in the input dataset are irrelevant,\n"
64 " except in that they are zero or nonzero (UNLESS the -preserve\n"
65 " option is used).\n"
66 "\n"
67 "The purpose of this program is to allow the resampling of a mask\n"
68 "dataset (the input) from a fine grid to a coarse grid (defined by\n"
69 "the template). When you are using the output, you will probably\n"
70 "want to threshold the mask so that voxels with a tiny occupancy\n"
71 "fraction aren't used. This can be done in 3dmaskave, by using\n"
72 "3calc, or with the '-clip' option below.\n"
73 "\n"
74 "Options are [the first 2 are 'mandatory options']:\n"
75 " -template tset = Use dataset 'tset' as a template for the output.\n"
76 " The output dataset will be on the same grid as\n"
77 " this dataset.\n"
78 "\n"
79 " -input iset = Use dataset 'iset' for the input.\n"
80 " Only the sub-brick #0 of the input is used.\n"
81 " You can use the sub-brick selection technique\n"
82 " described in '3dcalc -help' to choose the\n"
83 " desired sub-brick from a multi-brick dataset.\n"
84 "\n"
85 " -prefix ppp = Use 'ppp' for the prefix of the output.\n"
86 " [default prefix = 'fractionize']\n"
87 "\n"
88 " -clip fff = Clip off voxels that are less than 'fff' occupied.\n"
89 " 'fff' can be a number between 0.0 and 1.0, meaning\n"
90 " the fraction occupied, can be a number between 1.0\n"
91 " and 100.0, meaning the percent occupied, or can be\n"
92 " a number between 100.0 and 10000.0, meaning the\n"
93 " direct output value to use as a clip level.\n"
94 " ** Some sort of clipping is desirable; otherwise,\n"
95 " an output voxel that is barely overlapped by a\n"
96 " single nonzero input voxel will enter the mask.\n"
97 " [default clip = 0.0]\n"
98 "\n"
99 " -warp wset = If this option is used, 'wset' is a dataset that\n"
100 " provides a transformation (warp) from +orig\n"
101 " coordinates to the coordinates of 'iset'.\n"
102 " In this case, the output dataset will be in\n"
103 " +orig coordinates rather than the coordinates\n"
104 " of 'iset'. With this option:\n"
105 " ** 'tset' must be in +orig coordinates\n"
106 " ** 'iset' must be in +acpc or +tlrc coordinates\n"
107 " ** 'wset' must be in the same coordinates as 'iset'\n"
108 "\n"
109 " -preserve = When this option is used, the program will copy\n"
110 " or the nonzero values of input voxels to the output\n"
111 " -vote dataset, rather than create a fractional mask.\n"
112 " Since each output voxel might be overlapped\n"
113 " by more than one input voxel, the program 'votes'\n"
114 " for which input value to preserve. For example,\n"
115 " if input voxels with value=1 occupy 10%% of an\n"
116 " output voxel, and inputs with value=2 occupy 20%%\n"
117 " of the same voxel, then the output value in that\n"
118 " voxel will be set to 2 (provided that 20%% is >=\n"
119 " to the clip fraction).\n"
120 " ** Voting can only be done on short-valued datasets,\n"
121 " or on byte-valued datasets.\n"
122 " ** Voting is a relatively time-consuming option,\n"
123 " since a separate loop is made through the\n"
124 " input dataset for each distinct value found.\n"
125 " ** Combining this with the -warp option does NOT\n"
126 " make a general +tlrc to +orig transformer!\n"
127 " This is because for any value to survive the\n"
128 " vote, its fraction in the output voxel must be\n"
129 " >= clip fraction, regardless of other values\n"
130 " present in the output voxel.\n"
131 "\n"
132 "Sample usage:\n"
133 "\n"
134 " 1. Compute the fraction of each voxel occupied by the warped input.\n"
135 "\n"
136 " 3dfractionize -template grid+orig -input data+tlrc \\\n"
137 " -warp anat+tlrc -clip 0.2\n"
138 "\n"
139 " 2. Apply the (inverse) -warp transformation to transform the -input\n"
140 " from +tlrc space to +orig space, storing it according to the grid\n"
141 " of the -template.\n"
142 " A voxel in the output dataset gets the value that occupies most of\n"
143 " its volume, providing that value occupies 20%% of the voxel.\n"
144 "\n"
145 " Note that the essential difference from above is '-preserve'.\n"
146 "\n"
147 " 3dfractionize -template grid+orig -input data+tlrc \\\n"
148 " -warp anat+tlrc -preserve -clip 0.2 \\\n"
149 " -prefix new_data\n"
150 "\n"
151 " Note that 3dAllineate can also be used to warp from +tlrc to +orig\n"
152 " space. In this case, data is computed through interpolation, rather\n"
153 " than voting based on the fraction of a voxel occupied by each data\n"
154 " value. The transformation comes from the WARP_DATA attribute directly.\n"
155 " Nearest neighbor interpolation is used in this 'mask' example.\n"
156 "\n"
157 " cat_matvec -ONELINE anat+tlrc::WARP_DATA > tlrc.aff12.1D\n"
158 " 3dAllineate -1Dmatrix_apply tlrc.aff12.1D -source group_mask+tlrc \\\n"
159 " -master subj_epi+orig -prefix subj_mask -final NN\n"
160 "\n"
161 "This program will also work in going from a coarse grid to a fine grid,\n"
162 "but it isn't clear that this capability has any purpose.\n"
163 "-- RWCox - February 1999\n"
164 " - October 1999: added -warp and -preserve options\n"
165 ) ;
166 PRINT_COMPILE_DATE ; exit(0) ;
167 }
168
169 mainENTRY("3dfractionize main"); machdep();
170 AFNI_logger("3dfractionize",argc,argv);
171 PRINT_VERSION("3dfractionize"); AUTHOR("RW Cox");
172
173 /*-- read command line args --*/
174
175 while( iarg < argc && argv[iarg][0] == '-' ){
176
177 if( strcmp(argv[iarg],"-clip") == 0 ){
178 clip = strtod( argv[++iarg] , NULL ) ;
179
180 if( clip <= 1.0 ) sclip = (short)(10000.0 * clip + 0.49999) ;
181 else if( clip <= 100.0 ) sclip = (short)( 100.0 * clip + 0.49999) ;
182 else if( clip <= 10000.0 ) sclip = (short)(clip) ;
183 else sclip = -1 ;
184
185 if( sclip < 0 || sclip > 10000 ){
186 fprintf(stderr,"** Illegal clip value!\n") ; exit(1) ;
187 }
188
189 fclip = 0.0001 * sclip ;
190 iarg++ ; continue ;
191 }
192
193 if( strcmp(argv[iarg],"-template") == 0 || strcmp(argv[iarg],"-tset") == 0 ){
194 if( tset != NULL ){
195 fprintf(stderr,"** Can't have more than one -template argument!\n") ;
196 exit(1) ;
197 }
198 tset = THD_open_one_dataset( argv[++iarg] ) ;
199 if( tset == NULL ){
200 fprintf(stderr,"** Can't open template %s\n",argv[iarg]) ;
201 exit(1) ;
202 }
203 iarg++ ; continue ;
204 THD_make_cardinal(tset); /* deoblique 21 Oct, 2011 [rickr] */
205 }
206
207 if( strcmp(argv[iarg],"-input") == 0 || strcmp(argv[iarg],"-iset") == 0 ){
208 if( iset != NULL ){
209 fprintf(stderr,"** Can't have more than one -input argument!\n") ;
210 exit(1) ;
211 }
212 iset = THD_open_dataset( argv[++iarg] ) ;
213 if( iset == NULL ){
214 fprintf(stderr,"** Can't open input %s\n",argv[iarg]) ;
215 exit(1) ;
216 }
217 THD_make_cardinal(iset); /* deoblique 21 Oct, 2011 [rickr] */
218 iarg++ ; continue ;
219 }
220
221 if( strcmp(argv[iarg],"-warp") == 0 || strcmp(argv[iarg],"-wset") == 0 ){
222 if( wset != NULL ){
223 fprintf(stderr,"** Can't have more than one -warp argument!\n") ;
224 exit(1) ;
225 }
226 wset = THD_open_dataset( argv[++iarg] ) ;
227 if( wset == NULL ){
228 fprintf(stderr,"** Can't open warp %s\n",argv[iarg]) ;
229 exit(1) ;
230 }
231 THD_make_cardinal(wset); /* deoblique 21 Oct, 2011 [rickr] */
232 iarg++ ; continue ;
233 }
234
235 if( strcmp(argv[iarg],"-prefix") == 0 ){
236 prefix = argv[++iarg] ;
237 iarg++ ; continue ;
238 }
239
240 if( strcmp(argv[iarg],"-preserve") == 0 || strcmp(argv[iarg],"-vote") == 0 ){
241 do_vote = 1 ;
242 iarg++ ; continue ;
243 }
244
245 fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
246 }
247
248 /*-- check inputs for sanity --*/
249
250 if( tset == NULL ){
251 fprintf(stderr,"** No template dataset?\n") ; exit(1) ;
252 }
253 if( iset == NULL ){
254 fprintf(stderr,"** No input dataset?\n") ; exit(1) ;
255 }
256 if( !THD_filename_ok(prefix) ){
257 fprintf(stderr,"** Illegal prefix?\n") ; exit(1) ;
258 }
259
260 /*- 14 Oct 1999: additional checks with the -warp option -*/
261
262 if( wset != NULL ){
263
264 if( tset->view_type != VIEW_ORIGINAL_TYPE ){
265 fprintf(stderr,"** Template is not in +orig view - this is illegal!\n");
266 exit(1);
267 }
268
269 if( iset->view_type == VIEW_ORIGINAL_TYPE ){
270 fprintf(stderr,"** Input is in +orig view - this is illegal!\n"); exit(1);
271 }
272
273 if( wset != NULL && wset->view_type != iset->view_type ){
274 fprintf(stderr,"** Warp and Input are not in same view - this is illegal!\n");
275 exit(1);
276 }
277
278 warp = wset->warp ;
279 if( warp == NULL ){
280 fprintf(stderr,"** Warp dataset doesn't contain a warp transformation!\n");
281 exit(1) ;
282 }
283 }
284
285 /*- 18 Oct 1999: check input for legality if we are voting --*/
286
287 if( do_vote ){
288 vtype = DSET_BRICK_TYPE(iset,0) ;
289 if( vtype != MRI_short && vtype != MRI_byte ){
290 fprintf(stderr,"** -preserve option requires short- or byte-valued input dataset!\n");
291 exit(1);
292 }
293 }
294
295 /*-- start to create output dataset --*/
296
297 dset = EDIT_empty_copy( tset ) ;
298
299 tross_Copy_History( iset , dset ) ; /* 14 Oct 1999 */
300 tross_Make_History( "3dfractionize" , argc,argv , dset ) ;
301
302 EDIT_dset_items( dset ,
303 ADN_prefix , prefix ,
304 ADN_nvals , 1 ,
305 ADN_ntt , 0 ,
306 ADN_brick_fac , NULL ,
307 ADN_datum_all , MRI_short ,
308 ADN_none ) ;
309
310 if( ISFUNC(dset) )
311 EDIT_dset_items( dset , ADN_func_type,FUNC_FIM_TYPE , ADN_none ) ;
312
313 if( THD_deathcon() && THD_is_file(dset->dblk->diskptr->header_name) ){
314 fprintf(stderr,
315 "** Output file %s already exists -- cannot continue!\n",
316 dset->dblk->diskptr->header_name ) ;
317 exit(1) ;
318 }
319
320 DSET_delete( tset ) ; /* don't need template any more */
321
322 /* load the input dataset */
323
324 DSET_load( iset ) ; CHECK_LOAD_ERROR(iset) ;
325 voxin = DSET_ARRAY(iset,0) ; vtype = DSET_BRICK_TYPE(iset,0) ;
326 switch( vtype ){
327 default: fprintf(stderr,"** Illegal brick type in input dataset!\n"); exit(1);
328
329 case MRI_byte: bin = (byte * ) voxin ; break ;
330 case MRI_short: sin = (short *) voxin ; break ;
331 case MRI_float: fin = (float *) voxin ; break ;
332 }
333
334 /*-- setup voxel index (iv,jv,kv) to coordinate (x,y,z) calculations --*/
335
336 nxin =iset->daxes->nxx ; nyin =iset->daxes->nyy ; nzin =iset->daxes->nzz ;
337 dxin =iset->daxes->xxdel; dyin =iset->daxes->yydel; dzin =iset->daxes->zzdel;
338 xorgin =iset->daxes->xxorg; yorgin =iset->daxes->yyorg; zorgin =iset->daxes->zzorg;
339
340 nxout =dset->daxes->nxx ; nyout =dset->daxes->nyy ; nzout =dset->daxes->nzz ;
341 dxout =dset->daxes->xxdel; dyout =dset->daxes->yydel; dzout =dset->daxes->zzdel;
342 xorgout=dset->daxes->xxorg; yorgout=dset->daxes->yyorg; zorgout=dset->daxes->zzorg;
343
344 /* voxel (i,j,k) center is at (xorg+i*dx,yorg+j*dy,zorg+k*dz) */
345
346 nxyout = nxout * nyout ; nvoxout = nxyout * nzout ;
347
348 voxout = (float *) malloc( sizeof(float) * nvoxout ) ;
349
350 nxyin = nxin * nyin ;
351
352 /*-- if voting, do some setup --*/
353
354 if( do_vote ){
355 int nvoxin = nxyin * nzin ;
356
357 /* 1: find all distinct nonzero values in dataset */
358
359 vote_val = (int *) malloc( sizeof(int) ) ;
360 nvote_val = 0 ;
361 for( iv=0 ; iv < nvoxin ; iv++ ){
362 kv = (vtype == MRI_short) ? sin[iv] : bin[iv] ; /* voxel value */
363 if( kv == 0 ) continue ; /* skip zeroes */
364 for( jv=0 ; jv < nvote_val ; jv++ ) /* find in list */
365 if( vote_val[jv] == kv ) break ;
366 if( jv < nvote_val ) continue ; /* skip if already in list */
367
368 vote_val = (int *) realloc( vote_val , sizeof(int)*(nvote_val+1) ) ;
369 vote_val[nvote_val++] = kv ;
370 }
371 if( nvote_val == 0 ){
372 fprintf(stderr,"** Input dataset is all zero!\n") ; exit(1) ;
373 }
374 fprintf(stderr,"++ Found %d distinct nonzero values in input.\n",nvote_val) ;
375
376 if( nvote_val > 1 ) /* arrange into ascending value */
377 qsort_int( nvote_val , vote_val ) ;
378
379 /* 2: make a receptacle for the voting results */
380
381 if( vtype == MRI_byte ){
382 vote_bout = (byte *) malloc( sizeof(byte) * nvoxout ) ;
383 memset( vote_bout , 0 , sizeof(byte) * nvoxout ) ;
384 } else {
385 vote_sout = (short *) malloc( sizeof(short) * nvoxout ) ;
386 memset( vote_sout , 0 , sizeof(short) * nvoxout ) ;
387 }
388
389 /* 3: make a place to hold the best fraction yet */
390
391 vote_best = (float *) malloc( sizeof(float) * nvoxout ) ;
392 for( i=0 ; i < nvoxout ; i++ ) vote_best[i] = 0.0 ;
393
394 /* 4: if needed, attach the correct brick factor to the output */
395
396 f = DSET_BRICK_FACTOR(iset,0) ;
397 if( f != 0.0 && f != 1.0 ) EDIT_BRICK_FACTOR(dset,0,f) ;
398
399 /* 5: if input is bytes, make the output be that too */
400
401 if( vtype == MRI_byte )
402 EDIT_dset_items( dset , ADN_datum_all , MRI_byte , ADN_none ) ;
403 }
404
405 /*-- loop over values to vote on --*/
406
407 ivote = 0 ;
408 do{
409
410 if( do_vote ){
411 voter = vote_val[ivote] ; /* who gets to vote this time */
412 if( ivote%10 == 9 ){
413 if( !vote_print ) fprintf(stderr,"++ Voting:") ;
414 fprintf(stderr,"%d",(ivote/10)%10) ;
415 vote_print++ ;
416 }
417 }
418
419 for( i=0 ; i < nvoxout ; i++ ) voxout[i] = 0.0 ; /* fractions */
420
421 /*-- loop over input voxels --*/
422
423 for( kv=0 ; kv < nzin ; kv++ ){
424
425 z1 = zorgin + dzin * (kv-0.5) ; z2 = zorgin + dzin * (kv+0.49999) ;
426
427 for( jv=0 ; jv < nyin ; jv++ ){
428
429 y1 = yorgin + dyin * (jv-0.5) ; y2 = yorgin + dyin * (jv+0.49999) ;
430
431 for( iv=0 ; iv < nxin ; iv++ ){
432
433 ijk = iv + jv*nxin + kv*nxyin ; /* 1D index of voxel (iv,jv,kv) */
434
435 /* decide if we use this voxel, based on its value */
436
437 if( do_vote ){
438 if( vtype == MRI_short && sin[ijk] != voter ) continue ;
439 else if( vtype == MRI_byte && bin[ijk] != voter ) continue ;
440 } else {
441 if( vtype == MRI_short && sin[ijk] == 0 ) continue ;
442 else if( vtype == MRI_byte && bin[ijk] == 0 ) continue ;
443 else if( vtype == MRI_float && fin[ijk] == 0 ) continue ;
444 }
445
446 x1 = xorgin + dxin * (iv-0.5) ; x2 = xorgin + dxin * (iv+0.49999) ;
447
448 /* input voxel (iv,jv,kv) spans coordinates [x1,x2] X [y1,y2] X [z1,z2] */
449
450 /* transform these corner coordinates to output dataset grid coordinates */
451
452 LOAD_FVEC3(vv , x1,y1,z1) ;
453 vv = THD_3dmm_to_dicomm( iset , vv ); /* transpose from iset to Dicom */
454 vv = AFNI_backward_warp_vector( warp , vv ); /* transform to +orig ???*/
455 vv = THD_dicomm_to_3dmm( dset , vv ); /* and back from Dicom to dset */
456 UNLOAD_FVEC3(vv , xx1,yy1,zz1) ;
457
458 LOAD_FVEC3(vv , x2,y2,z2) ;
459 vv = THD_3dmm_to_dicomm( iset , vv ) ;
460 vv = AFNI_backward_warp_vector( warp , vv ) ;
461 vv = THD_dicomm_to_3dmm( dset , vv ) ;
462 UNLOAD_FVEC3(vv , xx2,yy2,zz2) ;
463
464 /* [xx1,xx2] X [yy1,yy2] X [zz1,zz2] is now in coordinates of output dataset */
465
466 /* compute indices into output dataset voxel (keeping fractions) */
467
468 f1 = (xx1-xorgout)/dxout + 0.49999 ; f2 = (xx2-xorgout)/dxout + 0.49999 ;
469 if( f1 > f2 ){ tx = f1 ; f1 = f2 ; f2 = tx ; }
470 if( f1 >= nxout || f2 <= 0.0 ) continue ;
471 if( f1 < 0.0 ) f1 = 0.0 ; if( f2 >= nxout ) f2 = nxout - 0.001 ;
472
473 g1 = (yy1-yorgout)/dyout + 0.49999 ; g2 = (yy2-yorgout)/dyout + 0.49999 ;
474 if( g1 > g2 ){ ty = g1 ; g1 = g2 ; g2 = ty ; }
475 if( g1 >= nyout || g2 <= 0.0 ) continue ;
476 if( g1 < 0.0 ) g1 = 0.0 ; if( g2 >= nyout ) g2 = nyout - 0.001 ;
477
478 h1 = (zz1-zorgout)/dzout + 0.49999 ; h2 = (zz2-zorgout)/dzout + 0.49999 ;
479 if( h1 > h2 ){ tz = h1 ; h1 = h2 ; h2 = tz ; }
480 if( h1 >= nzout || h2 <= 0.0 ) continue ;
481 if( h1 < 0.0 ) h1 = 0.0 ; if( h2 >= nzout ) h2 = nzout - 0.001 ;
482
483 /* input voxel covers voxels [f1,f2] X [g1,g2] X [h1,h2] in the output */
484
485 /* For example, [6.3,7.2] X [9.3,9.6] X [11.7,13.4], which must be */
486 /* distributed into these voxels: */
487 /* (6,9,11), (7,9,11), (6,9,12), (7,9,12), (6,9,13), and (7,9,13) */
488
489 for( f=f1 ; f < f2 ; f = ip ){
490 i = (int) f ; ip = i+1 ; tx = MIN(ip,f2) ; sx = tx - f ;
491 for( g=g1 ; g < g2 ; g = jp ){
492 j = (int) g ; jp = j+1 ; ty = MIN(jp,g2) ; sy = ty - g ;
493 for( h=h1 ; h < h2 ; h = kp ){
494 k = (int) h ; kp = k+1 ; tz = MIN(kp,h2) ; sz = tz - h ;
495 voxout[ i + j*nxout + k * nxyout ] += sx * sy * sz ;
496
497 }
498 }
499 }
500
501 }}} /* end of loops over voxels */
502
503 /* at this point, voxout[i] contains the fraction that output voxel [i]
504 is occupied by input voxels that were allowed to vote (or were nonzero) */
505
506 /* if not voting, we are done; otherwise, we must tally the count so far */
507
508 if( do_vote ){
509 for( i=0 ; i < nvoxout ; i++ ){
510 if( voxout[i] > vote_best[i] && voxout[i] >= fclip ){ /* wins */
511 vote_best[i] = voxout[i] ;
512 if( vtype == MRI_byte ) vote_bout[i] = (byte) voter ;
513 else if( vtype == MRI_short ) vote_sout[i] = (short) voter ;
514 }
515 }
516 }
517
518 /* if we are voting, loop back for the next ballot */
519
520 } while( do_vote && ++ivote < nvote_val ) ; /* end of voting loop */
521
522 /*- throw out some trash -*/
523
524 DSET_delete(iset) ;
525 if( wset != NULL ) DSET_delete(wset) ;
526
527 if( do_vote ){
528 free(voxout) ; free(vote_best) ; free(vote_val) ;
529 if( vote_print ) fprintf(stderr,"\n") ;
530 }
531
532 /*-- if not voting, compute the output brick --*/
533
534 if( ! do_vote ){
535 sin = (short *) malloc( sizeof(short) * nvoxout ) ;
536 if( sin == NULL ){
537 fprintf(stderr,"** Can't malloc output brick!\n") ; exit(1) ;
538 }
539
540 for( i=ijk=0 ; i < nvoxout ; i++ ){
541 sin[i] = (short) (10000.0 * voxout[i] + 0.49999) ;
542 if( sin[i] < sclip ) sin[i] = 0 ;
543 else if( sin[i] > 10000 ) sin[i] = 10000 ;
544
545 if( sin[i] != 0 ) ijk++ ;
546 }
547
548 free(voxout) ;
549 EDIT_substitute_brick( dset , 0 , MRI_short , sin ) ;
550
551 fprintf(stderr,"-- Writing %d nonzero mask voxels to dataset %s\n",
552 ijk , DSET_BRIKNAME(dset) ) ;
553
554 /*-- if voting, output brick will be in vote_bout or vote_sout --*/
555
556 } else {
557 if( vtype == MRI_byte ){
558 EDIT_substitute_brick( dset , 0 , MRI_byte , vote_bout ) ;
559 for( i=ijk=0 ; i < nvoxout ; i++ )
560 if( vote_bout[i] != 0 ) ijk++ ;
561 } else {
562 EDIT_substitute_brick( dset , 0 , MRI_short , vote_sout ) ;
563 for( i=ijk=0 ; i < nvoxout ; i++ )
564 if( vote_sout[i] != 0 ) ijk++ ;
565 }
566
567 fprintf(stderr,"-- Writing %d nonzero voted voxels to dataset %s\n",
568 ijk , DSET_BRIKNAME(dset) ) ;
569 }
570
571 DSET_write(dset) ;
572 exit(0) ;
573 }
574
575 #if 0 /* Now in libmri.a ZSS Feb 2006 */
576 /*------------------------------------------------------------------------
577 Backward transform a vector following a warp
578 - lifted from afni.c
579 - note that a NULL warp is equivalent to the identity
580 --------------------------------------------------------------------------*/
581
582 THD_fvec3 AFNI_backward_warp_vector( THD_warp * warp , THD_fvec3 old_fv )
583 {
584 THD_fvec3 new_fv ;
585
586 if( warp == NULL ) return old_fv ;
587
588 switch( warp->type ){
589
590 default: new_fv = old_fv ; break ;
591
592 case WARP_TALAIRACH_12_TYPE:{
593 THD_linear_mapping map ;
594 int iw ;
595
596 /* test if input is in bot..top of each defined map */
597
598 for( iw=0 ; iw < 12 ; iw++ ){
599 map = warp->tal_12.warp[iw] ;
600
601 if( old_fv.xyz[0] >= map.bot.xyz[0] &&
602 old_fv.xyz[1] >= map.bot.xyz[1] &&
603 old_fv.xyz[2] >= map.bot.xyz[2] &&
604 old_fv.xyz[0] <= map.top.xyz[0] &&
605 old_fv.xyz[1] <= map.top.xyz[1] &&
606 old_fv.xyz[2] <= map.top.xyz[2] ) break ; /* leave loop */
607 }
608 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
609 }
610 break ;
611
612 case WARP_AFFINE_TYPE:{
613 THD_linear_mapping map = warp->rig_bod.warp ;
614 new_fv = MATVEC_SUB(map.mbac,old_fv,map.svec) ;
615 }
616 break ;
617
618 }
619 return new_fv ;
620 }
621 #endif
622