1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 Program to resample a 'data parent' dataset to the grid defined by an
10 'anat parent' dataset.
11
12 File: adwarp.c
13 Author: B. Douglas Ward
14 Date: 02 April 1999
15
16 Mod: Added changes for incorporating History notes.
17 Date: 10 September 1999
18
19 Mod: Added changes for proper byte ordering on output.
20 Date: 23 Nov 1999 - RW Cox
21
22 Mod: Added -force option, and some checks related to it.
23 Date: 13 Dec 1999 - RW Cox
24
25 Mod: Added 'data parent' sub-brick selector. Also, allow operator to
26 specify separate resampling modes for the functional sub-bricks
27 and the threshold sub-bricks.
28 Date: 25 February 2000
29
30 Mod: Added call to AFNI_logger.
31 Date: 15 August 2001
32
33 */
34
35 /*---------------------------------------------------------------------------*/
36
37 #define PROGRAM_NAME "adwarp.c" /* name of this program */
38 #define PROGRAM_AUTHOR "R. W. Cox and B. D. Ward" /* program author */
39 #define PROGRAM_INITIAL "02 April 1999" /* date of initial program release */
40 #define PROGRAM_LATEST "15 August 2001" /* date of last program revision */
41
42 /*---------------------------------------------------------------------------*/
43 #include "afni.h"
44 #include "afni_warp.h"
45
46 #define MAIN
47
48 /*---------------------------------------------------------------------------*/
49
AFNI_copy_statistics(THD_3dim_dataset * dsold,THD_3dim_dataset * dsnew)50 void AFNI_copy_statistics( THD_3dim_dataset *dsold , THD_3dim_dataset *dsnew )
51 {
52 int ibr , nvold , nvnew ;
53 THD_statistics *stold , *stnew ;
54
55 ENTRY("AFNI_copy_statistics") ;
56
57 if( !ISVALID_3DIM_DATASET(dsold) || !ISVALID_3DIM_DATASET(dsnew) ) EXRETURN ;
58
59 nvold = dsold->dblk->nvals ;
60 nvnew = dsnew->dblk->nvals ;
61 stold = dsold->stats ;
62 stnew = dsnew->stats ;
63 if( !ISVALID_STATISTIC(stold) ) EXRETURN ;
64
65 if( stnew == NULL ){
66 dsnew->stats = stnew = myXtNew( THD_statistics ) ;
67 stnew->type = STATISTICS_TYPE ;
68 stnew->nbstat = nvnew ;
69 stnew->bstat = (THD_brick_stats *)
70 XtMalloc( sizeof(THD_brick_stats) * nvnew ) ;
71 ADDTO_KILL(dsnew->kl,stnew) ;
72 stnew->parent = (XtPointer) dsnew ;
73 } else {
74 stnew->nbstat = nvnew ;
75 stnew->bstat = (THD_brick_stats *)
76 XtRealloc( (char *) stnew->bstat ,
77 sizeof(THD_brick_stats) * nvnew ) ;
78 }
79
80 for( ibr=0 ; ibr < nvnew ; ibr++ ){
81 if( ibr < nvold )
82 stnew->bstat[ibr] = stold->bstat[ibr] ;
83 else
84 INVALIDATE_BSTAT(stnew->bstat[ibr]) ;
85 }
86
87 EXRETURN ;
88 }
89
90 /*---------------------------------------------------------------------------*/
91
92 #define PARENTIZE(ds,par) \
93 if( ISVALID_3DIM_DATASET((ds)) ) (ds)->parent = (XtPointer) (par)
94
95 /** macro to test a malloc-ed pointer for validity **/
96 #define MTEST(ptr) \
97 if((ptr)==NULL) \
98 ( AW_error ("Cannot allocate memory") )
99
100 /*---------------------------------------------------------------------------*/
101
102 typedef struct adwarp_options
103 {
104 THD_3dim_dataset *aset; /* anat parent dataset */
105 THD_3dim_dataset *dset; /* data parent dataset */
106
107 char *prefix; /* prefix for output file */
108
109 float dxyz; /* grid spacing in output dataset (mm) */
110
111 int anat_resam_mode; /* anatomical dataset resampling mode */
112 int thr_resam_mode; /* threshold sub-brick resampling mode */
113 int func_resam_mode; /* functional sub-brick resampling mode */
114
115 int verbose ; /* RWCox: 06 Apr 1999 */
116 int force ; /* RWCox: 13 Dec 1999 */
117
118 } adwarp_options;
119
120
121 /*---------------------------------------------------------------------------*/
122 /*
123 Print error message and stop.
124 */
125
AW_error(char * message)126 void AW_error (char *message)
127 {
128 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
129 exit(1);
130 }
131
132
133 /*---------------------------------------------------------------------------*/
134 /*
135 Routine to display adwarp help menu.
136 */
137
display_help_menu()138 void display_help_menu()
139 {
140 int ii ;
141
142
143 printf
144 ("Usage: adwarp [options]\n"
145 "Resamples a 'data parent' dataset to the grid defined by an\n"
146 "'anat parent' dataset. The anat parent dataset must contain\n"
147 "in its .HEAD file the coordinate transformation (warp) needed\n"
148 "to bring the data parent dataset to the output grid. This\n"
149 "program provides a batch implementation of the interactive\n"
150 "AFNI 'Write' buttons, one dataset at a time.\n"
151 "\n"
152 " Example: create dataset func+tlrc (.HEAD and .BRIK) by applying\n"
153 " the orig->tlrc transformation from the anat.\n"
154 "\n"
155 " adwarp -apar anat+tlrc -dpar func+orig\n"
156 "\n"
157 " Example: in the case of a manual tlrc transformation, maybe the\n"
158 " anat+tlrc.BRIK does not exist (just the .HEAD file does).\n"
159 " In such a case on might apply the anat+tlrc transformation\n"
160 " to the anat+orig dataset. But since the anat+tlrc.HEAD\n"
161 " file already exists, the -overwrite option is needed.\n"
162 "\n"
163 " adwarp -apar anat+tlrc -dpar anat+orig -overwrite\n"
164 "\n"
165 "Options (so to speak):\n"
166 "----------------------\n"
167 "-apar aset = Set the anat parent dataset to 'aset'. This\n"
168 " is a nonoptional option (must be present).\n"
169 "\n"
170 "-dpar dset = Set the data parent dataset to 'dset'. This\n"
171 " is a nonoptional option (must be present).\n"
172 " Note: dset may contain a sub-brick selector,\n"
173 " e.g., -dpar 'dset+orig[2,5,7]' \n"
174 "\n"
175 "-prefix ppp = Set the prefix for the output dataset to 'ppp'.\n"
176 " The default is the prefix of 'dset'.\n"
177 "\n"
178 "-dxyz ddd = Set the grid spacing in the output datset to\n"
179 " 'ddd' mm. The default is 1 mm.\n"
180 "\n"
181 "-verbose = Print out progress reports.\n"
182 "-force = Write out result even if it means deleting\n"
183 " an existing dataset. The default is not\n"
184 " to overwrite.\n"
185 "\n"
186 "-resam rrr = Set resampling mode to 'rrr' for all sub-bricks\n"
187 " --- OR --- \n"
188 "-thr rrr = Set resampling mode to 'rrr' for threshold sub-bricks\n"
189 "-func rrr = Set resampling mode to 'rrr' for functional sub-bricks\n"
190 "\n"
191 "The resampling mode 'rrr' must be one of the following:\n"
192 ) ;
193
194 for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
195 printf(" %s = %s\n", RESAM_shortstr[ii],RESAM_typestr[ii] ) ;
196
197
198 printf
199 (
200 "\n"
201 "NOTE: The default resampling mode is %s for all sub-bricks. \n" ,
202 RESAM_shortstr[1]
203 ) ;
204
205 PRINT_COMPILE_DATE ; exit(0) ;
206 }
207
208
209 /*---------------------------------------------------------------------------*/
210 /*
211 Routine to initialize the input options.
212 */
213
214
initialize_options(adwarp_options * option_data)215 void initialize_options (adwarp_options *option_data)
216 {
217 option_data->aset = NULL;
218 option_data->dset = NULL;
219 option_data->prefix = NULL;
220
221 option_data->dxyz = 1.0;
222
223 option_data->anat_resam_mode = 1;
224 option_data->thr_resam_mode = 1;
225 option_data->func_resam_mode = 1;
226
227 option_data->verbose = 0 ; /* 06 Apr 1999 */
228 option_data->force = 0 ; /* 13 Dec 1999 */
229 }
230
231
232 /*---------------------------------------------------------------------------*/
233 /*
234 Routine to get user specified input options.
235 */
236
get_options(int argc,char ** argv,adwarp_options * option_data)237 void get_options
238 (
239 int argc, /* number of input arguments */
240 char ** argv, /* array of input arguments */
241 adwarp_options *option_data /* adwarp program options */
242 )
243
244 {
245 int nopt = 1; /* input option argument counter */
246 int ival; /* integer input */
247 float fval; /* float input */
248 int ii; /* index */
249 char message[THD_MAX_NAME]; /* error message */
250
251
252 /*----- does user request help menu? -----*/
253 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
254
255 AFNI_logger(PROGRAM_NAME,argc,argv) ;
256
257
258 /*----- initialize the input options -----*/
259 initialize_options (option_data);
260
261
262 /*----- main loop over input options -----*/
263 while (nopt < argc )
264 {
265
266 /*----- -apar aset -----*/
267 if (strncmp(argv[nopt], "-apar", 5) == 0)
268 {
269 nopt++;
270 if (nopt >= argc) AW_error ("need argument after -apar ");
271 option_data->aset = THD_open_one_dataset (argv[nopt]);
272 if (! ISVALID_3DIM_DATASET(option_data->aset))
273 AW_error ("Cannot read anat parent dataset.\n") ;
274 THD_make_cardinal(option_data->aset); /* 21 Oct, 2011 [rickr] */
275
276 nopt++;
277 continue;
278 }
279
280 /*----- -dpar dset -----*/
281 if (strncmp(argv[nopt], "-dpar", 5) == 0)
282 {
283 nopt++;
284 if (nopt >= argc) AW_error ("need argument after -dpar ");
285 option_data->dset = THD_open_dataset (argv[nopt]);
286 if (! ISVALID_3DIM_DATASET(option_data->dset))
287 AW_error ("Cannot read data parent dataset.\n") ;
288 THD_make_cardinal(option_data->dset); /* 21 Oct, 2011 [rickr] */
289 nopt++;
290 continue;
291 }
292
293 /*----- -prefix ppp -----*/
294 if (strncmp(argv[nopt], "-prefix", 7) == 0)
295 {
296 nopt++;
297 if (nopt >= argc) AW_error ("need argument after -prefix ");
298 option_data->prefix = AFMALL(char,sizeof(char)*THD_MAX_PREFIX);
299 MTEST (option_data->prefix);
300 MCW_strncpy (option_data->prefix, argv[nopt], THD_MAX_PREFIX);
301
302 if( strstr(option_data->prefix,".nii") != NULL || /* 06 Apr 2005 */
303 strstr(option_data->prefix,".hdr") != NULL ){ /* 11 Oct 2005 */
304 fprintf(stderr,"** You can't use adwarp to create a NIfTI file!\n") ;
305 exit(1) ;
306 }
307
308 nopt++;
309 continue;
310 }
311
312 /*----- -verbose -----*/
313 if( strncmp(argv[nopt],"-verbose",5) == 0 ){ /* 06 Apr 1999 */
314 option_data->verbose = 1 ;
315 nopt++ ; continue ;
316 }
317
318 /*----- -force -----*/
319 if( strncmp(argv[nopt],"-force",5) == 0 ){ /* 06 Apr 1999 */
320 option_data->force = 1 ;
321 nopt++ ; continue ;
322 }
323
324 /*----- -dxyz ddd -----*/
325 if (strncmp(argv[nopt], "-dxyz", 5) == 0)
326 {
327 nopt++;
328 if (nopt >= argc) AW_error ("need argument after -dxyz ");
329 sscanf (argv[nopt], "%f", &fval);
330 if (fval <= 0.0) AW_error ("Illegal argument after -dxyz ");
331 option_data->dxyz = fval;
332 nopt++;
333 continue;
334 }
335
336 /*----- -resam rrr -----*/
337 if (strncmp(argv[nopt], "-resam", 6) == 0)
338 {
339 nopt++;
340 if (nopt >= argc) AW_error ("need argument after -resam ");
341 ival = -1;
342 for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
343 if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
344 ival = ii;
345 if ((ival < 0) || (ival > LAST_RESAM_TYPE))
346 AW_error ("illegal argument after -resam ");
347 option_data->anat_resam_mode = ival;
348 option_data->thr_resam_mode = ival;
349 option_data->func_resam_mode = ival;
350 nopt++;
351 continue;
352 }
353
354
355 /*----- -thr rrr -----*/
356 if (strncmp(argv[nopt], "-thr", 4) == 0)
357 {
358 nopt++;
359 if (nopt >= argc) AW_error ("need argument after -thr ");
360 ival = -1;
361 for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
362 if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
363 ival = ii;
364 if ((ival < 0) || (ival > LAST_RESAM_TYPE))
365 AW_error ("illegal argument after -thr ");
366 option_data->thr_resam_mode = ival;
367 nopt++;
368 continue;
369 }
370
371
372 /*----- -func rrr -----*/
373 if (strncmp(argv[nopt], "-func", 5) == 0)
374 {
375 nopt++;
376 if (nopt >= argc) AW_error ("need argument after -func ");
377 ival = -1;
378 for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
379 if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
380 ival = ii;
381 if ((ival < 0) || (ival > LAST_RESAM_TYPE))
382 AW_error ("illegal argument after -func ");
383 option_data->func_resam_mode = ival;
384 nopt++;
385 continue;
386 }
387
388
389 /*----- unknown command -----*/
390 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
391 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
392 suggest_best_prog_option(argv[0], argv[nopt]);
393 exit(1);
394
395
396 }
397
398
399 /*----- Check for required inputs -----*/
400 if (option_data->aset == NULL)
401 AW_error ("Must specify anat parent dataset");
402 if (option_data->dset == NULL)
403 AW_error ("Must specify data parent dataset");
404
405 /*-- 13 Dec 1999: check if datasets have the same view --*/
406
407 if( option_data->aset->view_type == option_data->dset->view_type ){
408 if( option_data->force ){
409 if( option_data->prefix == NULL ){
410 fprintf(stderr,
411 "** Error: -apar & -dpar are in same +view!\n"
412 " This is illegal without the -prefix option!\n") ;
413 exit(1) ;
414 } else {
415 fprintf(stderr,
416 "++ Warning: -apar & -dpar are in same +view!\n") ;
417 }
418 } else {
419 fprintf(stderr,
420 "** Error: -apar & -dpar are in same +view!\n"
421 " If this is OK, use -force and -prefix options.\n" ) ;
422 exit(1) ;
423 }
424 }
425
426
427 }
428
429
430 /*---------------------------------------------------------------------------*/
431 /*
432 Make a warped dataset whose grid corresponds to the anat_parent and
433 whose data comes from the data_parent.
434 Note that the assumption is made that data_parent and the warp parent
435 of the anat_parent are both in the same coordinate system (up to the
436 to_dicomm transformation of their dataxes structs).
437
438 This routine is adapted from AFNI_follower_dataset.
439 */
440
adwarp_follower_dataset(adwarp_options * option_data,THD_3dim_dataset * anat_parent,THD_3dim_dataset * data_parent)441 THD_3dim_dataset * adwarp_follower_dataset
442 (
443 adwarp_options *option_data, /* adwarp program options */
444 THD_3dim_dataset *anat_parent, /* dataset containing warp information */
445 THD_3dim_dataset *data_parent /* dataset to be warped */
446 )
447 {
448 THD_3dim_dataset *new_dset ;
449 int ii ;
450
451 ENTRY("adwarp_follower_dataset") ;
452
453 /* sanity checks */
454
455 if( ! ISVALID_3DIM_DATASET(anat_parent) ||
456 ! ISVALID_3DIM_DATASET(data_parent) ) RETURN(NULL) ;
457
458 /* make new dataset, copying appropriate fields from its various parents */
459
460 new_dset = myXtNew( THD_3dim_dataset ) ; INIT_KILL( new_dset->kl ) ;
461
462 new_dset->type = data_parent->type; /* same data type */
463 new_dset->func_type = data_parent->func_type;
464 new_dset->view_type = anat_parent->view_type; /* but different view type */
465 /* use template space of parent to mark as TLRC/MNI/... */
466 MCW_strncpy( new_dset->atlas_space ,
467 anat_parent->atlas_space , THD_MAX_NAME ) ;
468 new_dset->anat_parent = anat_parent; /* what else makes sense? */
469
470 new_dset->tagset = NULL ; /* Oct 1998 */
471 new_dset->Label_Dtable = NULL; /* ZSS Feb 26 2010 */
472
473 MCW_strncpy( new_dset->anat_parent_name ,
474 anat_parent->self_name , THD_MAX_NAME ) ;
475
476 new_dset->anat_parent_idcode = anat_parent->idcode ;
477
478 /* 11/09/94 addition: the data_parent may itself be a warp;
479 in this case, we want the true warp parent to be the original data */
480
481 new_dset->warp_parent = (data_parent->warp_parent != NULL)
482 ? (data_parent->warp_parent) : (data_parent) ;
483
484 MCW_strncpy( new_dset->warp_parent_name ,
485 new_dset->warp_parent->self_name , THD_MAX_NAME ) ;
486
487 new_dset->warp_parent_idcode = new_dset->warp_parent->idcode ;
488
489 new_dset->idcode = MCW_new_idcode() ;
490
491 /* make the actual warp from the warp_parent to this dataset */
492
493 new_dset->vox_warp = NULL ;
494 new_dset->warp = myXtNew( THD_warp ) ;
495 *(new_dset->warp) = IDENTITY_WARP ; /* start with (Dicom) identity */
496
497 new_dset->self_warp = NULL ; /* 26 Aug 2002 */
498
499 /* follow the links backward from desired view to original view */
500
501 AFNI_concatenate_warp( new_dset->warp , anat_parent->warp ) ;
502 AFNI_concatenate_warp( new_dset->warp , data_parent->warp ) ;
503
504 /* reset the bounds in the new warp to be the same as in the anat_parent */
505
506 if( ISVALID_WARP(anat_parent->warp) &&
507 anat_parent->warp->type == new_dset->warp->type ){
508
509 switch( anat_parent->warp->type ){
510
511 case WARP_AFFINE_TYPE:
512 COPY_LMAP_BOUNDS( new_dset->warp->rig_bod.warp ,
513 anat_parent->warp->rig_bod.warp ) ;
514 break ;
515
516 case WARP_TALAIRACH_12_TYPE:
517 for( ii=0 ; ii < 12 ; ii++ )
518 COPY_LMAP_BOUNDS( new_dset->warp->tal_12.warp[ii] ,
519 anat_parent->warp->tal_12.warp[ii] ) ;
520 break ;
521 }
522 }
523
524 /* make up some names for this new dataset */
525
526 MCW_strncpy( new_dset->self_name ,
527 new_dset->warp_parent->self_name , THD_MAX_NAME ) ;
528 ii = strlen( new_dset->self_name ) ;
529 new_dset->self_name[ii++] = '@' ;
530 MCW_strncpy( &(new_dset->self_name[ii]) ,
531 VIEW_typestr[new_dset->view_type] , THD_MAX_NAME-ii ) ;
532
533 MCW_strncpy( new_dset->label1 , data_parent->label1 , THD_MAX_LABEL ) ;
534 MCW_strncpy( new_dset->label2 , data_parent->label2 , THD_MAX_LABEL ) ;
535
536 /* set the axes for this new dataset
537 (same as anatomy parent, since that's the meaning of this routine) */
538
539 new_dset->daxes = myXtNew( THD_dataxes ) ; /* copy data axes of */
540 *(new_dset->daxes) = *(anat_parent->daxes) ; /* anatomy parent */
541
542 new_dset->wod_daxes = NULL ;
543 new_dset->wod_flag = True ;
544
545 /* 06 Aug 1996: added ability to use 3D+t datasets here */
546
547 if( DSET_NUM_TIMES(data_parent) < 2 ){
548 new_dset->taxis = NULL ;
549 } else {
550 new_dset->taxis = myXtNew( THD_timeaxis ) ; /* new */
551 *(new_dset->taxis) = *(data_parent->taxis) ; /* copy insides */
552
553 new_dset->taxis->nsl = 0 ; /* no slice stuff */
554 new_dset->taxis->toff_sl = NULL ;
555 new_dset->taxis->zorg_sl = 0.0 ;
556 new_dset->taxis->dz_sl = 0.0 ;
557 }
558
559 /* create a datablock and diskptr, in case the data is ever
560 filled into memory (instead of wod) and written to disk */
561
562 new_dset->dblk = myXtNew( THD_datablock ) ; INIT_KILL( new_dset->dblk->kl ) ;
563
564 new_dset->dblk->type = DATABLOCK_TYPE ;
565 new_dset->dblk->nvals = data_parent->dblk->nvals ;
566 new_dset->dblk->malloc_type = DATABLOCK_MEM_UNDEFINED ;
567 new_dset->dblk->natr = new_dset->dblk->natr_alloc = 0 ;
568 new_dset->dblk->atr = NULL ;
569 new_dset->dblk->parent = (XtPointer) new_dset ;
570
571 if( data_parent->dblk->brick_lab == NULL ){
572 THD_init_datablock_labels( new_dset->dblk ) ; /* 30 Nov 1997 */
573 } else {
574 THD_copy_datablock_auxdata( data_parent->dblk , new_dset->dblk ) ;
575 }
576
577 DSET_unlock(new_dset) ; /* Feb 1998 */
578
579 new_dset->dblk->diskptr = myXtNew( THD_diskptr ) ;
580 new_dset->dblk->diskptr->type = DISKPTR_TYPE ;
581 new_dset->dblk->diskptr->nvals = data_parent->dblk->nvals ;
582 new_dset->dblk->diskptr->rank = 3 ;
583 new_dset->dblk->diskptr->storage_mode = STORAGE_UNDEFINED ;
584 new_dset->dblk->diskptr->byte_order = THD_get_write_order() ;
585 /* 25 April 1998 */
586 new_dset->dblk->diskptr->dimsizes[0] = new_dset->daxes->nxx ;
587 new_dset->dblk->diskptr->dimsizes[1] = new_dset->daxes->nyy ;
588 new_dset->dblk->diskptr->dimsizes[2] = new_dset->daxes->nzz ;
589
590 new_dset->dblk->brick_fac = NULL ; /* initialized below */
591 new_dset->dblk->brick_bytes = NULL ;
592 new_dset->dblk->brick = NULL ;
593 THD_init_datablock_brick( new_dset->dblk , -1 , data_parent->dblk ) ;
594
595 new_dset->dblk->master_nvals = 0 ; /* 11 Jan 1999 */
596 new_dset->dblk->master_ival = NULL ;
597 new_dset->dblk->master_bytes = NULL ;
598
599 /* create the names for storage on disk (if ever)
600 -- note we put it in the same directory as the data_parent */
601
602 if (option_data->prefix == NULL)
603 THD_init_diskptr_names (new_dset->dblk->diskptr,
604 data_parent->dblk->diskptr->directory_name, NULL,
605 data_parent->dblk->diskptr->prefix,
606 new_dset->view_type, True);
607 else
608 THD_init_diskptr_names (new_dset->dblk->diskptr,
609 data_parent->dblk->diskptr->directory_name, NULL,
610 option_data->prefix,
611 new_dset->view_type, True);
612
613
614 ADDTO_KILL( new_dset->dblk->kl , new_dset->dblk->diskptr ) ;
615
616 /* oh yeah, set the new_dset kill list,
617 copy statistics if available, and NULL out any unused stuff */
618
619 ADDTO_KILL( new_dset->kl , new_dset->warp ) ;
620 ADDTO_KILL( new_dset->kl , new_dset->daxes ) ;
621 ADDTO_KILL( new_dset->kl , new_dset->dblk ) ;
622
623 new_dset->stats = NULL ;
624 AFNI_copy_statistics( data_parent , new_dset ) ;
625
626 INIT_STAT_AUX( new_dset , MAX_STAT_AUX , data_parent->stat_aux ) ;
627
628 new_dset->markers = NULL ; /* no markers */
629 new_dset->death_mark = 0 ; /* don't kill me! */
630 #ifdef ALLOW_DATASET_VLIST
631 new_dset->pts = NULL ;
632 #endif
633
634 PARENTIZE(new_dset,data_parent->parent) ;
635
636 new_dset->tcat_list = NULL ; /* 03 Aug 2004 */
637 new_dset->tcat_num = 0 ;
638 new_dset->tcat_len = NULL ;
639
640 RETURN(new_dset) ;
641 }
642
643
644 /*---------------------------------------------------------------------------*/
645 /*
646 (re)compute and (re)write a dataset to disk, with the indicated
647 geometric parameters. Note that the filenames for the save should
648 already be initialized in the dset->dblk->diskptr. Also, note
649 that the dataset's "permanent" dataxes will be remade to fit the
650 new geometry.
651
652 This routine is adapted from AFNI_refashion_dataset.
653 */
654
655
adwarp_refashion_dataset(adwarp_options * option_data,THD_3dim_dataset * dset,THD_dataxes * daxes)656 RwcBoolean adwarp_refashion_dataset
657 (
658 adwarp_options *option_data, /* adwarp program options */
659 THD_3dim_dataset *dset, /* new (output) dataset */
660 THD_dataxes *daxes /* new dataset axes */
661 )
662 {
663 THD_datablock *dblk = dset->dblk ;
664 THD_diskptr *dkptr = dset->dblk->diskptr ;
665 RwcBoolean good ;
666 int npix , nx,ny,nz,nv , kk , ival , code , nzv , dsiz , isfunc , cmode ;
667 MRI_IMAGE *im ;
668 void *imar ;
669 FILE *far ;
670 float brfac_save ;
671 int resam_mode = 0;
672
673 int native_order , save_order ; /* 23 Nov 1999 */
674
675 ENTRY("adwarp_refashion_dataset") ;
676
677 /* set up for warp-on-demand */
678
679 dset->wod_daxes = myXtNew(THD_dataxes) ; /* 02 Nov 1996 */
680 dset->wod_daxes->type = DATAXES_TYPE ; /* 02 Nov 1996 */
681 dset->vox_warp = myXtNew(THD_warp) ; /* 02 Nov 1996 */
682
683 dset->self_warp = NULL ; /* 26 Aug 2002 */
684
685 *(dset->wod_daxes) = *daxes ; /* copy insides of daxes */
686 dset->wod_flag = True ; /* mark for warp-on-demand */
687 dset->vox_warp->type = ILLEGAL_TYPE ; /* mark for recomputation */
688
689 /* copy the new geometric information into various places */
690
691 *(dset->daxes) = *(daxes) ; /* Make daxes permanent */
692 dkptr->dimsizes[0] = dset->daxes->nxx ; /* Will cause trouble */
693 dkptr->dimsizes[1] = dset->daxes->nyy ; /* if diskptr and */
694 dkptr->dimsizes[2] = dset->daxes->nzz ; /* daxes don't match! */
695
696 /* write the header out */
697
698 good = THD_write_3dim_dataset( NULL,NULL , dset , False ) ;
699 if( !good ){
700 fprintf(stderr,"\a\n*** cannot write dataset header ***\n") ;
701
702 RETURN(False) ;
703 }
704 STATUS("wrote output header file") ;
705
706 /* at this point in time, the output dataset is set and the header
707 * has been written - so we must overwrite from now on */
708 putenv("AFNI_DECONFLICT=OVERWRITE") ; /* 19 Nov 2007 */
709
710 /* purge the datablock that now exists,
711 then delete the file on disk that now exists (if any) */
712
713 DSET_unlock( dset ) ; /* Feb 1998 */
714 PURGE_DSET( dset ) ;
715 COMPRESS_unlink(dkptr->brick_name) ;
716
717 /* refashion its brick data structure,
718 which requires first saving in a temporary
719 array the datum type for each sub-brick */
720
721 { int ibr ; int *typ ;
722 typ = (int *) XtMalloc( sizeof(int) * dblk->nvals ) ;
723 for( ibr=0 ; ibr < dblk->nvals ; ibr++ )
724 typ[ibr] = DBLK_BRICK_TYPE(dblk,ibr) ;
725 THD_init_datablock_brick( dblk , dblk->nvals , typ ) ;
726 myXtFree( typ ) ;
727 }
728
729 /*-- 13 Mar 2006: check for free disk space --*/
730
731 { int mm = THD_freemegabytes(dkptr->header_name) ;
732 int rr = (int)(dblk->total_bytes/(1024*1024)) ;
733 if( rr >= 666 )
734 fprintf(stderr,"++ WARNING: output filesize %s will be %d Mbytes!\n"
735 "++ SUGGEST: increase voxel size to save disk space.\n",
736 dkptr->brick_name , rr ) ;
737 if( mm >= 0 && mm <= rr )
738 WARNING_message("Disk space: writing file %s (%d MB),"
739 " but only %d free MB on disk" ,
740 dkptr->brick_name , rr , mm ) ;
741 }
742
743
744 dkptr->storage_mode = STORAGE_UNDEFINED ; /* just for now */
745 dblk->malloc_type = DATABLOCK_MEM_UNDEFINED ;
746
747 /*--- open the output file
748 (N.B.: much of the following code is from THD_write_datablock) ---*/
749
750 /*-- create directory if necessary --*/
751
752 if( ! THD_is_directory(dkptr->directory_name) ){
753 kk = mkdir( dkptr->directory_name , THD_MKDIR_MODE ) ;
754 if( kk != 0 ){
755 fprintf(stderr,
756 "\a\n*** cannot mkdir new directory: %s\n",dkptr->directory_name) ;
757
758 RETURN(False) ;
759 }
760 STATUS("created subdirectory") ;
761 }
762
763 /*-- open output file --*/
764
765 if( option_data->verbose )
766 fprintf(stderr,"++ Opening output file %s\n",dkptr->brick_name) ;
767
768 cmode = THD_get_write_compression() ;
769 far = COMPRESS_fopen_write( dkptr->brick_name , cmode ) ;
770 if( far == NULL ){
771 fprintf(stderr,
772 "\a\n*** cannot open output file %s\n",dkptr->brick_name) ;
773
774 RETURN(False) ;
775 }
776 STATUS("created output brick file") ;
777
778 /*--------- now, create each slice and write it out ----------*/
779
780 nx = dset->daxes->nxx ;
781 ny = dset->daxes->nyy ; npix = nx*ny ;
782 nz = dset->daxes->nzz ;
783 nv = dkptr->nvals ; nzv = nz*nv ;
784
785 isfunc = ISFUNC(dset) ; /* 09 Dec 1997 */
786
787 if( ! isfunc )
788 resam_mode = option_data->anat_resam_mode ;
789
790 native_order = mri_short_order() ; /* 23 Nov 1999 */
791 save_order = (dkptr->byte_order > 0) ? dkptr->byte_order
792 : THD_get_write_order() ;
793
794 for( ival=0 ; ival < nv ; ival++ ){ /* for each sub-brick */
795
796 if( option_data->verbose )
797 fprintf(stderr,"++ Start sub-brick %d",ival) ;
798
799 dsiz = mri_datum_size( DSET_BRICK_TYPE(dset,ival) ) ;
800
801 /** force return of unscaled slices for output **/
802
803 brfac_save = DBLK_BRICK_FACTOR(dblk,ival) ;
804 DBLK_BRICK_FACTOR(dblk,ival) = 0.0 ;
805
806 if( isfunc )
807 resam_mode = (DSET_BRICK_STATCODE(dset,ival) > 0) /* 09 Dec 1997 */
808 ? option_data->thr_resam_mode
809 : option_data->func_resam_mode ;
810
811 for( kk=0 ; kk < nz ; kk++ ){ /* for each slice */
812
813 im = AFNI_dataset_slice( dset , 3 , kk , ival , resam_mode ) ;
814 STATUS("have new image") ;
815
816 if( option_data->verbose && kk%7==0 ) fprintf(stderr,".");
817
818 if( im == NULL ){
819 fprintf(stderr,"\a\n*** failure to compute dataset slice %d\n",kk) ;
820 COMPRESS_fclose(far) ;
821 COMPRESS_unlink( dkptr->brick_name ) ;
822
823 RETURN(False) ;
824 }
825
826 #ifdef AFNI_DEBUG
827 { char str[256] ;
828 sprintf(str,"writing slice %d: type=%s nx=%d ny=%d\n",
829 kk,MRI_TYPE_NAME(im) , im->nx,im->ny ) ;
830 STATUS(str) ; }
831 #endif
832
833 imar = mri_data_pointer(im) ;
834 if( save_order != native_order ){ /* 23 Nov 1999 */
835 switch( im->kind ){
836 default: break ;
837 case MRI_short: mri_swap2( npix,imar) ; break ;
838 case MRI_float:
839 case MRI_int: mri_swap4( npix,imar) ; break ;
840 case MRI_complex: mri_swap4(2*npix,imar) ; break ;
841 }
842 }
843 code = fwrite( imar , dsiz , npix , far ) ;
844 mri_free(im) ;
845
846 if( code != npix ){
847 fprintf(stderr,
848 "\a\n*** failure to write dataset slice %d (is disk full?)\n",kk) ;
849 COMPRESS_fclose(far) ;
850 COMPRESS_unlink( dkptr->brick_name ) ;
851
852 RETURN(False) ;
853 }
854
855 } /* end of loop over kk (z-direction) */
856
857 if( option_data->verbose ) fprintf(stderr,"\n");
858
859 /* restore the correct scaling of this sub-brick */
860
861 DBLK_BRICK_FACTOR(dblk,ival) = brfac_save ;
862
863 } /* end of loop over iv (nvals direction) */
864 STATUS("all slices written") ;
865
866 /*--------------------- done!!! ---------------------*/
867
868 COMPRESS_fclose(far) ;
869 STATUS("output file closed") ;
870
871 /*--- do a little surgery on the dataset's storage flags ---*/
872
873 dkptr->storage_mode = STORAGE_BY_BRICK ;
874 #if MMAP_THRESHOLD > 0
875 dblk->malloc_type = (dblk->total_bytes > MMAP_THRESHOLD)
876 ? DATABLOCK_MEM_MMAP : DATABLOCK_MEM_MALLOC ;
877
878 if( cmode >= 0 ) dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
879 DBLK_mmapfix(dblk) ; /* 28 Mar 2005 */
880 #else
881 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
882 #endif
883
884 /*--- recompute the statistics and rewrite the header to hold them ---*/
885
886 STATUS("recomputing statistics") ;
887
888 if( option_data->verbose ) fprintf(stderr,"++ Computing statistics\n");
889 THD_load_statistics( dset ) ;
890
891 STATUS("rewriting header") ;
892 (void) THD_write_3dim_dataset( NULL,NULL , dset , False ) ;
893
894 STATUS("purging datablock") ;
895
896 PURGE_DSET( dset ) ;
897
898 myXtFree(dset->wod_daxes) ; myXtFree(dset->vox_warp) ; /* 02 Nov 1996 */
899
900 RETURN(True) ;
901 }
902
903
904 /*---------------------------------------------------------------------------*/
905
main(int argc,char * argv[])906 int main( int argc , char *argv[] )
907 {
908 adwarp_options *option_data; /* adwarp program options */
909 THD_3dim_dataset *new_dset = NULL; /* new (output) dataset */
910 THD_dataxes new_daxes; /* new dataset axes */
911
912 mainENTRY("adwarp main") ; machdep() ; PRINT_VERSION("adwarp") ;
913 AUTHOR(PROGRAM_AUTHOR) ;
914
915 #if 0
916 /*----- Identify software -----*/
917 printf ("\n\n");
918 printf ("Program: %s \n", PROGRAM_NAME);
919 printf ("Author: %s \n", PROGRAM_AUTHOR);
920 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
921 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
922 printf ("\n");
923 #endif
924
925 /*----- Allocate memory -----*/
926 option_data = (adwarp_options *) malloc (sizeof(adwarp_options));
927
928 /*----- Get user inputs -----*/
929 get_options (argc, argv, option_data);
930
931 /*----- Create empty shell of dataset from anat parent and data parent ----*/
932 new_dset = adwarp_follower_dataset (option_data, option_data->aset,
933 option_data->dset);
934
935 /*--- 13 Dec 1999: check if output dataset already exists on disk ---*/
936
937 if( THD_is_file(DSET_HEADNAME(new_dset)) ||
938 THD_is_file(DSET_BRIKNAME(new_dset)) ){
939
940 if( option_data->force ){
941 fprintf(stderr,
942 "++ Warning: overwriting dataset %s and %s\n",
943 DSET_HEADNAME(new_dset), DSET_BRIKNAME(new_dset) ) ;
944 putenv("AFNI_DECONFLICT=OVERWRITE") ; /* 12 Nov 2007 */
945 }
946 }
947
948 /* aside from -force, let the user's AFNI_DECONFLICT variable control this */
949 /* 19 Nov 2007 [rickr] */
950 #if 0
951 else if( THD_deathcon() ){
952 fprintf(stderr,
953 "** Error: can't overwrite dataset %s and %s\n"
954 " unless you use the -force option!\n" ,
955 DSET_HEADNAME(new_dset), DSET_BRIKNAME(new_dset) ) ;
956 exit(1) ;
957 } else {
958 putenv("AFNI_DECONFLICT=YES") ; /* 12 Nov 2007 */
959 THD_deconflict_prefix(new_dset) ;
960 WARNING_message("Changed dataset name to '%s' to avoid conflict",
961 DSET_BRIKNAME(new_dset) ) ;
962 }
963 #endif
964
965 /*----- Record history of dataset -----*/
966 tross_Copy_History( option_data->dset , new_dset ) ;
967 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
968
969 /*----- Allow for resampling to a new voxel size -----*/
970 new_daxes.type = DATAXES_TYPE;
971 THD_edit_dataxes (option_data->dxyz, option_data->aset->daxes, &new_daxes);
972
973
974 /*----- Fill in the dataset and write out to disk -----*/
975 adwarp_refashion_dataset (option_data, new_dset, &new_daxes);
976
977 exit(0) ;
978 }
979