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 #include <string.h>
9
10 /******* global data *******/
11
12 /** define results of scanning the command line **/
13
14 static char RG_out_prefix[256] = "reg." ;
15 static char RG_out_suffix[256] = "\0" ;
16 static int RG_out_start = 1 ;
17 static int RG_out_step = 1 ;
18 static int RG_out_mode = -666 ;
19 static int RG_keepsize = 0 ;
20 static char RG_dout_prefix[256] = "\0" ;
21
22 static int RG_meth = ALIGN_DFSPACE_TYPE ;
23 static int RG_methcode = 0 ;
24 static int RG_verbose = 1 ;
25 static int RG_skip_output = 0 ;
26 static int RG_debug = 0 ;
27
28 static int RG_almode_coarse = MRI_BICUBIC ;
29 static int RG_almode_fine = MRI_BICUBIC ;
30 static int RG_almode_reg = MRI_BICUBIC ;
31
32 static int RG_use_cmass = 0 ; /* 12 Nov 2001 */
33
34 static MRI_IMAGE * RG_imwt = NULL ;
35
36 static int Iarg = 1 ;
37 static int Argc ;
38 static char ** Argv ;
39
40 static MRI_IMAGE * RG_baseimage = NULL ;
41 static MRI_IMARR * RG_imseq = NULL ;
42
43 /******* prototypes *******/
44
45 void REG_syntax(void) ;
46 void REG_command_line(void) ;
47
48 /******* the program! *****/
49
50 /*-- 07 Apr 1998: for testing the new mri_2dalign_* routines --*/
51
52 #undef USE_2DALIGN
53 #ifdef USE_2DALIGN
54 # undef ALLOW_DFTIME
55 #endif
56
main(int argc,char * argv[])57 int main( int argc , char *argv[] )
58 {
59 MRI_IMAGE * qim , * tim=NULL ;
60 MRI_IMARR * regar = NULL ;
61 float * dx , * dy , * phi ;
62 int kim , imcount ;
63 char fname[666] ;
64 float dxtop , dytop , phitop ;
65 FILE * dxfil , * dyfil , * phifil ;
66
67 /** handle command line options **/
68
69 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){ REG_syntax() ; exit(0); }
70
71 machdep() ;
72
73 Argc = argc ;
74 Argv = argv ;
75 Iarg = 1 ;
76
77 REG_command_line() ;
78
79 imcount = RG_imseq->num ;
80 dx = (float *) malloc( sizeof(float) * imcount ) ;
81 dy = (float *) malloc( sizeof(float) * imcount ) ;
82 phi = (float *) malloc( sizeof(float) * imcount ) ;
83 if( dx == NULL || dy == NULL || phi == NULL ){
84 fprintf(stderr,"** malloc failure for mri_rotate parameters!\a\n") ; exit(1) ;
85 }
86
87 /* 12 Nov 2001: shift input images to align centers of mass */
88
89 if( RG_use_cmass ){
90 float xbase,ybase , xcm,ycm ;
91 int ii,jj,kk , di,dj , nx,ny , sj,si ;
92 MRI_IMAGE *fim , *nim ;
93 float *far , *nar ;
94
95 if( RG_verbose ) printf("-- doing -cmass adjustments\n") ;
96
97 mri_get_cmass_2D( RG_baseimage , &xbase , &ybase ) ;
98
99 for( kk=0 ; kk < imcount ; kk++ ){
100 fim = mri_to_float( IMARR_SUBIM(RG_imseq,kk) ) ; /* copy image */
101 mri_get_cmass_2D( fim , &xcm , &ycm ) ;
102 di = (int) rint(xbase-xcm) ; /* integer shifts */
103 dj = (int) rint(ybase-ycm) ;
104 if( di != 0 || dj != 0 ){ /* shift input image */
105
106 #if 0
107 fprintf(stderr,"Image %d: xbase=%g ybase=%g xcm=%g ycm=%g di=%d dj=%d\n",
108 kk,xbase,ybase,xcm,ycm,di,dj ) ;
109 #endif
110
111 nim = mri_new_conforming( fim , MRI_float ); /* zero filled */
112 far = MRI_FLOAT_PTR(fim) ;
113 nar = MRI_FLOAT_PTR(nim) ;
114 nx = fim->nx ; ny = fim->ny ;
115 for( jj=0 ; jj < ny ; jj++ ){ /* copy fim into nim */
116 sj = jj + dj ;
117 if( sj < 0 || sj >= ny ) continue ;
118 for( ii=0 ; ii < nx ; ii++ ){
119 si = ii + di ;
120 if( si < 0 || si >= nx ) continue ;
121 nar[si+nx*sj] = far[ii+nx*jj] ;
122 }
123 }
124 mri_free( IMARR_SUBIM(RG_imseq,kk) ) ; /* replace image */
125 IMARR_SUBIM(RG_imseq,kk) = nim ; /* in the array */
126 }
127 mri_free(fim) ;
128 }
129 } /* end of RG_use_cmass */
130
131 /* do the actual iterative alignments */
132
133 if( RG_verbose ) printf("-- beginning alignment\n") ;
134 switch( RG_meth ){
135
136 default:
137 case ALIGN_DFSPACE_TYPE:
138 #ifdef USE_2DALIGN
139 regar = mri_2dalign_many( RG_baseimage,RG_imwt, RG_imseq, dx,dy,phi ) ;
140 #else
141 if( ! RG_skip_output ) RG_methcode |= ALIGN_REGISTER_CODE ;
142 regar = mri_align_dfspace( RG_baseimage,RG_imwt, RG_imseq, RG_methcode, dx,dy,phi ) ;
143 #endif
144 break ;
145
146 #ifdef ALLOW_DFTIME
147 case ALIGN_DFTIME_TYPE:
148 if( ! RG_skip_output ) RG_methcode |= ALIGN_REGISTER_CODE ;
149 regar = mri_align_dftime( RG_baseimage,RG_imwt, RG_imseq, RG_methcode, dx,dy,phi ) ;
150 break ;
151 #endif
152 }
153
154 dxtop = dytop = phitop = 0.0 ;
155
156 if( strlen(RG_dout_prefix) > 0 ){
157 sprintf( fname , "%s%s" , RG_dout_prefix , "dx" ) ; dxfil = fopen( fname , "w" ) ;
158 sprintf( fname , "%s%s" , RG_dout_prefix , "dy" ) ; dyfil = fopen( fname , "w" ) ;
159 sprintf( fname , "%s%s" , RG_dout_prefix , "phi" ) ; phifil = fopen( fname , "w" ) ;
160 } else {
161 dxfil = dyfil = phifil = NULL ;
162 }
163
164 for( kim=0 ; kim < imcount ; kim++ ){
165
166 if( fabs(dx[kim]) > fabs(dxtop) ) dxtop = dx[kim] ;
167 if( fabs(dy[kim]) > fabs(dytop) ) dytop = dy[kim] ;
168 if( fabs(phi[kim]) > fabs(phitop) ) phitop = phi[kim] ;
169
170 if( ! RG_skip_output ){
171 sprintf( fname , "%s%04d%s" , RG_out_prefix ,
172 RG_out_start + kim*RG_out_step , RG_out_suffix ) ;
173
174 if( regar == NULL ){
175 tim = mri_rota_variable( RG_almode_reg ,
176 RG_imseq->imarr[kim] ,
177 dx[kim],dy[kim],phi[kim] ) ;
178 } else {
179 tim = regar->imarr[kim] ;
180 }
181 } else {
182 sprintf( fname , "%4d" , kim+1 ) ;
183 }
184
185 if( RG_verbose )
186 printf("-- image %s: dx = %6.3f dy = %6.3f phi = %6.3f\n" ,
187 fname,dx[kim],dy[kim],(180.0/PI)*phi[kim] ) ;
188
189 if( dxfil != NULL ) fprintf( dxfil , "%6.3f\n" , dx[kim] ) ;
190 if( dyfil != NULL ) fprintf( dyfil , "%6.3f\n" , dy[kim] ) ;
191 if( phifil != NULL ) fprintf( phifil , "%6.3f\n" , (180.0/PI)*phi[kim] ) ;
192
193 if (RG_keepsize) {
194 int kkk;
195 MRI_IMARR * trimarr = mri_uncat2D( RG_baseimage->nx , RG_baseimage->ny , tim );
196 if( RG_verbose )
197 printf("-- image cropped back to original size of %dx%d\n", RG_baseimage->nx , RG_baseimage->ny);
198 mri_free(tim); tim = IMARR_SUBIMAGE(trimarr,0) ;
199 for (kkk=1; kkk<trimarr->num; ++kkk) mri_free(IMARR_SUBIMAGE(trimarr,kkk));
200 }
201
202 if( ! RG_skip_output ){
203 switch( RG_out_mode ){
204
205 default:
206 mri_write( fname , tim ) ;
207 break ;
208
209 case MRI_short:
210 qim = mri_to_short( 1.0 , tim ) ;
211 mri_write( fname , qim ) ; mri_free( qim ) ;
212 break ;
213
214 case MRI_byte:
215 qim = mri_to_byte( tim ) ;
216 mri_write( fname , qim ) ; mri_free( qim ) ;
217 break ;
218 }
219 mri_free( tim ) ; mri_free( RG_imseq->imarr[kim] ) ;
220 }
221 }
222
223 if( RG_verbose )
224 printf("-- MAX: %*s dx = %6.3f dy = %6.3f phi = %6.3f\n" ,
225 (int)strlen(fname) , " " , dxtop,dytop,phitop*(180.0/PI) ) ;
226
227 exit(0) ;
228 }
229
230 /*---------------------------------------------------------------------*/
231
REG_syntax(void)232 void REG_syntax(void)
233 {
234 printf(
235 "Usage: imreg [options] base_image image_sequence ...\n"
236 " * Registers each 2D image in 'image_sequence' to 'base_image'.\n"
237 " * If 'base_image' = '+AVER', will compute the base image as\n"
238 " the average of the images in 'image_sequence'.\n"
239 " * If 'base_image' = '+count', will use the count-th image in the\n"
240 " sequence as the base image. Here, count is 1,2,3, ....\n"
241 "\n"
242 "OUTPUT OPTIONS:\n"
243 " -nowrite Don't write outputs, just print progress reports.\n"
244 " -prefix pname The output files will be named in the format\n"
245 " -suffix sname 'pname.index.sname' where 'pname' and 'sname'\n"
246 " -start si are strings given by the first 2 options.\n"
247 " -step ss 'index' is a number, given by 'si+(i-1)*ss'\n"
248 " for the i-th output file, for i=1,2,...\n"
249 " *** Default pname = 'reg.'\n"
250 " *** Default sname = nothing at all\n"
251 " *** Default si = 1\n"
252 " *** Default ss = 1\n"
253 "\n"
254 " -flim Write output in mrilib floating point format\n"
255 " (which can be converted to shorts using program ftosh).\n"
256 " *** Default is to write images in format of first\n"
257 " input file in the image_sequence.\n"
258 " -keepsize Preserve the original image size on output.\n"
259 " Default is without this option, and results\n"
260 " in images that are padded to be square.\n"
261 "\n"
262 " -quiet Don't write progress report messages.\n"
263 " -debug Write lots of debugging output!\n"
264 "\n"
265 " -dprefix dname Write files 'dname'.dx, 'dname'.dy, 'dname'.phi\n"
266 " for use in time series analysis.\n"
267 "\n"
268 "ALIGNMENT ALGORITHMS:\n"
269 " -bilinear Uses bilinear interpolation during the iterative\n"
270 " adjustment procedure, rather than the default\n"
271 " bicubic interpolation. NOT RECOMMENDED!\n"
272 " -modes c f r Uses interpolation modes 'c', 'f', and 'r' during\n"
273 " the coarse, fine, and registration phases of the\n"
274 " algorithm, respectively. The modes can be selected\n"
275 " from 'bilinear', 'bicubic', and 'Fourier'. The\n"
276 " default is '-modes bicubic bicubic bicubic'.\n"
277 " -mlcF Equivalent to '-modes bilinear bicubic Fourier'.\n"
278 "\n"
279 " -wtim filename Uses the image in 'filename' as a weighting factor\n"
280 " for each voxel (the larger the value the more\n"
281 " importance is given to that voxel).\n"
282 "\n"
283 " -dfspace[:0] Uses the 'iterated diffential spatial' method to\n"
284 " align the images. The optional :0 indicates to\n"
285 " skip the iteration of the method, and to use the\n"
286 " simpler linear differential spatial alignment method.\n"
287 " ACCURACY: displacments of at most a few pixels.\n"
288 " *** This is the default method (without the :0).\n"
289 "\n"
290
291 #ifdef ALLOW_DFTIME
292 " -dftime[:0] Similar to the above, but after determining the\n"
293 " displacements, it modifies the images by using a\n"
294 " local fit in each voxel. The optional :0 has the\n"
295 " same meaning as for -dfspace.\n"
296 #if 0
297 " The optional :d means to also remove the mean and\n"
298 " linear trend from each pixel in the image time series.\n"
299 #endif
300 " ACCURACY: unevaluated. This option is intended\n"
301 " for FMRI use only!\n"
302 "\n"
303 " -dfspacetime[:0] Apply both algorithms: dfspace, then dftime.\n"
304 "\n"
305 #endif /* ALLOW_DFTIME */
306
307 " -cmass Initialize the translation estimate by aligning\n"
308 " the centers of mass of the images.\n"
309 " N.B.: The reported shifts from the registration algorithm\n"
310 " do NOT include the shifts due to this initial step.\n"
311 "\n"
312 "The new two options are used to play with the -dfspace algorithm,\n"
313 "which has a 'coarse' fit phase and a 'fine' fit phase:\n"
314 "\n"
315 " -fine blur dxy dphi Set the 3 'fine' fit parameters:\n"
316 " blur = FWHM of image blur prior to registration,\n"
317 " in pixels [must be > 0];\n"
318 " dxy = convergence tolerance for translations,\n"
319 " in pixels;\n"
320 " dphi = convergence tolerance for rotations,\n"
321 " in degrees.\n"
322 "\n"
323 " -nofine Turn off the 'fine' fit algorithm. By default, the\n"
324 " algorithm is on, with parameters 1.0, 0.07, 0.21.\n"
325 ) ;
326
327 return ;
328 }
329
330 /*---------------------------------------------------------------------*/
331
332 #define BASE_INPUT -1
333 #define BASE_AVER -2
334
REG_command_line(void)335 void REG_command_line(void)
336 {
337 int ii , nxbase=0 , nybase=0 , nerr , basecode ;
338 MRI_IMAGE * tim ;
339 MRI_IMARR * tarr ;
340
341 /*** look for options ***/
342
343 while( Iarg < Argc-2 && Argv[Iarg][0] == '-' ){
344
345 /** -cmass [12 Nov 2001] **/
346
347 if( strcmp(Argv[Iarg],"-cmass") == 0 ){
348 RG_use_cmass = 1 ;
349 Iarg++ ; continue ;
350 }
351
352 /** -nowrite **/
353
354 if( strncmp(Argv[Iarg],"-nowrite",5) == 0 ){
355 RG_skip_output = 1 ;
356 Iarg++ ; continue ;
357 }
358
359 if( strncmp(Argv[Iarg],"-keepsize",5) == 0 ){
360 RG_keepsize = 1 ;
361 Iarg++ ; continue ;
362 }
363
364 /** -debug **/
365
366 if( strncmp(Argv[Iarg],"-debug",5) == 0 ){
367 RG_debug = 1 ;
368 Iarg++ ; continue ;
369 }
370
371 /** -quiet **/
372
373 if( strncmp(Argv[Iarg],"-quiet",5) == 0 ){
374 RG_verbose = 0 ;
375 Iarg++ ; continue ;
376 }
377
378 /** -flim **/
379
380 if( strncmp(Argv[Iarg],"-flim",5) == 0 ){
381 RG_out_mode = MRI_float ;
382 Iarg++ ; continue ;
383 }
384
385 /** -wtim **/
386
387 if( strncmp(Argv[Iarg],"-wtim",5) == 0 ){
388 RG_imwt = mri_read_just_one( Argv[++Iarg] ) ;
389 Iarg++ ; continue ;
390 }
391
392 /** -prefix **/
393
394 if( strncmp(Argv[Iarg],"-prefix",5) == 0 ){
395 strcpy( RG_out_prefix , Argv[++Iarg] ) ;
396 ii = strlen(RG_out_prefix) ;
397 if( ii > 0 && RG_out_prefix[ii-1] != '.' ){
398 RG_out_prefix[ii] = '.' ;
399 RG_out_prefix[ii+1] = '\0' ;
400 }
401 Iarg++ ; continue ;
402 }
403
404 /** -dprefix **/
405
406 if( strncmp(Argv[Iarg],"-dprefix",5) == 0 ){
407 strcpy( RG_dout_prefix , Argv[++Iarg] ) ;
408 ii = strlen(RG_dout_prefix) ;
409 if( ii > 0 && RG_dout_prefix[ii-1] != '.' ){
410 RG_dout_prefix[ii] = '.' ;
411 RG_dout_prefix[ii+1] = '\0' ;
412 }
413 Iarg++ ; continue ;
414 }
415
416 /** -suffix **/
417
418 if( strncmp(Argv[Iarg],"-suffix",5) == 0 ){
419 Iarg++ ;
420 if( Argv[Iarg][0] == '.' ){
421 strcpy( RG_out_suffix , Argv[Iarg] ) ;
422 } else {
423 RG_out_suffix[0] = '.' ;
424 strcpy( RG_out_suffix + 1 , Argv[Iarg] ) ;
425 }
426 Iarg++ ; continue ;
427 }
428
429 /** -start **/
430
431 if( strncmp(Argv[Iarg],"-start",5) == 0 ){
432 char * ch ;
433 RG_out_start = strtol( Argv[++Iarg] , &ch , 10 ) ;
434 if( *ch != '\0' ){
435 fprintf(stderr,"** value after -start is illegal!\a\n") ;
436 exit(1) ;
437 }
438 Iarg++ ; continue ;
439 }
440
441 /** -step **/
442
443 if( strncmp(Argv[Iarg],"-step",5) == 0 ){
444 char * ch ;
445 RG_out_step = strtol( Argv[++Iarg] , &ch , 10 ) ;
446 if( *ch != '\0' ){
447 fprintf(stderr,"** value after -step is illegal!\a\n") ;
448 exit(1) ;
449 }
450 Iarg++ ; continue ;
451 }
452
453 /** -bilinear **/
454
455 if( strncmp(Argv[Iarg],"-bilinear",5) == 0 ){
456 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_BILINEAR ;
457 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
458 Iarg++ ; continue ;
459 }
460
461 if( strncmp(Argv[Iarg],"-Fourier",5) == 0 ){
462 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_FOURIER ;
463 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
464 Iarg++ ; continue ;
465 }
466
467 if( strncmp(Argv[Iarg],"-bicubic",5) == 0 ){
468 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_BICUBIC ;
469 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
470 Iarg++ ; continue ;
471 }
472
473 if( strncmp(Argv[Iarg],"-modes",5) == 0 ){ /* 1 Oct 1998 */
474 char * cpt ;
475
476 cpt = Argv[++Iarg] ;
477 if( strlen(cpt) > 2 ){
478 switch( cpt[2] ){
479 case 'l': RG_almode_coarse = MRI_BILINEAR ; break ;
480 case 'c': RG_almode_coarse = MRI_BICUBIC ; break ;
481 case 'u': RG_almode_coarse = MRI_FOURIER ; break ;
482 }
483 }
484
485 cpt = Argv[++Iarg] ;
486 if( strlen(cpt) > 2 ){
487 switch( cpt[2] ){
488 case 'l': RG_almode_fine = MRI_BILINEAR ; break ;
489 case 'c': RG_almode_fine = MRI_BICUBIC ; break ;
490 case 'u': RG_almode_fine = MRI_FOURIER ; break ;
491 }
492 }
493
494 cpt = Argv[++Iarg] ;
495 if( strlen(cpt) > 2 ){
496 switch( cpt[2] ){
497 case 'l': RG_almode_reg = MRI_BILINEAR ; break ;
498 case 'c': RG_almode_reg = MRI_BICUBIC ; break ;
499 case 'u': RG_almode_reg = MRI_FOURIER ; break ;
500 }
501 }
502
503 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
504 Iarg++ ; continue ;
505 }
506
507 if( strncmp(Argv[Iarg],"-mlcF",5) == 0 ){ /* 1 Oct 1998 */
508 RG_almode_coarse = MRI_BILINEAR ;
509 RG_almode_fine = MRI_BICUBIC ;
510 RG_almode_reg = MRI_FOURIER ;
511 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
512 Iarg++ ; continue ;
513 }
514
515 /** 05 Nov 1997: -params maxite sig dxy dph fsig fdxy fdph **/
516
517 if( strncmp(Argv[Iarg],"-params",6) == 0 ){
518 int maxite ;
519 float sig,dxy,dph,fsig,fdxy,fdph ;
520 maxite = strtol( Argv[++Iarg] , NULL , 10 ) ;
521 sig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
522 dxy = strtod( Argv[++Iarg] , NULL ) ;
523 dph = strtod( Argv[++Iarg] , NULL ) ;
524 fsig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
525 fdxy = strtod( Argv[++Iarg] , NULL ) ;
526 fdph = strtod( Argv[++Iarg] , NULL ) ;
527 #ifdef USE_2DALIGN
528 mri_2dalign_params( maxite,sig,dxy,dph,fsig,fdxy,fdph ) ;
529 #else
530 mri_align_params( maxite,sig,dxy,dph,fsig,fdxy,fdph ) ;
531 #endif
532 Iarg++ ; continue ;
533 }
534
535 if( strncmp(Argv[Iarg],"-nofine",6) == 0 ){
536 #ifdef USE_2DALIGN
537 mri_2dalign_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
538 #else
539 mri_align_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
540 #endif
541 Iarg++ ; continue ;
542 }
543
544 if( strncmp(Argv[Iarg],"-fine",4) == 0 ){
545 float fsig,fdxy,fdph ;
546 fsig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
547 fdxy = strtod( Argv[++Iarg] , NULL ) ;
548 fdph = strtod( Argv[++Iarg] , NULL ) ;
549 #ifdef USE_2DALIGN
550 mri_2dalign_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
551 #else
552 mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
553 #endif
554 Iarg++ ; continue ;
555 }
556
557 /** ALGORITHM: -dfspace **/
558
559 if( strstr(Argv[Iarg],"-dfspace") != NULL && strstr(Argv[Iarg],"-dfspacetime") == NULL ){
560 RG_meth = ALIGN_DFSPACE_TYPE ;
561 RG_methcode = 0 ;
562 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
563 Iarg++ ; continue ;
564 }
565
566 #ifdef ALLOW_DFTIME
567 /** -dftime **/
568
569 if( strstr(Argv[Iarg],"-dftime") != NULL ){
570 RG_meth = ALIGN_DFTIME_TYPE ;
571 RG_methcode = 0 ;
572 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
573 if( strstr(Argv[Iarg],":d") != NULL ) RG_methcode |= ALIGN_DETREND_CODE ;
574 Iarg++ ; continue ;
575 }
576
577 /** ALGORITHM: -dfspacetime **/
578
579 if( strstr(Argv[Iarg],"-dfspacetime") != NULL ){
580 RG_meth = ALIGN_DFTIME_TYPE ;
581 RG_methcode = ALIGN_DOBOTH_CODE ;
582 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
583 if( strstr(Argv[Iarg],":d") != NULL ) RG_methcode |= ALIGN_DETREND_CODE ;
584 Iarg++ ; continue ;
585 }
586 #endif
587
588 /** get to here is bad news **/
589
590 fprintf(stderr,"** Unknown option: %s\a\n",Argv[Iarg]) ;
591 exit(1) ;
592
593 }
594
595 if( RG_verbose ) RG_methcode |= ALIGN_VERBOSE_CODE ;
596 if( RG_debug ) RG_methcode |= ALIGN_DEBUG_CODE ;
597
598 /*** All options have been loaded. Next, read base_image, if available ***/
599
600 if( Iarg >= Argc ){
601 fprintf(stderr,"** No baseimage specified on command line!\a\n") ;
602 exit(1) ;
603 }
604
605 if( Argv[Iarg][0] != '+' ){
606 tim = mri_read_just_one( Argv[Iarg] ) ;
607 if( tim == NULL ){
608 fprintf(stderr,"** Can't read base_image -- end of run!\a\n") ; exit(1) ;
609 } else if ( ! MRI_IS_2D(tim) ){
610 fprintf(stderr,"** Base image is not 2D!\a\n") ; exit(1) ;
611 }
612 nxbase = tim->nx ;
613 nybase = tim->ny ;
614 if( tim->kind == MRI_float ) RG_baseimage = tim ;
615 else { RG_baseimage = mri_to_float(tim) ; mri_free(tim) ; }
616 basecode = BASE_INPUT ;
617 if( RG_verbose ) printf("-- read base image: file %s\n",Argv[Iarg]) ;
618 } else if( strcmp(Argv[Iarg],"+AVER")==0 || strcmp(Argv[Iarg],"+aver")==0 ){
619 basecode = BASE_AVER ;
620 if( RG_verbose ) printf("-- will set base image to AVER\n") ;
621 } else {
622 char * cp ;
623 basecode = strtol( Argv[Iarg]+1 , &cp , 10 ) ;
624 if( *cp != '\0' || basecode < 1 ){
625 fprintf(stderr,"** Can't interpret '+count' base_image input: %s\a\n",Argv[Iarg]) ;
626 exit(1) ;
627 }
628 if( RG_verbose ) printf("-- will set base image to input # %d\n",basecode) ;
629 }
630
631 /*** Read entire image sequence ***/
632
633 Iarg++ ;
634 if( Iarg >= Argc ){
635 fprintf(stderr,"** No image sequence specified on command line!\a\n") ;
636 exit(1) ;
637 }
638
639 INIT_IMARR(RG_imseq) ;
640
641 for( ; Iarg < Argc ; Iarg++ ){
642
643 tarr = mri_read_file( Argv[Iarg] ) ;
644
645 if( tarr == NULL || tarr->num == 0 ){
646 fprintf(stderr,
647 "** Can't read image(s) from file %s -- end of run!\a\n", Argv[Iarg]) ;
648 exit(1) ;
649 }
650
651 if( RG_out_mode < 0 ) RG_out_mode = tarr->imarr[0]->kind ;
652
653 for( ii=0 ; ii < tarr->num ; ii++ ){
654 if( ! MRI_IS_2D(tarr->imarr[ii]) ){
655 fprintf(stderr,"** Some input image is not 2D\a\n") ; exit(1) ;
656 }
657 if( tarr->imarr[ii]->kind == MRI_float ){
658 ADDTO_IMARR( RG_imseq , tarr->imarr[ii] ) ;
659 } else {
660 tim = mri_to_float( tarr->imarr[ii] ) ;
661 ADDTO_IMARR( RG_imseq , tim ) ;
662 mri_free( tarr->imarr[ii] ) ;
663 }
664 }
665
666 FREE_IMARR( tarr ) ; /* not DESTROY */
667 }
668
669 /*** if no base image, get dimensions from 1st sequence image ***/
670
671 if( RG_baseimage == NULL ){
672 if( RG_imseq->num <= 1 ){
673 fprintf(stderr,
674 "** No base_image supplied and only 1 image in sequence?\n"
675 "** Makes no sense! End of run!\a\n" ) ;
676 exit(1) ;
677 }
678 nxbase = RG_imseq->imarr[0]->nx ;
679 nybase = RG_imseq->imarr[0]->ny ;
680 }
681
682 /*** for each image in the sequence:
683 check for conformant dimensions ***/
684
685 nerr = 0 ;
686 for( ii=0 ; ii < RG_imseq->num ; ii++ ){
687 tim = RG_imseq->imarr[ii] ;
688 if( tim->nx != nxbase || tim->ny != nybase ){
689 fprintf(stderr,"** Image %d dimensions do not match base image!\a\n",ii+1) ;
690 nerr++ ;
691 }
692 }
693 if( nerr > 0 ) exit(1) ;
694
695 /*** if needed, create base image from image sequence ***/
696
697 if( RG_baseimage == NULL ){
698 if( basecode == BASE_AVER ){
699 register int pp , npix ;
700 register float * flar , * flfl ;
701 register float fac ;
702
703 if( RG_verbose ) printf("-- computing AVER image now\n") ;
704
705 RG_baseimage = mri_new( nxbase , nybase , MRI_float ) ;
706 flar = mri_data_pointer( RG_baseimage ) ;
707 npix = nxbase * nybase ;
708 for( pp=0 ; pp < npix ; pp++ ) { flar[pp] = 0.0 ; }
709
710 for( ii=0 ; ii < RG_imseq->num ; ii++ ){
711 flfl = mri_data_pointer( RG_imseq->imarr[ii] ) ;
712 for( pp=0 ; pp < npix ; pp++ ) flar[pp] += flfl[pp] ;
713 }
714
715 fac = 1.0 / RG_imseq->num ;
716 for( pp=0 ; pp < npix ; pp++ ) flar[pp] *= fac ;
717
718 } else if( basecode > 0 && basecode <= RG_imseq->num ){
719 if( RG_verbose ) printf("-- setting base image now\n") ;
720 RG_baseimage = mri_to_float( RG_imseq->imarr[basecode-1] ) ; /* copy it */
721 } else {
722 fprintf(stderr,"** Can't make baseimage as specified!\a\n") ;
723 exit(1) ;
724 }
725 }
726
727 /*** done (I hope) ***/
728
729 return ;
730 }
731