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 /*---------------------------------------------------------------------------*/
8 /*
9 This program performs 2d image registration of slices contained in an AFNI
10 3d+time dataset. This program was adapted from plug_imreg.c and imreg.c.
11
12
13 File: 2dImReg.c
14 Author: B. Douglas Ward
15 Date: 04 February 1998
16
17
18 Mod: Added routines to write the registration parameters, and the RMS
19 error, to user specified ASCII files.
20 Date: 20 March 1998
21
22 Mod: Added option to change dx and dy output format from pixels to mm.
23 Date: 24 March 1998
24
25 Mod: No longer assume that input and base datasets have same length.
26 This problem was reported by Paul Reber.
27 Date: 3 June 1998
28
29 Mod: Routine eval_registration extended to include byte and float datum
30 types for base and input datasets.
31 Date: 3 June 1998
32
33 Mod: Corrected problem with base image memory deallocation.
34 Date: 20 July 1998
35
36 Mod: Added changes for incorporating History notes.
37 Date: 10 September 1999
38
39 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
40 Date: 02 December 2002
41
42 Mod: If nsl == 0, set num_slices equal to nz. [4 occurrences]
43 Date: 06 October 2003 [rickr]
44 */
45
46 /*---------------------------------------------------------------------------*/
47
48 #define PROGRAM_NAME "2dImReg" /* name of this program */
49 #define PROGRAM_INITIAL "04 Feb 1998" /* date of initial program release */
50 #define PROGRAM_LATEST "02 Dec 2002" /* date of latest program revision */
51
52 /*---------------------------------------------------------------------------*/
53
54 #include "mrilib.h"
55 #include "matrix.h"
56
57 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
58 #define STATE_DIM 4 /* number of registration parameters */
59
60
61 /*----- Global variables -----*/
62 int * t_to_z = NULL; /* convert time-order to z-order of slices */
63 int * z_to_t = NULL; /* convert z-order to time-order of slices */
64
65
66 static char * commandline = NULL ; /* command line for history notes */
67
68
69 typedef struct IR_options /* user input options */
70 {
71 char * input_filename; /* file name for input 3d+time dataset */
72 char * base_filename; /* file name for reference (base) volume */
73 int base_vol_index; /* image number for base volume */
74 int nofine; /* boolean for no fine fit */
75 float blur; /* FWHM of blurring prior to registration */
76 float dxy; /* convergence tolerance for translations */
77 float dphi; /* convergence tolerance for rotations */
78 char * new_prefix; /* prefix name for registered dataset */
79 char * dprefix; /* prefix name for registration parameters */
80 int dmm; /* change dx and dy output from pixels to mm */
81 char * rprefix; /* prefix name for volume RMS error */
82 int debug; /* write additional output to screen */
83 } IR_options;
84
85
86
87 /*---------------------------------------------------------------------------*/
88 /*
89 Routine to display 2dImReg help menu.
90 */
91
display_help_menu()92 void display_help_menu()
93 {
94 printf
95 (
96 "This program performs 2d image registration. Image alignment is \n"
97 "performed on a slice-by-slice basis for the input 3d+time dataset, \n"
98 "relative to a user specified base image. \n"
99 "\n"
100 " ** Note that the script @2dwarper.Allin can do similar things, **\n"
101 " ** with nonlinear (polynomial) warping on a slice-wise basis. **\n"
102 " \n"
103 "Usage: \n"
104 "2dImReg \n"
105 "-input fname Filename of input 3d+time dataset to process \n"
106 "-basefile fname Filename of 3d+time dataset for base image \n"
107 " (default = current input dataset) \n"
108 "-base num Time index for base image (0 <= num) \n"
109 " (default: num = 3) \n"
110 "-nofine Deactivate fine fit phase of image registration\n"
111 " (default: fine fit is active) \n"
112 "-fine blur dxy dphi Set fine fit parameters \n"
113 " where: \n"
114 " blur = FWHM of blurring prior to registration (in pixels) \n"
115 " (default: blur = 1.0) \n"
116 " dxy = Convergence tolerance for translations (in pixels) \n"
117 " (default: dxy = 0.07) \n"
118 " dphi = Convergence tolerance for rotations (in degrees) \n"
119 " (default: dphi = 0.21) \n"
120 " \n"
121 "-prefix pname Prefix name for output 3d+time dataset \n"
122 " \n"
123 "-dprefix dname Write files 'dname'.dx, 'dname'.dy, 'dname'.psi \n"
124 " containing the registration parameters for each \n"
125 " slice in chronological order. \n"
126 " File formats: \n"
127 " 'dname'.dx: time(sec) dx(pixels) \n"
128 " 'dname'.dy: time(sec) dy(pixels) \n"
129 " 'dname'.psi: time(sec) psi(degrees) \n"
130 "-dmm Change dx and dy output format from pixels to mm \n"
131 " \n"
132 "-rprefix rname Write files 'rname'.oldrms and 'rname'.newrms \n"
133 " containing the volume RMS error for the original \n"
134 " and the registered datasets, respectively. \n"
135 " File formats: \n"
136 " 'rname'.oldrms: volume(number) rms_error \n"
137 " 'rname'.newrms: volume(number) rms_error \n"
138 " \n"
139 "-debug Lots of additional output to screen \n"
140 );
141
142 exit(0);
143 }
144
145
146 /*---------------------------------------------------------------------------*/
147 /*
148 Print error message and stop.
149 */
150
IR_error(char * message)151 void IR_error
152 (
153 char * message /* error message to be displayed */
154 )
155
156 {
157 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
158 exit(1);
159 }
160
161
162 /*---------------------------------------------------------------------------*/
163 /*
164 Routine to initialize the input options.
165 */
166
initialize_options(IR_options ** opt)167 void initialize_options
168 (
169 IR_options ** opt /* user input options */
170 )
171
172 {
173 (*opt) = (IR_options *) malloc (sizeof(IR_options));
174
175 (*opt)->input_filename = NULL;
176 (*opt)->base_filename = NULL;
177 (*opt)->base_vol_index = -3;
178 (*opt)->nofine = 0;
179 (*opt)->blur = 1.0;
180 (*opt)->dxy = 0.07;
181 (*opt)->dphi = 0.21;
182 (*opt)->new_prefix = NULL;
183 (*opt)->dprefix = NULL;
184 (*opt)->rprefix = NULL;
185 (*opt)->dmm = 0;
186 (*opt)->debug = 0;
187 }
188
189
190 /*---------------------------------------------------------------------------*/
191 /*
192 Routine to get user specified input options.
193 */
194
get_user_inputs(int argc,char ** argv,IR_options ** option_data)195 void get_user_inputs
196 (
197 int argc, /* number of input arguments */
198 char ** argv, /* array of input arguments */
199 IR_options ** option_data /* user input options */
200 )
201
202 {
203 int nopt = 1; /* input option argument counter */
204 int ival;
205 float fval;
206 char message[80]; /* error message */
207
208
209 /*----- does user request help menu? -----*/
210 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
211
212
213 /*----- main loop over input options -----*/
214 while (nopt < argc )
215 {
216
217 /*----- -input filename -----*/
218 if (strncmp(argv[nopt], "-input", 6) == 0)
219 {
220 nopt++;
221 if (nopt >= argc) IR_error ("Need argument after -input ");
222 (*option_data)->input_filename
223 = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
224 strcpy ((*option_data)->input_filename, argv[nopt]);
225 nopt++;
226 continue;
227 }
228
229
230 /*----- -basefile filename -----*/
231 if (strncmp(argv[nopt], "-basefile", 9) == 0)
232 {
233 nopt++;
234 if (nopt >= argc) IR_error ("Need argument after -basefile ");
235 (*option_data)->base_filename
236 = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
237 strcpy ((*option_data)->base_filename, argv[nopt]);
238 nopt++;
239 continue;
240 }
241
242
243 /*----- -base num -----*/
244 if (strncmp(argv[nopt], "-base", 5) == 0)
245 {
246 nopt++;
247 if (nopt >= argc) IR_error ("Need argument after -base ");
248 sscanf (argv[nopt], "%d", &ival);
249 if (ival < 0)
250 IR_error ("Illegal argument after -base ( must be >= 0 ) ");
251 (*option_data)->base_vol_index = ival;
252 nopt++;
253 continue;
254 }
255
256
257 /*----- -prefix pname -----*/
258 if (strncmp(argv[nopt], "-prefix", 7) == 0)
259 {
260 nopt++;
261 if (nopt >= argc) IR_error ("Need argument after -prefix ");
262 (*option_data)->new_prefix
263 = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
264 strcpy ((*option_data)->new_prefix, argv[nopt]);
265
266 nopt++;
267 continue;
268 }
269
270
271 /*----- -dprefix dname -----*/
272 if (strncmp(argv[nopt], "-dprefix", 8) == 0)
273 {
274 nopt++;
275 if (nopt >= argc) IR_error ("Need argument after -dprefix ");
276 (*option_data)->dprefix
277 = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
278 strcpy ((*option_data)->dprefix, argv[nopt]);
279
280 nopt++;
281 continue;
282 }
283
284
285 /*----- -rprefix rname -----*/
286 if (strncmp(argv[nopt], "-rprefix", 8) == 0)
287 {
288 nopt++;
289 if (nopt >= argc) IR_error ("Need argument after -rprefix ");
290 (*option_data)->rprefix
291 = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
292 strcpy ((*option_data)->rprefix, argv[nopt]);
293
294 nopt++;
295 continue;
296 }
297
298
299 /*----- -nofine -----*/
300 if (strncmp(argv[nopt], "-nofine", 7) == 0)
301 {
302 (*option_data)->nofine = 1;
303 nopt++;
304 continue;
305 }
306
307
308 /*----- -fine blur dxy dphi -----*/
309 if (strncmp(argv[nopt], "-fine", 5) == 0)
310 {
311 nopt++;
312 if (nopt+2 >= argc) IR_error ("Need 3 arguments after -fine ");
313 (*option_data)->nofine = 0;
314
315 sscanf (argv[nopt], "%f", &fval);
316 if (fval <= 0.0)
317 IR_error ("Illegal argument for blur ( must be > 0 ) ");
318 (*option_data)->blur = fval;
319 nopt++;
320
321 sscanf (argv[nopt], "%f", &fval);
322 if (fval <= 0.0)
323 IR_error ("Illegal argument for dxy ( must be > 0 ) ");
324 (*option_data)->dxy = fval;
325 nopt++;
326
327 sscanf (argv[nopt], "%f", &fval);
328 if (fval <= 0.0)
329 IR_error ("Illegal argument for dphi ( must be > 0 ) ");
330 (*option_data)->dphi = fval;
331 nopt++;
332
333 continue;
334 }
335
336
337 /*----- -dmm -----*/
338 if (strncmp(argv[nopt], "-dmm", 4) == 0)
339 {
340 (*option_data)->dmm = 1;
341 nopt++;
342 continue;
343 }
344
345
346 /*----- -debug -----*/
347 if (strncmp(argv[nopt], "-debug", 6) == 0)
348 {
349 (*option_data)->debug = 1;
350 nopt++;
351 continue;
352 }
353
354
355 /*----- unknown command -----*/
356 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
357 IR_error (message);
358
359
360 }
361
362
363 }
364
365
366 /*---------------------------------------------------------------------------*/
367 /*
368 Routine to read a 3d+time dataset.
369 */
370
read_dataset(char * filename,THD_3dim_dataset ** dset)371 void read_dataset
372 (
373 char * filename, /* file name of 3d+time dataset */
374 THD_3dim_dataset ** dset /* pointer to 3d+time dataset */
375 )
376
377 {
378 char message[80]; /* error message */
379
380
381 /*----- Open the 3d+time dataset -----*/
382 *dset = THD_open_one_dataset (filename);
383 if (*dset == NULL)
384 {
385 sprintf (message,
386 "Unable to open data file: %s", filename);
387 IR_error (message);
388 }
389
390 }
391
392
393 /*---------------------------------------------------------------------------*/
394 /*
395 Initialize the slice sequence arrays.
396 */
397
initialize_slice_sequence(IR_options * option_data,THD_3dim_dataset * dset)398 void initialize_slice_sequence
399 (
400 IR_options * option_data, /* user input options */
401 THD_3dim_dataset * dset /* pointer to 3d+time dataset */
402 )
403
404 {
405 int num_slices; /* number of slices per volume */
406 int ivolume; /* volume index number */
407 int itemp; /* temporary variable */
408 float ttemp; /* temporary variable */
409 float * time_array = NULL; /* array of slice acquisition times */
410 int iz, i, j; /* index numbers */
411 float z; /* slice z location */
412
413
414 if (!dset->taxis) { /* Jan 07 [ZSS] */
415 IR_error ("NULL taxis, should not be here.");
416 }
417
418 /*----- Initialize local variables -----*/
419 num_slices = dset->taxis->nsl;
420
421 ivolume = 0;
422
423 if ( num_slices <= 0 ) /* 06 Oct 2003 [rickr] */
424 num_slices = dset->daxes->nzz;
425
426 /*----- Allocate memory for arrays -----*/
427 t_to_z = (int *) malloc (sizeof(int) * num_slices);
428 z_to_t = (int *) malloc (sizeof(int) * num_slices);
429 time_array = (float *) malloc (sizeof(float) * num_slices);
430
431
432 /*----- Initialize array of slice acquisition times -----*/
433 for (iz = 0; iz < num_slices; iz++)
434 {
435 z = iz * dset->taxis->dz_sl + dset->taxis->zorg_sl;
436 time_array[iz] = THD_timeof (ivolume, z, dset->taxis);
437 t_to_z[iz] = iz;
438 }
439
440
441 /*----- Sort slice z-indices by increasing time -----*/
442 for (i = 0; i < num_slices-1; i++)
443 for (j = i+1; j < num_slices; j++)
444 if (time_array[j] < time_array[i])
445 {
446 itemp = t_to_z[i];
447 t_to_z[i] = t_to_z[j];
448 t_to_z[j] = itemp;
449
450 ttemp = time_array[i];
451 time_array[i] = time_array[j];
452 time_array[j] = ttemp;
453 }
454
455
456 /*----- Sort slice time-indices by increasing z index -----*/
457 for (i = 0; i < num_slices; i++)
458 {
459 j = t_to_z[i];
460 z_to_t[j] = i;
461 }
462
463
464 /*----- Write out the slice ordering arrays -----*/
465 if (option_data->debug)
466 for (i = 0; i < num_slices; i++)
467 printf ("time[%2d] = %12.3f t_to_z[%2d] = %2d z_to_t[%2d] = %2d\n",
468 i, time_array[i], i, t_to_z[i], i, z_to_t[i]);
469
470
471 /*----- Release memory -----*/
472 free (time_array); time_array = NULL;
473 }
474
475
476 /*---------------------------------------------------------------------------*/
477 /*
478 Routine to initialize the array of state vectors.
479 */
480
initialize_state_history(THD_3dim_dataset * dset,vector ** state_history)481 void initialize_state_history
482 (
483 THD_3dim_dataset * dset, /* pointer to input 3d+time dataset */
484 vector ** state_history /* time series of state vectors */
485 )
486
487 {
488 int num_slices; /* number of slices per volume */
489 int ts_length; /* number of volumes */
490 int num_vectors; /* total number of state vectors */
491 int i; /* state vector index */
492
493
494 /*----- Initialize local variables -----*/
495 num_slices = dset->taxis->nsl;
496
497 if ( num_slices <= 0 ) /* 06 Oct 2003 [rickr] */
498 num_slices = dset->daxes->nzz;
499
500 ts_length = DSET_NUM_TIMES(dset);
501 num_vectors = ts_length * num_slices;
502
503
504 /*----- Allocate memory for array of state vectors -----*/
505 *state_history = (vector *) malloc (sizeof(vector) * num_vectors);
506
507
508 /*----- Initialize array of state vectors -----*/
509 for (i = 0; i < num_vectors; i++)
510 {
511 vector_initialize (&((*state_history)[i]));
512 vector_create (STATE_DIM, &((*state_history)[i]));
513 }
514 }
515
516
517 /*---------------------------------------------------------------------------*/
518 /*
519 Routine to initialize the RMS error arrays.
520 */
521
initialize_rms_arrays(THD_3dim_dataset * dset,float ** old_rms_array,float ** new_rms_array)522 void initialize_rms_arrays
523 (
524 THD_3dim_dataset * dset, /* pointer to input 3d+time dataset */
525 float ** old_rms_array, /* volume RMS error for input dataset */
526 float ** new_rms_array /* volume RMS error for registered dataset */
527 )
528
529 {
530 int ts_length; /* number of volumes */
531
532
533 ts_length = DSET_NUM_TIMES(dset);
534
535 /*----- Allocate space for RMS error arrays -----*/
536 *old_rms_array = (float *) malloc (sizeof(float) * ts_length);
537 *new_rms_array = (float *) malloc (sizeof(float) * ts_length);
538
539
540 }
541
542
543 /*---------------------------------------------------------------------------*/
544 /*
545 Routine to perform all program initialization.
546 */
547
initialize_program(int argc,char ** argv,IR_options ** option_data,vector ** state_history,float ** old_rms_array,float ** new_rms_array)548 void initialize_program
549 (
550 int argc, /* number of input arguments */
551 char ** argv, /* array of input arguments */
552 IR_options ** option_data, /* user input options */
553 vector ** state_history, /* time series of state vectors */
554 float ** old_rms_array, /* volume RMS error for input dataset */
555 float ** new_rms_array /* volume RMS error for registered dataset */
556 )
557
558 {
559 THD_3dim_dataset * dset = NULL;
560
561
562 /*----- save command line for history notes -----*/
563 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
564
565
566 /*----- Initialize input options -----*/
567 initialize_options (option_data);
568
569
570 /*----- Get user inputs -----*/
571 get_user_inputs (argc, argv, option_data);
572
573 /*----- Read the input 3d+time dataset to be registered -----*/
574 read_dataset ((*option_data)->input_filename, &dset);
575
576
577 /*----- Initialize the z-slice time order arrays -----*/
578 if (dset->taxis) initialize_slice_sequence (*option_data, dset);
579
580
581 /*----- Initialize the array of state vectors -----*/
582 if ((*option_data)->dprefix != NULL)
583 initialize_state_history (dset, state_history);
584
585
586 /*----- Allocate space for RMS error arrays -----*/
587 if ((*option_data)->rprefix != NULL)
588 initialize_rms_arrays (dset, old_rms_array, new_rms_array);
589
590
591 /*----- Release memory -----*/
592 THD_delete_3dim_dataset (dset, False); dset = NULL;
593
594 }
595
596
597 /*---------------------------------------------------------------------------*/
598 /*
599 Routine to make a copy of a dataset, with data attached.
600 This routine is copied directly from afni_plugin.c
601 */
602
copy_dset(THD_3dim_dataset * dset,char * new_prefix)603 THD_3dim_dataset * copy_dset( THD_3dim_dataset * dset , char * new_prefix )
604 {
605 THD_3dim_dataset * new_dset ;
606 int ival , ityp , nbytes , nvals ;
607 void * new_brick , * old_brick ;
608
609 /*-- sanity check --*/
610
611 if( ! ISVALID_3DIM_DATASET(dset) ) return NULL ;
612
613 /*-- make the empty copy --*/
614
615 new_dset = EDIT_empty_copy( dset ) ;
616
617 /*-- change its name? --*/
618
619 if( new_prefix != NULL )
620 EDIT_dset_items( new_dset ,
621 ADN_prefix , new_prefix ,
622 ADN_label1 , new_prefix ,
623 ADN_none ) ;
624
625 /*-- make brick(s) for this dataset --*/
626
627 THD_load_datablock( dset->dblk ) ; /* make sure old one is in memory */
628
629 nvals = DSET_NVALS(dset) ;
630
631 for( ival=0 ; ival < nvals ; ival++ ){
632 ityp = DSET_BRICK_TYPE(new_dset,ival) ; /* type of data */
633 nbytes = DSET_BRICK_BYTES(new_dset,ival) ; /* how much data */
634 new_brick = malloc( nbytes ) ; /* make room */
635
636 if( new_brick == NULL ){
637 THD_delete_3dim_dataset( new_dset , False ) ;
638 return NULL ;
639 }
640
641 EDIT_substitute_brick( new_dset , ival , ityp , new_brick ) ;
642
643 /*-- copy data from old brick to new brick --*/
644
645 old_brick = DSET_BRICK_ARRAY(dset,ival) ;
646
647 if( old_brick == NULL ){
648 THD_delete_3dim_dataset( new_dset , False ) ;
649 return NULL ;
650 }
651
652 memcpy( new_brick , old_brick , nbytes ) ;
653 }
654
655 return new_dset ;
656 }
657
658
659 /*---------------------------------------------------------------------------*/
660 /*
661 Check whether output file already exists.
662 */
663
check_output_file(THD_3dim_dataset * dset,char * filename)664 void check_output_file
665 (
666 THD_3dim_dataset * dset, /* input afni data set pointer */
667 char * filename /* output file name */
668 )
669
670 {
671 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
672 int ierror; /* number of errors in editing data */
673
674
675 /*-- make an empty copy of the input dataset --*/
676 new_dset = EDIT_empty_copy( dset ) ;
677
678
679 ierror = EDIT_dset_items( new_dset ,
680 ADN_prefix , filename ,
681 ADN_label1 , filename ,
682 ADN_self_name , filename ,
683 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
684 GEN_FUNC_TYPE ,
685 ADN_none ) ;
686
687 if( ierror > 0 ){
688 fprintf(stderr,
689 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
690 exit(1) ;
691 }
692
693 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
694 fprintf(stderr,
695 "*** Output dataset file %s already exists--cannot continue!\a\n",
696 new_dset->dblk->diskptr->header_name ) ;
697 exit(1) ;
698 }
699
700 /*----- deallocate memory -----*/
701 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
702
703 }
704
705
706 /*---------------------------------------------------------------------------*/
707 /*
708 Evaluate accuracy of registration.
709 */
710
eval_registration(IR_options * option_data,THD_3dim_dataset * old_dset,THD_3dim_dataset * new_dset,THD_3dim_dataset * base_dset,int base,float * old_rms_array,float * new_rms_array)711 void eval_registration
712 (
713 IR_options * option_data, /* user input options */
714 THD_3dim_dataset * old_dset, /* input dataset */
715 THD_3dim_dataset * new_dset, /* output datasets */
716 THD_3dim_dataset * base_dset, /* base image dataset */
717 int base, /* base volume index */
718 float * old_rms_array, /* volume RMS error for input dataset */
719 float * new_rms_array /* volume RMS error for registered dataset */
720 )
721
722 {
723 int nx, ny, nz, nxyz;
724 int ix, jy, kz;
725 int ixyz;
726 int datum, base_datum;
727 float old_sse, new_sse;
728 float diff;
729 float old_rmse, new_rmse;
730 int ivolume, num_volumes;
731 byte * bold = NULL, * bnew = NULL, * bbase = NULL;
732 short * sold = NULL, * snew = NULL, * sbase = NULL;
733 float * fold = NULL, * fnew = NULL, * fbase = NULL;
734 float float_base=0.0, float_old=0.0, float_new=0.0;
735
736
737 /*----- Initialize local variables -----*/
738 nx = old_dset->daxes->nxx;
739 ny = old_dset->daxes->nyy;
740 nz = old_dset->daxes->nzz;
741 nxyz = nx * ny * nz;
742 num_volumes = DSET_NUM_TIMES (old_dset);
743 datum = DSET_BRICK_TYPE (new_dset,0) ;
744 base_datum = DSET_BRICK_TYPE (base_dset,0);
745
746
747 /*----- Set base dataset pointer depending on base datum type -----*/
748 switch ( base_datum )
749 {
750 case MRI_byte:
751 bbase = (byte *) DSET_ARRAY(base_dset,base); break;
752
753 case MRI_short:
754 sbase = (short *) DSET_ARRAY(base_dset,base); break;
755
756 case MRI_float:
757 fbase = (float *) DSET_ARRAY(base_dset,base); break;
758 }
759
760
761 for (ivolume = 0; ivolume < num_volumes; ivolume++)
762 {
763 old_sse = 0.0;
764 new_sse = 0.0;
765
766 /*----- Set old and new dataset pointers depending on datum type -----*/
767 switch ( datum )
768 {
769 case MRI_byte:
770 bold = (byte *) DSET_ARRAY(old_dset,ivolume);
771 bnew = (byte *) DSET_ARRAY(new_dset,ivolume); break;
772
773 case MRI_short:
774 sold = (short *) DSET_ARRAY(old_dset,ivolume);
775 snew = (short *) DSET_ARRAY(new_dset,ivolume); break;
776
777 case MRI_float:
778 fold = (float *) DSET_ARRAY(old_dset,ivolume);
779 fnew = (float *) DSET_ARRAY(new_dset,ivolume); break;
780 }
781
782
783 for (kz = 0; kz < nz; kz++)
784 for (jy = 0; jy < ny; jy++)
785 for (ix = 0; ix < nx; ix++)
786 {
787 ixyz = ix + jy*nx + kz*nx*ny;
788
789 /*----- Get base voxel floating point value -----*/
790 switch ( base_datum )
791 {
792 case MRI_byte: float_base = (float) bbase[ixyz]; break;
793
794 case MRI_short: float_base = (float) sbase[ixyz]; break;
795
796 case MRI_float: float_base = fbase[ixyz]; break;
797 }
798
799 /*----- Get old and new voxel floating point value -----*/
800 switch ( datum )
801 {
802 case MRI_byte:
803 float_old = (float) bold[ixyz];
804 float_new = (float) bnew[ixyz]; break;
805
806 case MRI_short:
807 float_old = (float) sold[ixyz];
808 float_new = (float) snew[ixyz]; break;
809
810 case MRI_float:
811 float_old = fold[ixyz];
812 float_new = fnew[ixyz]; break;
813 }
814
815 diff = float_old - float_base;
816 old_sse += diff*diff;
817 diff = float_new - float_base;
818 new_sse += diff*diff;
819 }
820
821 old_rmse = sqrt (old_sse / nxyz);
822 new_rmse = sqrt (new_sse / nxyz);
823
824 if (option_data->debug)
825 printf ("Volume = %d OLD RMSE = %f NEW RMSE = %f \n",
826 ivolume, old_rmse, new_rmse);
827
828 old_rms_array[ivolume] = old_rmse;
829 new_rms_array[ivolume] = new_rmse;
830
831 }
832
833 }
834
835
836 /*---------------------------------------------------------------------------*/
837 /*
838 Main routine for this program.
839 If the return string is not NULL, some error transpired, and
840 the program will display the error message.
841 */
842
843 #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
844
845 #define FREE_WORKSPACE \
846 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
847 FREEUP(bout) ; FREEUP(sout) ; FREEUP(fout) ; \
848 FREEUP(dxar) ; FREEUP(dyar) ; FREEUP(phiar); \
849 } while(0) ;
850
851
IMREG_main(IR_options * opt,vector * state_history,float * old_rms_array,float * new_rms_array)852 char * IMREG_main
853 (
854 IR_options * opt,
855 vector * state_history,
856 float * old_rms_array,
857 float * new_rms_array
858 )
859 {
860 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
861 THD_3dim_dataset * base_dset; /* base image dataset */
862 char * new_prefix ; /* string from user */
863 int base , ntime , datum , nx,ny,nz , ii,kk , npix ;
864 float dx,dy,dz ;
865 int base_datum, failed;
866 MRI_IMARR * ims_in , * ims_out ;
867 MRI_IMAGE * im , * imbase ;
868
869 byte ** bptr = NULL , * bbase = NULL, ** bout = NULL ;
870 short ** sptr = NULL , * sbase = NULL, ** sout = NULL ;
871 float ** fptr = NULL , * fbase = NULL, ** fout = NULL ;
872
873 float * dxar = NULL , * dyar = NULL , * phiar = NULL ;
874
875 int it;
876
877 /*--------------------------------------------------------------------*/
878 /*----- Check batch command inputs to see if they are reasonable -----*/
879
880 old_dset = THD_open_one_dataset(opt->input_filename) ;
881 /* get ptr to dataset */
882 if( old_dset == NULL )
883 return "*************************\n"
884 "Cannot find Input Dataset\n"
885 "*************************" ;
886
887 ntime = DSET_NUM_TIMES(old_dset) ;
888 if( ntime < 2 && !opt->base_filename)
889 return "*****************************\n"
890 "Dataset has only 1 time point\n"
891 " and no -basefile specified.\n"
892 "*****************************" ;
893
894 ii = DSET_NVALS_PER_TIME(old_dset) ;
895 if( ii > 1 )
896 return "************************************\n"
897 "Dataset has > 1 value per time point\n"
898 "************************************" ;
899
900 nx = old_dset->daxes->nxx ; dx = old_dset->daxes->xxdel ;
901 ny = old_dset->daxes->nyy ; dy = old_dset->daxes->yydel ; npix = nx*ny ;
902 nz = old_dset->daxes->nzz ; dz = old_dset->daxes->zzdel ;
903
904 if( nx != ny || fabs(dx) != fabs(dy) ){
905
906 /* No need to quit, works fine. ZSS 07
907 * Only if nx >= ny (so fix might be easy). 12 Jan 2010 [rickr] */
908 if (opt->debug)
909 fprintf(stderr,"\nNotice 2dImreg: nx=%d ny=%d nz=%d dx=%f dy=%f dz=%f\n",
910 nx,ny,nz,dx,dy,dz ) ;
911
912 return "***********************************\n"
913 "Dataset does not have square slices\n"
914 "***********************************" ;
915 }
916
917 new_prefix = opt->new_prefix; /* get string item (the output prefix) */
918 if (new_prefix == NULL) /* check if it is OK */
919 return "************************\n"
920 "Output Prefix is illegal\n"
921 "************************" ;
922
923 /*----- Check whether output file already exists -----*/
924 check_output_file (old_dset, new_prefix);
925
926
927 /*--------- go to "base" input option ---------*/
928
929 if (opt->base_filename == NULL)
930 base_dset = old_dset;
931 else
932 {
933 base_dset = THD_open_one_dataset(opt->base_filename) ;
934 /* get ptr to base dataset */
935 if( base_dset == NULL )
936 return "************************\n"
937 "Cannot find Base Dataset\n"
938 "************************" ;
939
940 if ( (nx != base_dset->daxes->nxx) || (dx != base_dset->daxes->xxdel)
941 || (ny != base_dset->daxes->nyy) || (dy != base_dset->daxes->yydel)
942 || (nz != base_dset->daxes->nzz) || (dz != base_dset->daxes->zzdel) )
943 {
944 if (opt->debug)
945 {
946 fprintf(stderr,"\nIMREG: Input Dataset:\n");
947 fprintf(stderr,"nx=%d ny=%d nz=%d dx=%f dy=%f dz=%f\n",
948 nx,ny,nz,dx,dy,dz ) ;
949
950 fprintf(stderr,"\nIMREG: Base Dataset:\n");
951 fprintf(stderr,"nx=%d ny=%d nz=%d dx=%f dy=%f dz=%f\n",
952 base_dset->daxes->nxx, base_dset->daxes->nyy,
953 base_dset->daxes->nzz, base_dset->daxes->xxdel,
954 base_dset->daxes->yydel, base_dset->daxes->zzdel) ;
955 }
956 return "*************************************************\n"
957 "Base Dataset is not compatible with Input Dataset\n"
958 "*************************************************" ;
959 }
960 }
961
962 base_datum = DSET_BRICK_TYPE(base_dset,0);
963
964 base = opt->base_vol_index;
965 if (base < 0) { /* default */
966 if (-base >= DSET_NUM_TIMES(base_dset)) {
967 base = DSET_NUM_TIMES(base_dset)-1;
968 } else {
969 base = -base;
970 }
971 }
972
973 if( base >= DSET_NUM_TIMES(base_dset) )
974 return "******************************\n"
975 "Base image number is too large\n"
976 "******************************" ;
977
978
979 /*--------- see if the "fine" input option is present --------*/
980 if (opt->nofine)
981 mri_align_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
982 else{
983 float fsig , fdxy , fdph ;
984 fsig = opt->blur * 0.42466090;
985 fdxy = opt->dxy;
986 fdph = opt->dphi;
987 mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
988
989 if (opt->debug)
990 fprintf(stderr,"Set fine params = %f %f %f\n",fsig,fdxy,fdph) ;
991 }
992
993 /*------------- ready to compute new dataset -----------*/
994
995 if (opt->debug)
996 fprintf(stderr,"IMREG: loading dataset\n") ;
997
998
999 DSET_load( old_dset ) ;
1000
1001 if (opt->base_filename != NULL)
1002 DSET_load (base_dset);
1003
1004 /*** 1) Copy the dataset in toto ***/
1005
1006 if (opt->debug)
1007 fprintf(stderr,"IMREG: Copying dataset\n") ;
1008
1009
1010 new_dset = copy_dset( old_dset , new_prefix ) ;
1011 if( new_dset == NULL )
1012 return "****************************\n"
1013 "Failed to copy input dataset\n"
1014 "****************************" ;
1015
1016 /*** 2) Make an array of empty images ***/
1017
1018 if (opt->debug)
1019 fprintf(stderr,"IMREG: making empty images\n") ;
1020
1021
1022 datum = DSET_BRICK_TYPE(new_dset,0) ;
1023
1024 INIT_IMARR(ims_in) ;
1025 for( ii=0 ; ii < ntime ; ii++ ){
1026 im = mri_new_vol_empty( nx , ny , 1 , datum ) ;
1027 ADDTO_IMARR(ims_in,im) ;
1028 }
1029
1030 imbase = mri_new_vol_empty( nx , ny , 1 , base_datum ) ;
1031
1032 dxar = (float *) malloc( sizeof(float) * ntime ) ;
1033 dyar = (float *) malloc( sizeof(float) * ntime ) ;
1034 phiar = (float *) malloc( sizeof(float) * ntime ) ;
1035
1036 /*** 3) Get pointers to sub-bricks in old, base, and new datasets ***/
1037
1038 if (opt->debug)
1039 fprintf(stderr,"IMREG: getting input brick pointers\n") ;
1040
1041
1042 switch( datum ){ /* pointer type depends on input datum type */
1043 case MRI_byte:
1044 bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
1045 bout = (byte **) malloc( sizeof(byte *) * ntime ) ;
1046 for( ii=0 ; ii < ntime ; ii++ ){
1047 bptr[ii] = (byte *) DSET_ARRAY(old_dset,ii) ;
1048 bout[ii] = (byte *) DSET_ARRAY(new_dset,ii) ;
1049 }
1050 break ;
1051
1052 case MRI_short:
1053 sptr = (short **) malloc( sizeof(short *) * ntime ) ;
1054 sout = (short **) malloc( sizeof(short *) * ntime ) ;
1055 for( ii=0 ; ii < ntime ; ii++ ){
1056 sptr[ii] = (short *) DSET_ARRAY(old_dset,ii) ;
1057 sout[ii] = (short *) DSET_ARRAY(new_dset,ii) ;
1058 }
1059
1060 if (opt->debug)
1061 fprintf(stderr,"IMREG: sptr[0] = %p sout[0] = %p\n",
1062 sptr[0],sout[0]) ;
1063
1064 break ;
1065
1066 case MRI_float:
1067 fptr = (float **) malloc( sizeof(float *) * ntime ) ;
1068 fout = (float **) malloc( sizeof(float *) * ntime ) ;
1069 for( ii=0 ; ii < ntime ; ii++ ){
1070 fptr[ii] = (float *) DSET_ARRAY(old_dset,ii) ;
1071 fout[ii] = (float *) DSET_ARRAY(new_dset,ii) ;
1072 }
1073 break ;
1074 }
1075
1076 switch( base_datum ){ /* pointer type depends on base datum type */
1077 case MRI_byte:
1078 bbase = (byte *) DSET_ARRAY(base_dset,base) ;
1079 break ;
1080
1081 case MRI_short:
1082 sbase = (short *) DSET_ARRAY(base_dset,base) ;
1083 if (opt->debug)
1084 fprintf(stderr,"IMREG: sbase[%d] = %p \n", base, sbase) ;
1085 break ;
1086
1087 case MRI_float:
1088 fbase = (float *) DSET_ARRAY(base_dset,base) ;
1089 break ;
1090 }
1091
1092 /*** 4) Loop over slices ***/
1093
1094 for( kk=0 ; kk < nz ; kk++ ){
1095
1096 /*** 4a) Setup ims_in images to point to input slices ***/
1097
1098 if (opt->debug)
1099 fprintf(stderr,"IMREG: slice %d -- setup input images\n",kk) ;
1100
1101
1102 for( ii=0 ; ii < ntime ; ii++ ){
1103 im = IMARR_SUBIMAGE(ims_in,ii) ;
1104 switch( datum ){
1105 case MRI_byte:
1106 mri_fix_data_pointer( bptr[ii] + kk*npix, im ) ; break ;
1107 case MRI_short:
1108 mri_fix_data_pointer( sptr[ii] + kk*npix, im ) ; break ;
1109 case MRI_float:
1110 mri_fix_data_pointer( fptr[ii] + kk*npix, im ) ; break ;
1111 }
1112 }
1113
1114 /*** 4b) Setup imbase to point to base image ***/
1115
1116 if (opt->debug)
1117 fprintf(stderr,"IMREG: slice %d -- setup base image\n",kk) ;
1118
1119
1120 switch( base_datum ){
1121 case MRI_byte:
1122 mri_fix_data_pointer( bbase + kk*npix, imbase ) ; break ;
1123 case MRI_short:
1124 mri_fix_data_pointer( sbase + kk*npix, imbase ) ; break ;
1125 case MRI_float:
1126 mri_fix_data_pointer( fbase + kk*npix, imbase ) ; break ;
1127 }
1128
1129 /*** 4c) Register this slice at all times ***/
1130
1131 if (opt->debug)
1132 fprintf(stderr,"IMREG: slice %d -- register\n",kk) ;
1133
1134
1135 ims_out = mri_align_dfspace( imbase , NULL , ims_in ,
1136 ALIGN_REGISTER_CODE , dxar,dyar,phiar ) ;
1137
1138 /* hack of a fix for failure case 08 Apr 2008 [rickr/dglen] */
1139 failed = 0;
1140 if( ims_out == NULL ) {
1141 fprintf(stderr,"IMREG failure (zero slice?), copying input\n");
1142 failed = 1;
1143 ims_out = ims_in;
1144 }
1145
1146
1147 /*----- Store the registration parameters -----*/
1148 if (opt->dprefix != NULL)
1149 for (ii = 0; ii < ntime; ii++)
1150 {
1151 it = ii*nz + z_to_t[kk];
1152 if (opt->dmm)
1153 {
1154 state_history[it].elts[1] = dxar[ii] * dx;
1155 state_history[it].elts[2] = dyar[ii] * dy;
1156 }
1157 else
1158 {
1159 state_history[it].elts[1] = dxar[ii];
1160 state_history[it].elts[2] = dyar[ii];
1161 }
1162
1163 state_history[it].elts[3] = (180.0/PI)*phiar[ii];
1164 }
1165
1166
1167 /*** 4d) Put the output back in on top of the input;
1168 note that the output is always in MRI_float format ***/
1169
1170 if (opt->debug)
1171 fprintf(stderr,"IMREG: slice %d -- put output back into dataset\n",kk);
1172
1173
1174 for( ii=0 ; ii < ntime ; ii++ ){
1175 switch( datum ){
1176 case MRI_byte:
1177 im = mri_to_mri( MRI_byte , IMARR_SUBIMAGE(ims_out,ii) ) ;
1178 memcpy( bout[ii] + kk*npix , MRI_BYTE_PTR(im) ,
1179 sizeof(byte)*npix ) ;
1180 mri_free(im) ;
1181 break ;
1182
1183 case MRI_short:
1184
1185 if (opt->debug)
1186 if( ii==0 )
1187 fprintf(stderr,"IMREG: conversion to short at ii=%d\n",ii) ;
1188
1189 im = mri_to_mri( MRI_short , IMARR_SUBIMAGE(ims_out,ii) ) ;
1190
1191 if (opt->debug)
1192 if( ii==0 )
1193 fprintf(stderr,"IMREG: copying to %p from %p\n",
1194 sout[ii] + kk*npix,MRI_SHORT_PTR(im)) ;
1195
1196
1197 memcpy( sout[ii] + kk*npix , MRI_SHORT_PTR(im) ,
1198 sizeof(short)*npix ) ;
1199
1200 if (opt->debug)
1201 if( ii==0 )
1202 fprintf(stderr,"IMREG: freeing\n") ;
1203
1204
1205 mri_free(im) ;
1206 break ;
1207
1208 case MRI_float:
1209 im = IMARR_SUBIMAGE(ims_out,ii) ;
1210 memcpy( fout[ii] + kk*npix , MRI_FLOAT_PTR(im) ,
1211 sizeof(float)*npix ) ;
1212 break ;
1213 }
1214 }
1215
1216 /*** 4e) Destroy the output images ***/
1217
1218 if (opt->debug)
1219 fprintf(stderr,"IMREG: destroying aligned output\n") ;
1220
1221
1222 if( !failed ) DESTROY_IMARR( ims_out ) ; /* 08 Apr 2008 */
1223 }
1224
1225 /*** 5) Destroy the empty images and other workspaces ***/
1226
1227 if (opt->debug)
1228 fprintf(stderr,"IMREG: destroy workspaces\n") ;
1229
1230
1231 mri_clear_data_pointer(imbase) ; mri_free(imbase) ;
1232 for( ii=0 ; ii < ntime ; ii++ ){
1233 im = IMARR_SUBIMAGE(ims_in,ii) ;
1234 mri_clear_data_pointer(im) ;
1235 }
1236 DESTROY_IMARR(ims_in) ;
1237 FREE_WORKSPACE ;
1238
1239 /*------------- write out the new dataset ------------*/
1240
1241 if (opt->debug)
1242 fprintf(stderr,"IMREG: write new dataset to disk\n") ;
1243
1244
1245 /*----- Record history of dataset -----*/
1246 tross_Copy_History( old_dset , new_dset ) ;
1247 if( commandline != NULL )
1248 tross_Append_History( new_dset , commandline ) ;
1249
1250
1251 THD_load_statistics( new_dset ) ;
1252 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
1253
1254
1255 /*----- evaluate results -----*/
1256 if (opt->rprefix != NULL)
1257 eval_registration (opt, old_dset, new_dset, base_dset, base,
1258 old_rms_array, new_rms_array);
1259
1260
1261 /*----- deallocate memory -----*/
1262 THD_delete_3dim_dataset( old_dset , False ) ; old_dset = NULL ;
1263 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
1264 if (opt->base_filename != NULL)
1265 THD_delete_3dim_dataset( base_dset , False ) ; base_dset = NULL ;
1266
1267
1268
1269 return NULL ; /* null string returned means all was OK */
1270 }
1271
1272
1273 /*---------------------------------------------------------------------------*/
1274 /*
1275 Routine to write RMS error arrays.
1276 */
1277
1278
output_rms_arrays(IR_options * option_data,THD_3dim_dataset * dset,float * old_rms_array,float * new_rms_array)1279 void output_rms_arrays
1280 (
1281 IR_options * option_data, /* user input options */
1282 THD_3dim_dataset * dset, /* pointer to input 3d+time dataset */
1283 float * old_rms_array, /* volume RMS error for input dataset */
1284 float * new_rms_array /* volume RMS error for registered dataset */
1285 )
1286
1287 {
1288 int it; /* time index */
1289 int ts_length; /* length of time series data */
1290 char filename[MAX_NAME_LENGTH]; /* name of output file */
1291 FILE * old_rms_file , * new_rms_file; /* file for volume RMS error */
1292
1293
1294 /*----- Initialize local variables -----*/
1295 ts_length = DSET_NUM_TIMES(dset);
1296
1297
1298 /*----- Open output files -----*/
1299 strcpy (filename, option_data->rprefix);
1300 strcat (filename, ".oldrms");
1301 old_rms_file = fopen (filename, "w");
1302 strcpy (filename, option_data->rprefix);
1303 strcat (filename, ".newrms");
1304 new_rms_file = fopen (filename, "w");
1305
1306
1307 /*----- Write volume RMS error data -----*/
1308 for (it = 0; it < ts_length; it++)
1309 {
1310 fprintf (old_rms_file, "%d %f\n" , it, old_rms_array[it]);
1311 fprintf (new_rms_file, "%d %f\n" , it, new_rms_array[it]);
1312 }
1313
1314
1315 /*----- Close output files -----*/
1316 fclose (old_rms_file);
1317 fclose (new_rms_file);
1318
1319 }
1320
1321
1322 /*---------------------------------------------------------------------------*/
1323 /*
1324 Get the time corresponding to this particular slice.
1325 */
1326
get_time(int ivolume,int iz,THD_3dim_dataset * dset)1327 float get_time (int ivolume, int iz, THD_3dim_dataset * dset)
1328
1329 {
1330 float time;
1331 float z;
1332
1333
1334 z = iz * dset->taxis->dz_sl + dset->taxis->zorg_sl;
1335 time = THD_timeof (ivolume, z, dset->taxis);
1336
1337 if (dset->taxis->units_type == UNITS_MSEC_TYPE)
1338 time /= 1000.0;
1339
1340 return (time);
1341 }
1342
1343
1344 /*---------------------------------------------------------------------------*/
1345 /*
1346 Routine to write the registration parameter time series.
1347 */
1348
1349
output_state_history(IR_options * option_data,THD_3dim_dataset * dset,vector * state_history)1350 void output_state_history
1351 (
1352 IR_options * option_data, /* user input options */
1353 THD_3dim_dataset * dset,
1354 vector * state_history
1355 )
1356
1357 {
1358 int iv; /* vector index */
1359 int ts_length; /* length of time series */
1360 int num_slices; /* number of slices in each volume */
1361 int num_vectors; /* total number of state vectors */
1362 int ivolume; /* volume index */
1363 int iz; /* z slice index */
1364 float t; /* time for current vector */
1365 char filename[MAX_NAME_LENGTH]; /* string for output file name */
1366
1367 FILE * dx_file, * dy_file, * psi_file; /* output files */
1368
1369
1370 /*----- Initialize local variables -----*/
1371 num_slices = dset->taxis->nsl;
1372 ts_length = DSET_NUM_TIMES(dset);
1373
1374 if ( num_slices <= 0 ) /* 06 Oct 2003 [rickr] */
1375 num_slices = dset->daxes->nzz;
1376
1377
1378 /*----- Calculate total number of state vectors -----*/
1379 num_vectors = ts_length * num_slices;
1380
1381
1382 /*----- Open output files -----*/
1383 strcpy (filename, option_data->dprefix);
1384 strcat (filename, ".dx");
1385 dx_file = fopen (filename, "w");
1386
1387 strcpy (filename, option_data->dprefix);
1388 strcat (filename, ".dy");
1389 dy_file = fopen (filename, "w");
1390
1391 strcpy (filename, option_data->dprefix);
1392 strcat (filename, ".psi");
1393 psi_file = fopen (filename, "w");
1394
1395
1396 /*----- Write the registration parameters -----*/
1397 if (dset->taxis) {
1398 for (iv = 0; iv < num_vectors; iv++)
1399 {
1400 ivolume = iv / num_slices;
1401 iz = t_to_z[iv % num_slices];
1402 t = get_time (ivolume, iz, dset);
1403 fprintf (dx_file, "%f %f\n", t, state_history[iv].elts[1]);
1404 fprintf (dy_file, "%f %f\n", t, state_history[iv].elts[2]);
1405 fprintf (psi_file, "%f %f\n", t, state_history[iv].elts[3]);
1406 }
1407 } else {
1408 for (iv = 0; iv < num_vectors; iv++)
1409 {
1410 fprintf (dx_file, "%d %f\n", iv, state_history[iv].elts[1]);
1411 fprintf (dy_file, "%d %f\n", iv, state_history[iv].elts[2]);
1412 fprintf (psi_file, "%d %f\n", iv, state_history[iv].elts[3]);
1413 }
1414 }
1415
1416
1417 /*----- Close output files -----*/
1418 fclose (dx_file);
1419 fclose (dy_file);
1420 fclose (psi_file);
1421 }
1422
1423
1424 /*---------------------------------------------------------------------------*/
1425 /*
1426 Output results.
1427 */
1428
output_results(IR_options * option_data,vector * state_history,float * old_rms_array,float * new_rms_array)1429 void output_results
1430 (
1431 IR_options * option_data, /* user input options */
1432 vector * state_history, /* time series of state vectors */
1433 float * old_rms_array, /* volume RMS error for input dataset */
1434 float * new_rms_array /* volume RMS error for registered dataset */
1435 )
1436
1437 {
1438 THD_3dim_dataset * dset;
1439
1440 read_dataset (option_data->input_filename, &dset);
1441
1442 /*----- Write the time series of state parameters -----*/
1443 if (option_data->dprefix != NULL)
1444 output_state_history (option_data, dset, state_history);
1445
1446
1447 /*----- Write user specified auxiliary time series data -----*/
1448 if (option_data->rprefix != NULL)
1449 output_rms_arrays (option_data, dset, old_rms_array, new_rms_array);
1450
1451
1452 THD_delete_3dim_dataset (dset, False); dset = NULL;
1453
1454 }
1455
1456
1457 /*---------------------------------------------------------------------------*/
1458 /*
1459 Terminate the program.
1460 */
1461
terminate_program(IR_options ** option_data,vector ** state_history,float ** old_rms_array,float ** new_rms_array)1462 void terminate_program
1463 (
1464 IR_options ** option_data,
1465 vector ** state_history, /* time series of state vectors */
1466 float ** old_rms_array, /* volume RMS error for input dataset */
1467 float ** new_rms_array /* volume RMS error for registered dataset */
1468 )
1469
1470 {
1471 THD_3dim_dataset * dset = NULL; /* pointer to input 3d+time dataset */
1472 int num_slices; /* number of slices in each volume */
1473 int ts_length; /* length of time series */
1474 int num_vectors; /* total number of state vectors */
1475 int i; /* index */
1476
1477
1478 /*----- Initialize local variables -----*/
1479 read_dataset ((*option_data)->input_filename, &dset);
1480 if (dset->taxis) num_slices = dset->taxis->nsl;
1481 else num_slices = -1;
1482
1483 if ( num_slices <= 0 ) /* 06 Oct 2003 [rickr] */
1484 num_slices = dset->daxes->nzz;
1485
1486 ts_length = DSET_NUM_TIMES(dset);
1487 num_vectors = ts_length * num_slices;
1488 THD_delete_3dim_dataset (dset, False); dset = NULL;
1489
1490
1491 /*----- Release memory -----*/
1492 free (*option_data); *option_data = NULL;
1493 if (t_to_z) free (t_to_z); t_to_z = NULL;
1494 if (z_to_t) free (z_to_t); z_to_t = NULL;
1495
1496
1497 if (*old_rms_array != NULL)
1498 { free (*old_rms_array); *old_rms_array = NULL; }
1499 if (*new_rms_array != NULL)
1500 { free (*new_rms_array); *new_rms_array = NULL; }
1501
1502
1503 /*----- Deallocate memory for array of state vectors -----*/
1504 if (*state_history != NULL)
1505 {
1506 for (i = 0; i < num_vectors; i++)
1507 vector_destroy (&((*state_history)[i]));
1508 free (*state_history); *state_history = NULL;
1509 }
1510
1511 }
1512
1513
1514 /*---------------------------------------------------------------------------*/
1515
main(int argc,char ** argv)1516 int main
1517 (
1518 int argc, /* number of input arguments */
1519 char ** argv /* array of input arguments */
1520 )
1521
1522 {
1523 IR_options * option_data = NULL; /* user input options */
1524 char * chptr; /* error message from processing */
1525 vector * state_history = NULL; /* time series of state vectors */
1526 float * old_rms_array = NULL; /* original data volume RMS error */
1527 float * new_rms_array = NULL; /* registered data volume RMS error */
1528
1529
1530
1531 #if 0
1532 /*----- Identify software -----*/
1533 printf ("\n\n");
1534 printf ("Program: %s \n", PROGRAM_NAME);
1535 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
1536 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
1537 printf ("\n");
1538 #endif
1539
1540 machdep() ; mainENTRY("2dImReg main") ; PRINT_VERSION("2dImReg") ;
1541
1542 /*----- Program initialization -----*/
1543 initialize_program (argc, argv, &option_data, &state_history,
1544 &old_rms_array, &new_rms_array);
1545
1546
1547 /*----- Register all slices in the dataset -----*/
1548 chptr = IMREG_main (option_data, state_history,
1549 old_rms_array, new_rms_array);
1550
1551
1552 /*----- Check for processing errors -----*/
1553 if (chptr != NULL)
1554 {
1555 printf ("%s \n\n", chptr);
1556 exit(1);
1557 }
1558
1559
1560 /*----- Output the results -----*/
1561 output_results (option_data, state_history,
1562 old_rms_array, new_rms_array);
1563
1564
1565 /*----- Terminate the program -----*/
1566 terminate_program (&option_data, &state_history,
1567 &old_rms_array, &new_rms_array);
1568
1569 return 0;
1570 }
1571
1572