00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include <string.h>
00009
00010
00011
00012
00013
00014 static char RG_out_prefix[256] = "reg." ;
00015 static char RG_out_suffix[256] = "\0" ;
00016 static int RG_out_start = 1 ;
00017 static int RG_out_step = 1 ;
00018 static int RG_out_mode = -666 ;
00019
00020 static char RG_dout_prefix[256] = "\0" ;
00021
00022 static int RG_meth = ALIGN_DFSPACE_TYPE ;
00023 static int RG_methcode = 0 ;
00024 static int RG_verbose = 1 ;
00025 static int RG_skip_output = 0 ;
00026 static int RG_debug = 0 ;
00027
00028 static int RG_almode_coarse = MRI_BICUBIC ;
00029 static int RG_almode_fine = MRI_BICUBIC ;
00030 static int RG_almode_reg = MRI_BICUBIC ;
00031
00032 static int RG_use_cmass = 0 ;
00033
00034 static MRI_IMAGE * RG_imwt = NULL ;
00035
00036 static int Iarg = 1 ;
00037 static int Argc ;
00038 static char ** Argv ;
00039
00040 static MRI_IMAGE * RG_baseimage = NULL ;
00041 static MRI_IMARR * RG_imseq = NULL ;
00042
00043
00044
00045 void REG_syntax(void) ;
00046 void REG_command_line(void) ;
00047
00048
00049
00050
00051
00052 #undef USE_2DALIGN
00053 #ifdef USE_2DALIGN
00054 # undef ALLOW_DFTIME
00055 #endif
00056
00057 int main( int argc , char *argv[] )
00058 {
00059 MRI_IMAGE * qim , * tim ;
00060 MRI_IMARR * regar = NULL ;
00061 float * dx , * dy , * phi ;
00062 int kim , imcount ;
00063 char fname[666] ;
00064 float dxtop , dytop , phitop ;
00065 FILE * dxfil , * dyfil , * phifil ;
00066
00067
00068
00069 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){ REG_syntax() ; exit(0); }
00070
00071 machdep() ;
00072
00073 Argc = argc ;
00074 Argv = argv ;
00075 Iarg = 1 ;
00076
00077 REG_command_line() ;
00078
00079 imcount = RG_imseq->num ;
00080 dx = (float *) malloc( sizeof(float) * imcount ) ;
00081 dy = (float *) malloc( sizeof(float) * imcount ) ;
00082 phi = (float *) malloc( sizeof(float) * imcount ) ;
00083 if( dx == NULL || dy == NULL || phi == NULL ){
00084 fprintf(stderr,"** malloc failure for mri_rotate parameters!\a\n") ; exit(1) ;
00085 }
00086
00087
00088
00089 if( RG_use_cmass ){
00090 float xbase,ybase , xcm,ycm ;
00091 int ii,jj,kk , di,dj , nx,ny , sj,si ;
00092 MRI_IMAGE *fim , *nim ;
00093 float *far , *nar ;
00094
00095 if( RG_verbose ) printf("-- doing -cmass adjustments\n") ;
00096
00097 mri_get_cmass_2D( RG_baseimage , &xbase , &ybase ) ;
00098
00099 for( kk=0 ; kk < imcount ; kk++ ){
00100 fim = mri_to_float( IMARR_SUBIM(RG_imseq,kk) ) ;
00101 mri_get_cmass_2D( fim , &xcm , &ycm ) ;
00102 di = (int) rint(xbase-xcm) ;
00103 dj = (int) rint(ybase-ycm) ;
00104 if( di != 0 || dj != 0 ){
00105
00106 #if 0
00107 fprintf(stderr,"Image %d: xbase=%g ybase=%g xcm=%g ycm=%g di=%d dj=%d\n",
00108 kk,xbase,ybase,xcm,ycm,di,dj ) ;
00109 #endif
00110
00111 nim = mri_new_conforming( fim , MRI_float );
00112 far = MRI_FLOAT_PTR(fim) ;
00113 nar = MRI_FLOAT_PTR(nim) ;
00114 nx = fim->nx ; ny = fim->ny ;
00115 for( jj=0 ; jj < ny ; jj++ ){
00116 sj = jj + dj ;
00117 if( sj < 0 || sj >= ny ) continue ;
00118 for( ii=0 ; ii < nx ; ii++ ){
00119 si = ii + di ;
00120 if( si < 0 || si >= nx ) continue ;
00121 nar[si+nx*sj] = far[ii+nx*jj] ;
00122 }
00123 }
00124 mri_free( IMARR_SUBIM(RG_imseq,kk) ) ;
00125 IMARR_SUBIM(RG_imseq,kk) = nim ;
00126 }
00127 mri_free(fim) ;
00128 }
00129 }
00130
00131
00132
00133 if( RG_verbose ) printf("-- beginning alignment\n") ;
00134 switch( RG_meth ){
00135
00136 default:
00137 case ALIGN_DFSPACE_TYPE:
00138 #ifdef USE_2DALIGN
00139 regar = mri_2dalign_many( RG_baseimage,RG_imwt, RG_imseq, dx,dy,phi ) ;
00140 #else
00141 if( ! RG_skip_output ) RG_methcode |= ALIGN_REGISTER_CODE ;
00142 regar = mri_align_dfspace( RG_baseimage,RG_imwt, RG_imseq, RG_methcode, dx,dy,phi ) ;
00143 #endif
00144 break ;
00145
00146 #ifdef ALLOW_DFTIME
00147 case ALIGN_DFTIME_TYPE:
00148 if( ! RG_skip_output ) RG_methcode |= ALIGN_REGISTER_CODE ;
00149 regar = mri_align_dftime( RG_baseimage,RG_imwt, RG_imseq, RG_methcode, dx,dy,phi ) ;
00150 break ;
00151 #endif
00152 }
00153
00154 dxtop = dytop = phitop = 0.0 ;
00155
00156 if( strlen(RG_dout_prefix) > 0 ){
00157 sprintf( fname , "%s%s" , RG_dout_prefix , "dx" ) ; dxfil = fopen( fname , "w" ) ;
00158 sprintf( fname , "%s%s" , RG_dout_prefix , "dy" ) ; dyfil = fopen( fname , "w" ) ;
00159 sprintf( fname , "%s%s" , RG_dout_prefix , "phi" ) ; phifil = fopen( fname , "w" ) ;
00160 } else {
00161 dxfil = dyfil = phifil = NULL ;
00162 }
00163
00164 for( kim=0 ; kim < imcount ; kim++ ){
00165
00166 if( fabs(dx[kim]) > fabs(dxtop) ) dxtop = dx[kim] ;
00167 if( fabs(dy[kim]) > fabs(dytop) ) dytop = dy[kim] ;
00168 if( fabs(phi[kim]) > fabs(phitop) ) phitop = phi[kim] ;
00169
00170 if( ! RG_skip_output ){
00171 sprintf( fname , "%s%04d%s" , RG_out_prefix ,
00172 RG_out_start + kim*RG_out_step , RG_out_suffix ) ;
00173
00174 if( regar == NULL ){
00175 tim = mri_rota_variable( RG_almode_reg ,
00176 RG_imseq->imarr[kim] ,
00177 dx[kim],dy[kim],phi[kim] ) ;
00178 } else {
00179 tim = regar->imarr[kim] ;
00180 }
00181 } else {
00182 sprintf( fname , "%4d" , kim+1 ) ;
00183 }
00184
00185 if( RG_verbose )
00186 printf("-- image %s: dx = %6.3f dy = %6.3f phi = %6.3f\n" ,
00187 fname,dx[kim],dy[kim],(180.0/PI)*phi[kim] ) ;
00188
00189 if( dxfil != NULL ) fprintf( dxfil , "%6.3f\n" , dx[kim] ) ;
00190 if( dyfil != NULL ) fprintf( dyfil , "%6.3f\n" , dy[kim] ) ;
00191 if( phifil != NULL ) fprintf( phifil , "%6.3f\n" , (180.0/PI)*phi[kim] ) ;
00192
00193 if( ! RG_skip_output ){
00194 switch( RG_out_mode ){
00195
00196 default:
00197 mri_write( fname , tim ) ;
00198 break ;
00199
00200 case MRI_short:
00201 qim = mri_to_short( 1.0 , tim ) ;
00202 mri_write( fname , qim ) ; mri_free( qim ) ;
00203 break ;
00204
00205 case MRI_byte:
00206 qim = mri_to_byte( tim ) ;
00207 mri_write( fname , qim ) ; mri_free( qim ) ;
00208 break ;
00209 }
00210 mri_free( tim ) ; mri_free( RG_imseq->imarr[kim] ) ;
00211 }
00212 }
00213
00214 if( RG_verbose )
00215 printf("-- MAX: %*s dx = %6.3f dy = %6.3f phi = %6.3f\n" ,
00216 (int)strlen(fname) , " " , dxtop,dytop,phitop*(180.0/PI) ) ;
00217
00218 exit(0) ;
00219 }
00220
00221
00222
00223 void REG_syntax(void)
00224 {
00225 printf(
00226 "Usage: imreg [options] base_image image_sequence ...\n"
00227 " * Registers each 2D image in 'image_sequence' to 'base_image'.\n"
00228 " * If 'base_image' = '+AVER', will compute the base image as\n"
00229 " the average of the images in 'image_sequence'.\n"
00230 " * If 'base_image' = '+count', will use the count-th image in the\n"
00231 " sequence as the base image. Here, count is 1,2,3, ....\n"
00232 "\n"
00233 "OUTPUT OPTIONS:\n"
00234 " -nowrite Don't write outputs, just print progress reports.\n"
00235 " -prefix pname The output files will be named in the format\n"
00236 " -suffix sname 'pname.index.sname' where 'pname' and 'sname'\n"
00237 " -start si are strings given by the first 2 options.\n"
00238 " -step ss 'index' is a number, given by 'si+(i-1)*ss'\n"
00239 " for the i-th output file, for i=1,2,...\n"
00240 " *** Default pname = 'reg.'\n"
00241 " *** Default sname = nothing at all\n"
00242 " *** Default si = 1\n"
00243 " *** Default ss = 1\n"
00244 "\n"
00245 " -flim Write output in mrilib floating point format\n"
00246 " (which can be converted to shorts using program ftosh).\n"
00247 " *** Default is to write images in format of first\n"
00248 " input file in the image_sequence.\n"
00249 "\n"
00250 " -quiet Don't write progress report messages.\n"
00251 " -debug Write lots of debugging output!\n"
00252 "\n"
00253 " -dprefix dname Write files 'dname'.dx, 'dname'.dy, 'dname'.phi\n"
00254 " for use in time series analysis.\n"
00255 "\n"
00256 "ALIGNMENT ALGORITHMS:\n"
00257 " -bilinear Uses bilinear interpolation during the iterative\n"
00258 " adjustment procedure, rather than the default\n"
00259 " bicubic interpolation. NOT RECOMMENDED!\n"
00260 " -modes c f r Uses interpolation modes 'c', 'f', and 'r' during\n"
00261 " the coarse, fine, and registration phases of the\n"
00262 " algorithm, respectively. The modes can be selected\n"
00263 " from 'bilinear', 'bicubic', and 'Fourier'. The\n"
00264 " default is '-modes bicubic bicubic bicubic'.\n"
00265 " -mlcF Equivalent to '-modes bilinear bicubic Fourier'.\n"
00266 "\n"
00267 " -wtim filename Uses the image in 'filename' as a weighting factor\n"
00268 " for each voxel (the larger the value the more\n"
00269 " importance is given to that voxel).\n"
00270 "\n"
00271 " -dfspace[:0] Uses the 'iterated diffential spatial' method to\n"
00272 " align the images. The optional :0 indicates to\n"
00273 " skip the iteration of the method, and to use the\n"
00274 " simpler linear differential spatial alignment method.\n"
00275 " ACCURACY: displacments of at most a few pixels.\n"
00276 " *** This is the default method (without the :0).\n"
00277 "\n"
00278
00279 #ifdef ALLOW_DFTIME
00280 " -dftime[:0] Similar to the above, but after determining the\n"
00281 " displacements, it modifies the images by using a\n"
00282 " local fit in each voxel. The optional :0 has the\n"
00283 " same meaning as for -dfspace.\n"
00284 #if 0
00285 " The optional :d means to also remove the mean and\n"
00286 " linear trend from each pixel in the image time series.\n"
00287 #endif
00288 " ACCURACY: unevaluated. This option is intended\n"
00289 " for FMRI use only!\n"
00290 "\n"
00291 " -dfspacetime[:0] Apply both algorithms: dfspace, then dftime.\n"
00292 "\n"
00293 #endif
00294
00295 " -cmass Initialize the translation estimate by aligning\n"
00296 " the centers of mass of the images.\n"
00297 " N.B.: The reported shifts from the registration algorithm\n"
00298 " do NOT include the shifts due to this initial step.\n"
00299 "\n"
00300 "The new two options are used to play with the -dfspace algorithm,\n"
00301 "which has a 'coarse' fit phase and a 'fine' fit phase:\n"
00302 "\n"
00303 " -fine blur dxy dphi Set the 3 'fine' fit parameters:\n"
00304 " blur = FWHM of image blur prior to registration,\n"
00305 " in pixels [must be > 0];\n"
00306 " dxy = convergence tolerance for translations,\n"
00307 " in pixels;\n"
00308 " dphi = convergence tolerance for rotations,\n"
00309 " in degrees.\n"
00310 "\n"
00311 " -nofine Turn off the 'fine' fit algorithm. By default, the\n"
00312 " algorithm is on, with parameters 1.0, 0.07, 0.21.\n"
00313 ) ;
00314
00315 return ;
00316 }
00317
00318
00319
00320 #define BASE_INPUT -1
00321 #define BASE_AVER -2
00322
00323 void REG_command_line(void)
00324 {
00325 int ii , nxbase , nybase , nerr , basecode ;
00326 MRI_IMAGE * tim ;
00327 MRI_IMARR * tarr ;
00328
00329
00330
00331 while( Iarg < Argc-2 && Argv[Iarg][0] == '-' ){
00332
00333
00334
00335 if( strcmp(Argv[Iarg],"-cmass") == 0 ){
00336 RG_use_cmass = 1 ;
00337 Iarg++ ; continue ;
00338 }
00339
00340
00341
00342 if( strncmp(Argv[Iarg],"-nowrite",5) == 0 ){
00343 RG_skip_output = 1 ;
00344 Iarg++ ; continue ;
00345 }
00346
00347
00348
00349 if( strncmp(Argv[Iarg],"-debug",5) == 0 ){
00350 RG_debug = 1 ;
00351 Iarg++ ; continue ;
00352 }
00353
00354
00355
00356 if( strncmp(Argv[Iarg],"-quiet",5) == 0 ){
00357 RG_verbose = 0 ;
00358 Iarg++ ; continue ;
00359 }
00360
00361
00362
00363 if( strncmp(Argv[Iarg],"-flim",5) == 0 ){
00364 RG_out_mode = MRI_float ;
00365 Iarg++ ; continue ;
00366 }
00367
00368
00369
00370 if( strncmp(Argv[Iarg],"-wtim",5) == 0 ){
00371 RG_imwt = mri_read_just_one( Argv[++Iarg] ) ;
00372 Iarg++ ; continue ;
00373 }
00374
00375
00376
00377 if( strncmp(Argv[Iarg],"-prefix",5) == 0 ){
00378 strcpy( RG_out_prefix , Argv[++Iarg] ) ;
00379 ii = strlen(RG_out_prefix) ;
00380 if( ii > 0 && RG_out_prefix[ii-1] != '.' ){
00381 RG_out_prefix[ii] = '.' ;
00382 RG_out_prefix[ii+1] = '\0' ;
00383 }
00384 Iarg++ ; continue ;
00385 }
00386
00387
00388
00389 if( strncmp(Argv[Iarg],"-dprefix",5) == 0 ){
00390 strcpy( RG_dout_prefix , Argv[++Iarg] ) ;
00391 ii = strlen(RG_dout_prefix) ;
00392 if( ii > 0 && RG_dout_prefix[ii-1] != '.' ){
00393 RG_dout_prefix[ii] = '.' ;
00394 RG_dout_prefix[ii+1] = '\0' ;
00395 }
00396 Iarg++ ; continue ;
00397 }
00398
00399
00400
00401 if( strncmp(Argv[Iarg],"-suffix",5) == 0 ){
00402 Iarg++ ;
00403 if( Argv[Iarg][0] == '.' ){
00404 strcpy( RG_out_suffix , Argv[Iarg] ) ;
00405 } else {
00406 RG_out_suffix[0] = '.' ;
00407 strcpy( RG_out_suffix + 1 , Argv[Iarg] ) ;
00408 }
00409 Iarg++ ; continue ;
00410 }
00411
00412
00413
00414 if( strncmp(Argv[Iarg],"-start",5) == 0 ){
00415 char * ch ;
00416 RG_out_start = strtol( Argv[++Iarg] , &ch , 10 ) ;
00417 if( *ch != '\0' ){
00418 fprintf(stderr,"** value after -start is illegal!\a\n") ;
00419 exit(1) ;
00420 }
00421 Iarg++ ; continue ;
00422 }
00423
00424
00425
00426 if( strncmp(Argv[Iarg],"-step",5) == 0 ){
00427 char * ch ;
00428 RG_out_step = strtol( Argv[++Iarg] , &ch , 10 ) ;
00429 if( *ch != '\0' ){
00430 fprintf(stderr,"** value after -step is illegal!\a\n") ;
00431 exit(1) ;
00432 }
00433 Iarg++ ; continue ;
00434 }
00435
00436
00437
00438 if( strncmp(Argv[Iarg],"-bilinear",5) == 0 ){
00439 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_BILINEAR ;
00440 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
00441 Iarg++ ; continue ;
00442 }
00443
00444 if( strncmp(Argv[Iarg],"-Fourier",5) == 0 ){
00445 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_FOURIER ;
00446 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
00447 Iarg++ ; continue ;
00448 }
00449
00450 if( strncmp(Argv[Iarg],"-bicubic",5) == 0 ){
00451 RG_almode_coarse = RG_almode_fine = RG_almode_reg = MRI_BICUBIC ;
00452 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
00453 Iarg++ ; continue ;
00454 }
00455
00456 if( strncmp(Argv[Iarg],"-modes",5) == 0 ){
00457 char * cpt ;
00458
00459 cpt = Argv[++Iarg] ;
00460 if( strlen(cpt) > 2 ){
00461 switch( cpt[2] ){
00462 case 'l': RG_almode_coarse = MRI_BILINEAR ; break ;
00463 case 'c': RG_almode_coarse = MRI_BICUBIC ; break ;
00464 case 'u': RG_almode_coarse = MRI_FOURIER ; break ;
00465 }
00466 }
00467
00468 cpt = Argv[++Iarg] ;
00469 if( strlen(cpt) > 2 ){
00470 switch( cpt[2] ){
00471 case 'l': RG_almode_fine = MRI_BILINEAR ; break ;
00472 case 'c': RG_almode_fine = MRI_BICUBIC ; break ;
00473 case 'u': RG_almode_fine = MRI_FOURIER ; break ;
00474 }
00475 }
00476
00477 cpt = Argv[++Iarg] ;
00478 if( strlen(cpt) > 2 ){
00479 switch( cpt[2] ){
00480 case 'l': RG_almode_reg = MRI_BILINEAR ; break ;
00481 case 'c': RG_almode_reg = MRI_BICUBIC ; break ;
00482 case 'u': RG_almode_reg = MRI_FOURIER ; break ;
00483 }
00484 }
00485
00486 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
00487 Iarg++ ; continue ;
00488 }
00489
00490 if( strncmp(Argv[Iarg],"-mlcF",5) == 0 ){
00491 RG_almode_coarse = MRI_BILINEAR ;
00492 RG_almode_fine = MRI_BICUBIC ;
00493 RG_almode_reg = MRI_FOURIER ;
00494 mri_align_method( RG_almode_coarse,RG_almode_fine,RG_almode_reg ) ;
00495 Iarg++ ; continue ;
00496 }
00497
00498
00499
00500 if( strncmp(Argv[Iarg],"-params",6) == 0 ){
00501 int maxite ;
00502 float sig,dxy,dph,fsig,fdxy,fdph ;
00503 maxite = strtol( Argv[++Iarg] , NULL , 10 ) ;
00504 sig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
00505 dxy = strtod( Argv[++Iarg] , NULL ) ;
00506 dph = strtod( Argv[++Iarg] , NULL ) ;
00507 fsig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
00508 fdxy = strtod( Argv[++Iarg] , NULL ) ;
00509 fdph = strtod( Argv[++Iarg] , NULL ) ;
00510 #ifdef USE_2DALIGN
00511 mri_2dalign_params( maxite,sig,dxy,dph,fsig,fdxy,fdph ) ;
00512 #else
00513 mri_align_params( maxite,sig,dxy,dph,fsig,fdxy,fdph ) ;
00514 #endif
00515 Iarg++ ; continue ;
00516 }
00517
00518 if( strncmp(Argv[Iarg],"-nofine",6) == 0 ){
00519 #ifdef USE_2DALIGN
00520 mri_2dalign_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
00521 #else
00522 mri_align_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
00523 #endif
00524 Iarg++ ; continue ;
00525 }
00526
00527 if( strncmp(Argv[Iarg],"-fine",4) == 0 ){
00528 float fsig,fdxy,fdph ;
00529 fsig = strtod( Argv[++Iarg] , NULL ) * 0.42466090 ;
00530 fdxy = strtod( Argv[++Iarg] , NULL ) ;
00531 fdph = strtod( Argv[++Iarg] , NULL ) ;
00532 #ifdef USE_2DALIGN
00533 mri_2dalign_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
00534 #else
00535 mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
00536 #endif
00537 Iarg++ ; continue ;
00538 }
00539
00540
00541
00542 if( strstr(Argv[Iarg],"-dfspace") != NULL && strstr(Argv[Iarg],"-dfspacetime") == NULL ){
00543 RG_meth = ALIGN_DFSPACE_TYPE ;
00544 RG_methcode = 0 ;
00545 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
00546 Iarg++ ; continue ;
00547 }
00548
00549 #ifdef ALLOW_DFTIME
00550
00551
00552 if( strstr(Argv[Iarg],"-dftime") != NULL ){
00553 RG_meth = ALIGN_DFTIME_TYPE ;
00554 RG_methcode = 0 ;
00555 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
00556 if( strstr(Argv[Iarg],":d") != NULL ) RG_methcode |= ALIGN_DETREND_CODE ;
00557 Iarg++ ; continue ;
00558 }
00559
00560
00561
00562 if( strstr(Argv[Iarg],"-dfspacetime") != NULL ){
00563 RG_meth = ALIGN_DFTIME_TYPE ;
00564 RG_methcode = ALIGN_DOBOTH_CODE ;
00565 if( strstr(Argv[Iarg],":0") != NULL ) RG_methcode |= ALIGN_NOITER_CODE ;
00566 if( strstr(Argv[Iarg],":d") != NULL ) RG_methcode |= ALIGN_DETREND_CODE ;
00567 Iarg++ ; continue ;
00568 }
00569 #endif
00570
00571
00572
00573 fprintf(stderr,"** Unknown option: %s\a\n",Argv[Iarg]) ;
00574 exit(1) ;
00575
00576 }
00577
00578 if( RG_verbose ) RG_methcode |= ALIGN_VERBOSE_CODE ;
00579 if( RG_debug ) RG_methcode |= ALIGN_DEBUG_CODE ;
00580
00581
00582
00583 if( Iarg >= Argc ){
00584 fprintf(stderr,"** No baseimage specified on command line!\a\n") ;
00585 exit(1) ;
00586 }
00587
00588 if( Argv[Iarg][0] != '+' ){
00589 tim = mri_read_just_one( Argv[Iarg] ) ;
00590 if( tim == NULL ){
00591 fprintf(stderr,"** Can't read base_image -- end of run!\a\n") ; exit(1) ;
00592 } else if ( ! MRI_IS_2D(tim) ){
00593 fprintf(stderr,"** Base image is not 2D!\a\n") ; exit(1) ;
00594 }
00595 nxbase = tim->nx ;
00596 nybase = tim->ny ;
00597 if( tim->kind == MRI_float ) RG_baseimage = tim ;
00598 else { RG_baseimage = mri_to_float(tim) ; mri_free(tim) ; }
00599 basecode = BASE_INPUT ;
00600 if( RG_verbose ) printf("-- read base image: file %s\n",Argv[Iarg]) ;
00601 } else if( strcmp(Argv[Iarg],"+AVER")==0 || strcmp(Argv[Iarg],"+aver")==0 ){
00602 basecode = BASE_AVER ;
00603 if( RG_verbose ) printf("-- will set base image to AVER\n") ;
00604 } else {
00605 char * cp ;
00606 basecode = strtol( Argv[Iarg]+1 , &cp , 10 ) ;
00607 if( *cp != '\0' || basecode < 1 ){
00608 fprintf(stderr,"** Can't interpret '+count' base_image input: %s\a\n",Argv[Iarg]) ;
00609 exit(1) ;
00610 }
00611 if( RG_verbose ) printf("-- will set base image to input # %d\n",basecode) ;
00612 }
00613
00614
00615
00616 Iarg++ ;
00617 if( Iarg >= Argc ){
00618 fprintf(stderr,"** No image sequence specified on command line!\a\n") ;
00619 exit(1) ;
00620 }
00621
00622 INIT_IMARR(RG_imseq) ;
00623
00624 for( ; Iarg < Argc ; Iarg++ ){
00625
00626 tarr = mri_read_file( Argv[Iarg] ) ;
00627
00628 if( tarr == NULL || tarr->num == 0 ){
00629 fprintf(stderr,
00630 "** Can't read image(s) from file %s -- end of run!\a\n", Argv[Iarg]) ;
00631 exit(1) ;
00632 }
00633
00634 if( RG_out_mode < 0 ) RG_out_mode = tarr->imarr[0]->kind ;
00635
00636 for( ii=0 ; ii < tarr->num ; ii++ ){
00637 if( ! MRI_IS_2D(tarr->imarr[ii]) ){
00638 fprintf(stderr,"** Some input image is not 2D\a\n") ; exit(1) ;
00639 }
00640 if( tarr->imarr[ii]->kind == MRI_float ){
00641 ADDTO_IMARR( RG_imseq , tarr->imarr[ii] ) ;
00642 } else {
00643 tim = mri_to_float( tarr->imarr[ii] ) ;
00644 ADDTO_IMARR( RG_imseq , tim ) ;
00645 mri_free( tarr->imarr[ii] ) ;
00646 }
00647 }
00648
00649 FREE_IMARR( tarr ) ;
00650 }
00651
00652
00653
00654 if( RG_baseimage == NULL ){
00655 if( RG_imseq->num <= 1 ){
00656 fprintf(stderr,
00657 "** No base_image supplied and only 1 image in sequence?\n"
00658 "** Makes no sense! End of run!\a\n" ) ;
00659 exit(1) ;
00660 }
00661 nxbase = RG_imseq->imarr[0]->nx ;
00662 nybase = RG_imseq->imarr[0]->ny ;
00663 }
00664
00665
00666
00667
00668 nerr = 0 ;
00669 for( ii=0 ; ii < RG_imseq->num ; ii++ ){
00670 tim = RG_imseq->imarr[ii] ;
00671 if( tim->nx != nxbase || tim->ny != nybase ){
00672 fprintf(stderr,"** Image %d dimensions do not match base image!\a\n",ii+1) ;
00673 nerr++ ;
00674 }
00675 }
00676 if( nerr > 0 ) exit(1) ;
00677
00678
00679
00680 if( RG_baseimage == NULL ){
00681 if( basecode == BASE_AVER ){
00682 register int pp , npix ;
00683 register float * flar , * flfl ;
00684 register float fac ;
00685
00686 if( RG_verbose ) printf("-- computing AVER image now\n") ;
00687
00688 RG_baseimage = mri_new( nxbase , nybase , MRI_float ) ;
00689 flar = mri_data_pointer( RG_baseimage ) ;
00690 npix = nxbase * nybase ;
00691 for( pp=0 ; pp < npix ; pp++ ) { flar[pp] = 0.0 ; }
00692
00693 for( ii=0 ; ii < RG_imseq->num ; ii++ ){
00694 flfl = mri_data_pointer( RG_imseq->imarr[ii] ) ;
00695 for( pp=0 ; pp < npix ; pp++ ) flar[pp] += flfl[pp] ;
00696 }
00697
00698 fac = 1.0 / RG_imseq->num ;
00699 for( pp=0 ; pp < npix ; pp++ ) flar[pp] *= fac ;
00700
00701 } else if( basecode > 0 && basecode <= RG_imseq->num ){
00702 if( RG_verbose ) printf("-- setting base image now\n") ;
00703 RG_baseimage = mri_to_float( RG_imseq->imarr[basecode-1] ) ;
00704 } else {
00705 fprintf(stderr,"** Can't make baseimage as specified!\a\n") ;
00706 exit(1) ;
00707 }
00708 }
00709
00710
00711
00712 return ;
00713 }
