00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include <string.h>
00009
00010 #include "thd_shear3d.h"
00011
00012 #define ERREX(str) (fprintf(stderr,"** %s\n",str),exit(1))
00013
00014
00015
00016
00017
00018 static int Iarg = 1 ;
00019 static int Argc ;
00020 static char ** Argv ;
00021
00022 static int VL_nbase = 0 ;
00023 static int VL_intern = 1 ;
00024 static int VL_resam = MRI_FOURIER ;
00025 static int VL_final = -1 ;
00026 static int VL_clipit = 1 ;
00027 static MRI_IMAGE * VL_imbase = NULL ;
00028 static MRI_IMAGE * VL_imwt = NULL ;
00029 static int VL_wtinp = 0 ;
00030
00031 static int VL_zpad = 0 ;
00032
00033 static int VL_twopass= 0 ;
00034 static float VL_twoblur= 2.0 ;
00035 static int VL_twosum = 1 ;
00036 static int VL_twodup = 0 ;
00037 static float VL_bxorg, VL_byorg, VL_bzorg ;
00038
00039 static float VL_edging ;
00040 static int VL_edperc=-1 ;
00041
00042 static int VL_coarse_del=10 ;
00043 static int VL_coarse_num=2 ;
00044
00045 static THD_3dim_dataset * VL_dset = NULL ;
00046 static THD_3dim_dataset * VL_bset = NULL ;
00047
00048 static char VL_prefix[256] = "volreg" ;
00049 static int VL_verbose = 0 ;
00050 static char VL_dfile[256] = "\0" ;
00051 static char VL_1Dfile[256] = "\0" ;
00052
00053 static int VL_maxite = 19 ;
00054 static float VL_dxy = 0.02;
00055 static float VL_dph = 0.03 ;
00056 static float VL_del = 0.70 ;
00057
00058 static int VL_rotcom = 0 ;
00059
00060 static THD_3dim_dataset *VL_rotpar_dset =NULL ,
00061 *VL_gridpar_dset=NULL ;
00062
00063 static int VL_tshift = 0 ,
00064 VL_tshift_ignore = 0 ;
00065
00066 static int VL_sinit = 1 ;
00067
00068
00069
00070 void VL_syntax(void) ;
00071 void VL_command_line(void) ;
00072
00073 float voldif( int nx, int ny, int nz, float *b,
00074 int dx, int dy, int dz, float *v, int edge ) ;
00075
00076 void get_best_shift( int nx, int ny, int nz,
00077 float *b, float *v, int *dxp,int *dyp,int *dzp ) ;
00078
00079
00080
00081
00082 int main( int argc , char *argv[] )
00083 {
00084 MRI_3dalign_basis * albase ;
00085 THD_3dim_dataset * new_dset ;
00086 MRI_IMAGE * qim , * tim , * fim ;
00087 double cputim ;
00088 float *dx, *dy, *dz, *roll, *yaw, *pitch, *rmsnew, *rmsold, *imb, *tar ;
00089 float ddx,ddy,ddz , sum ;
00090 float dxtop,dytop,dztop , rolltop,yawtop,pitchtop ;
00091 float dxbot,dybot,dzbot , rollbot,yawbot,pitchbot ;
00092 float dxbar,dybar,dzbar , rollbar,yawbar,pitchbar ;
00093 int kim,ii , imcount , iha , ax1,ax2,ax3 , hax1,hax2,hax3 ;
00094
00095 float *dx_1,*dy_1,*dz_1, *roll_1,*yaw_1,*pitch_1 ;
00096 int nx,ny,nz ;
00097
00098 static char * modes[] = {
00099 "-NN" , "-linear" , "-cubic" , "-Fourier" , "-quintic" , "-heptic" } ;
00100
00101 #define MATVEC_DICOM 1
00102 #define MATVEC_ORDER 2
00103
00104 int matvec=0 ;
00105 THD_dmat33 rmat , pp,ppt ;
00106 THD_dfvec3 tvec ;
00107 int npad_neg=0 ,
00108 npad_pos=0 , npadd=0 ;
00109
00110
00111
00112 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ){ VL_syntax() ; exit(0); }
00113
00114 mainENTRY("3dvolreg main") ; machdep() ; AFNI_logger("3dvolreg",argc,argv) ;
00115 PRINT_VERSION("3dvolreg") ;
00116
00117
00118
00119 { int new_argc ; char ** new_argv ;
00120 addto_args( argc , argv , &new_argc , &new_argv ) ;
00121 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00122 }
00123
00124
00125 Argc = argc ; Argv = argv ; Iarg = 1 ;
00126 VL_command_line() ;
00127
00128 mri_3dalign_wtrimming(1) ;
00129
00130
00131
00132 imcount = DSET_NVALS( VL_dset ) ;
00133 dx = (float *) malloc( sizeof(float) * imcount ) ;
00134 dy = (float *) malloc( sizeof(float) * imcount ) ;
00135 dz = (float *) malloc( sizeof(float) * imcount ) ;
00136 roll = (float *) malloc( sizeof(float) * imcount ) ;
00137 pitch = (float *) malloc( sizeof(float) * imcount ) ;
00138 yaw = (float *) malloc( sizeof(float) * imcount ) ;
00139 rmsnew = (float *) malloc( sizeof(float) * imcount ) ;
00140 rmsold = (float *) malloc( sizeof(float) * imcount ) ;
00141
00142 iha = THD_handedness( VL_dset ) ;
00143 ax1 = THD_axcode( VL_dset , 'I' ) ; hax1 = ax1 * iha ;
00144 ax2 = THD_axcode( VL_dset , 'R' ) ; hax2 = ax2 * iha ;
00145 ax3 = THD_axcode( VL_dset , 'A' ) ; hax3 = ax3 * iha ;
00146
00147
00148
00149 new_dset = EDIT_empty_copy( VL_dset ) ;
00150
00151
00152
00153 if( VL_gridpar_dset != NULL ){
00154 int mm , nz_gp , nz_ds ;
00155
00156
00157
00158 mm = THD_dataset_mismatch( VL_gridpar_dset , VL_dset ) ;
00159 if( mm & (MISMATCH_DELTA | MISMATCH_ORIENT) ){
00160 fprintf(stderr,"** Fatal Error:\n"
00161 "** -gridparent dataset and input dataset don't\n"
00162 "** match in grid spacing and/or orientation!\n" ) ;
00163 exit(1) ;
00164 }
00165
00166 if( DSET_NX(VL_gridpar_dset) != DSET_NX(VL_dset) ||
00167 DSET_NY(VL_gridpar_dset) != DSET_NY(VL_dset) ){
00168
00169 fprintf(stderr,"** Fatal Error:\n"
00170 "** -gridparent and input datasets\n"
00171 "** don't match in x,y dimensions!\n" ) ;
00172 exit(1) ;
00173 }
00174
00175
00176
00177 nz_gp = DSET_NZ(VL_gridpar_dset) ; nz_ds = DSET_NZ(VL_dset) ;
00178
00179 if( nz_gp < nz_ds ){
00180 fprintf(stderr,"** Fatal Error:\n"
00181 "** -gridparent has fewer slices than input dataset!\n") ;
00182 exit(1) ;
00183 }
00184 if( nz_gp > nz_ds ){
00185 int npad1 = (nz_gp - nz_ds) / 2 ;
00186 int npad2 = (nz_gp - nz_ds) - npad1 ;
00187 int add_I=0, add_S=0, add_A=0, add_P=0, add_L=0, add_R=0 ;
00188 THD_3dim_dataset * pset ;
00189 char *sp1,*sp2 ;
00190
00191
00192
00193 switch( VL_dset->daxes->zzorient ){
00194 case ORI_R2L_TYPE:
00195 case ORI_L2R_TYPE: add_R=npad1; add_L=npad2; sp1="R"; sp2="L"; break;
00196
00197 case ORI_P2A_TYPE:
00198 case ORI_A2P_TYPE: add_A=npad1; add_P=npad2; sp1="A"; sp2="P"; break;
00199
00200 case ORI_I2S_TYPE:
00201 case ORI_S2I_TYPE: add_I=npad1; add_S=npad2; sp1="I"; sp2="S"; break;
00202 }
00203
00204
00205
00206 switch( ORIENT_sign[VL_dset->daxes->zzorient] ){
00207 default:
00208 case '+': npad_neg = npad1 ; npad_pos = npad2 ; break ;
00209 case '-': npad_neg = npad2 ; npad_pos = npad1 ; break ;
00210 }
00211 npadd = (npad_neg > 0 || npad_pos > 0 ) ;
00212
00213
00214
00215 if( VL_verbose )
00216 fprintf(stderr,"++ Zero padding to match -gridparent: -%s %d -%s %d\n",
00217 sp1,npad1,sp2,npad2 ) ;
00218
00219 pset = THD_zeropad( new_dset,
00220 add_I,add_S,add_A,add_P,add_L,add_R,
00221 NULL , ZPAD_EMPTY ) ;
00222
00223 if( pset == NULL ){
00224 fprintf(stderr,"** Fatal Error:\n"
00225 "** Can't properly zeropad output dataset!\n" ) ;
00226 exit(1) ;
00227 }
00228
00229
00230
00231 DSET_delete(new_dset); new_dset = pset;
00232 }
00233 }
00234
00235
00236
00237 EDIT_dset_items( new_dset , ADN_prefix , VL_prefix , ADN_none ) ;
00238 if( THD_is_file( DSET_HEADNAME(new_dset) ) ){
00239 fprintf(stderr,
00240 "** Output file %s already exists -- cannot continue!\n",
00241 DSET_HEADNAME(new_dset) ) ;
00242 exit(1) ;
00243 }
00244
00245 tross_Copy_History( VL_dset , new_dset ) ;
00246 tross_Make_History( "3dvolreg" , argc,argv , new_dset ) ;
00247
00248
00249
00250 if( VL_rotpar_dset != NULL ){
00251 ATR_float *atr ;
00252 float *matar , sum ;
00253 THD_fvec3 fv ;
00254 THD_dfvec3 dv,ev,qv , cv_e2, cv_e1, cv_s1, cv_s2 ;
00255
00256
00257
00258 atr = THD_find_float_atr( VL_rotpar_dset->dblk , "VOLREG_MATVEC_000000" );
00259 matar = atr->fl ;
00260 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
00261 matar[4],matar[5],matar[6],
00262 matar[8],matar[9],matar[10] ) ;
00263 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
00264
00265
00266
00267 pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ;
00268 sum = fabs(pp.mat[0][0]-1.)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0])
00269 +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.)+fabs(pp.mat[2][1])
00270 +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.) ;
00271 if( sum > 0.01 ) ERREX("-rotparent matrix not orthogonal!") ;
00272
00273
00274
00275
00276
00277
00278 fv = THD_dataset_center( new_dset ) ;
00279 FVEC3_TO_DFVEC3( fv , cv_e2 ) ;
00280
00281
00282
00283 if( VL_gridpar_dset != NULL ){
00284 fv = THD_dataset_center( VL_gridpar_dset ) ;
00285 FVEC3_TO_DFVEC3( fv , cv_e1 ) ;
00286 } else {
00287 cv_e1 = cv_e2 ;
00288 }
00289
00290
00291
00292 atr = THD_find_float_atr( VL_rotpar_dset->dblk , "VOLREG_CENTER_OLD" ) ;
00293 LOAD_DFVEC3( cv_s2 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
00294
00295
00296
00297 atr = THD_find_float_atr( VL_rotpar_dset->dblk , "VOLREG_CENTER_BASE" );
00298 LOAD_DFVEC3( cv_s1 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
00299
00300
00301
00302
00303
00304 dv = SUB_DFVEC3( cv_e2 , cv_s2 ) ;
00305 ev = DMATVEC( rmat , dv ) ;
00306
00307 dv = ev ;
00308
00309 ev = SUB_DFVEC3( cv_e1 , cv_s1 ) ;
00310
00311 qv = SUB_DFVEC3( dv , ev ) ;
00312
00313 tvec = ADD_DFVEC3( tvec , qv ) ;
00314
00315
00316
00317 pp = DBLE_mat_to_dicomm( new_dset ) ;
00318 ppt = TRANSPOSE_DMAT(pp);
00319 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp);
00320 tvec = DMATVEC(ppt,tvec);
00321
00322
00323
00324 if( VL_gridpar_dset != NULL ){
00325 new_dset->daxes->xxorg = VL_gridpar_dset->daxes->xxorg ;
00326 new_dset->daxes->yyorg = VL_gridpar_dset->daxes->yyorg ;
00327 new_dset->daxes->zzorg = VL_gridpar_dset->daxes->zzorg ;
00328
00329
00330
00331 if( new_dset->taxis != NULL && new_dset->taxis->nsl > 0 ){
00332 new_dset->taxis->zorg_sl = new_dset->daxes->zzorg ;
00333 }
00334 }
00335
00336 matvec = MATVEC_ORDER ;
00337 }
00338
00339
00340
00341 if( new_dset->taxis != NULL && new_dset->taxis->nsl > 0 && matvec ){
00342 int ndz ;
00343 int kk,jj , nsl = new_dset->taxis->nsl ;
00344
00345 ndz = (int)rint( tvec.xyz[2] / fabs(new_dset->daxes->zzdel) );
00346
00347 if( ndz != 0 ){
00348 float * tsl = (float *)malloc(sizeof(float)*nsl) ;
00349 for( kk=0 ; kk < nsl ; kk ++ ){
00350 jj = kk - ndz ;
00351 if( jj < 0 || jj >= nsl ) tsl[kk] = 0.0 ;
00352 else tsl[kk] = new_dset->taxis->toff_sl[jj] ;
00353 }
00354 EDIT_dset_items( new_dset , ADN_toff_sl , tsl , ADN_none ) ;
00355 free(tsl) ;
00356 if( VL_verbose )
00357 fprintf(stderr,"++ adjusting time-offsets by %d slices\n",ndz) ;
00358 }
00359 }
00360
00361
00362
00363 if( VL_verbose ){
00364 if( VL_tshift )
00365 fprintf(stderr,"++ Time shifting input dataset %s\n",
00366 DSET_BRIKNAME(VL_dset)) ;
00367 else
00368 fprintf(stderr,"++ Reading input dataset %s\n",
00369 DSET_BRIKNAME(VL_dset)) ;
00370 }
00371
00372 if( VL_tshift ){
00373 int eee = THD_dataset_tshift( VL_dset , VL_tshift_ignore ) ;
00374 if( eee )
00375 fprintf(stderr,"++ WARNING: some error during -tshift operation!\n") ;
00376 else
00377 EDIT_dset_items( new_dset ,
00378 ADN_nsl , 0 ,
00379 ADN_ttorg , 0.0 ,
00380 ADN_ttdur , 0.0 ,
00381 ADN_none ) ;
00382 }
00383 DSET_load( VL_dset ) ;
00384
00385
00386
00387 if( VL_imbase == NULL ){
00388 VL_imbase = mri_to_float(DSET_BRICK(VL_dset,VL_nbase)) ;
00389 } else {
00390 VL_nbase = -1 ;
00391 }
00392
00393 VL_imbase->dx = fabs( DSET_DX(VL_dset) ) ;
00394 VL_imbase->dy = fabs( DSET_DY(VL_dset) ) ;
00395 VL_imbase->dz = fabs( DSET_DZ(VL_dset) ) ;
00396 imb = MRI_FLOAT_PTR( VL_imbase ) ;
00397
00398 nx = DSET_NX(VL_dset) ; ny = DSET_NY(VL_dset) ; nz = DSET_NZ(VL_dset) ;
00399
00400
00401
00402 { int xf=0,yf=0,zf=0 ;
00403 switch( VL_edperc ){
00404 case 0:
00405 xf = (int)( MIN(0.25*nx,VL_edging) ) ;
00406 yf = (int)( MIN(0.25*ny,VL_edging) ) ;
00407 zf = (int)( MIN(0.25*nz,VL_edging) ) ;
00408 break ;
00409
00410 case 1:
00411 xf = (int)( 0.01*VL_edging*nx + 0.5 ) ;
00412 yf = (int)( 0.01*VL_edging*ny + 0.5 ) ;
00413 zf = (int)( 0.01*VL_edging*nz + 0.5 ) ;
00414 break ;
00415 }
00416 mri_3dalign_edging(xf,yf,zf) ;
00417 if( VL_verbose )
00418 fprintf(stderr,"++ Edging: x=%d y=%d z=%d\n",xf,yf,zf) ;
00419 }
00420
00421
00422
00423 if( VL_twopass ){
00424 MRI_IMAGE * tp_base ;
00425 int sx=66666,sy,sz ;
00426
00427 if( VL_verbose ){
00428 fprintf(stderr,"++ Start of first pass alignment on all sub-bricks\n");
00429 cputim = COX_cpu_time();
00430 }
00431
00432 tp_base = mri_to_float(VL_imbase) ;
00433
00434 EDIT_blur_volume_3d( nx,ny,nz , 1.0,1.0,1.0 ,
00435 MRI_float , MRI_FLOAT_PTR(tp_base) ,
00436 VL_twoblur,VL_twoblur,VL_twoblur ) ;
00437
00438 mri_3dalign_params( VL_maxite ,
00439 VL_twoblur*VL_dxy, VL_twoblur*VL_dph,
00440 VL_twoblur*VL_del,
00441 abs(ax1)-1 , abs(ax2)-1 , abs(ax3)-1 , -1 ) ;
00442
00443
00444
00445 mri_3dalign_method( MRI_LINEAR , (VL_verbose>1) , 1 , 0 ) ;
00446
00447
00448
00449
00450 if( VL_imwt != NULL || !VL_twosum || VL_imbase == DSET_BRICK(VL_dset,0) ){
00451
00452 albase = mri_3dalign_setup( tp_base , VL_imwt ) ;
00453
00454 } else {
00455
00456 float *far , *bar=MRI_FLOAT_PTR(tp_base) , *qar , clip ;
00457 int ii,jj,kk , nxy=nx*ny , nxyz=nxy*nz ;
00458 int nxbot,nxtop,nybot,nytop,nzbot,nztop , ee,fade,ff ;
00459
00460 if( VL_verbose )
00461 fprintf(stderr,
00462 "++ Computing first pass weight as sum of base and brick\n");
00463
00464 fim = mri_to_float( DSET_BRICK(VL_dset,0) ) ;
00465 far = MRI_FLOAT_PTR(fim) ;
00466 EDIT_blur_volume_3d( nx,ny,nz , 1.0,1.0,1.0 ,
00467 MRI_float , far ,
00468 VL_twoblur,VL_twoblur,VL_twoblur ) ;
00469
00470
00471
00472 if( VL_coarse_del > 0 && VL_coarse_num > 0 ){
00473 if( VL_verbose )
00474 fprintf(stderr,"++ Getting best coarse shift [0]:") ;
00475 get_best_shift( nx,ny,nz , bar,far , &sx,&sy,&sz ) ;
00476 if( VL_verbose ) fprintf(stderr," %d %d %d\n",sx,sy,sz) ;
00477 } else {
00478 sx = sy = sz = 0 ;
00479 }
00480
00481 #define BAR(i,j,k) bar[(i)+(j)*nx+(k)*nxy]
00482 #define QAR(i,j,k) qar[(i)+(j)*nx+(k)*nxy]
00483 #define FAR(i,j,k) far[(i)+(j)*nx+(k)*nxy]
00484
00485 qim = mri_copy(tp_base) ; qar = MRI_FLOAT_PTR(qim) ;
00486
00487 ee = abs(sx) ; nxbot = ee ; nxtop = nx-ee ;
00488 ee = abs(sy) ; nybot = ee ; nytop = ny-ee ;
00489 ee = abs(sz) ; nzbot = ee ; nztop = nz-ee ;
00490
00491 if( VL_sinit ){
00492 float sf=0.0,sq=0.0 ;
00493 for( kk=nzbot ; kk < nztop ; kk++ )
00494 for( jj=nybot ; jj < nytop ; jj++ )
00495 for( ii=nxbot ; ii < nxtop ; ii++ ){
00496 sf += FAR(ii-sx,jj-sy,kk-sz) ; sq += QAR(ii,jj,kk) ;
00497 }
00498 if( sq > 0.0 ){
00499 sf = sf / sq ;
00500 if( sf > 0.005 && sf < 2000.0 ){
00501 mri_3dalign_scaleinit(sf) ;
00502 if (sf < 200.0) {
00503 if (VL_verbose) fprintf(stderr,"++ Scale init = %g\n",sf) ;
00504 } else {
00505 fprintf(stderr,"++ Warning: Scale init = %g is large. Check output.\n",sf) ;
00506 }
00507 } else {
00508 fprintf(stderr,"-- Large scale difference between datasets.\n"
00509 " Scale init = %g\n"
00510 " 3dvolreg might not converge.",sf) ;
00511 }
00512 }
00513 }
00514
00515
00516
00517 for( kk=nzbot ; kk < nztop ; kk++ )
00518 for( jj=nybot ; jj < nytop ; jj++ )
00519 for( ii=nxbot ; ii < nxtop ; ii++ )
00520 QAR(ii,jj,kk) += FAR(ii-sx,jj-sy,kk-sz) ;
00521
00522 mri_free(fim) ;
00523
00524
00525
00526 if( VL_verbose )
00527 fprintf(stderr,"++ Blurring first pass weight\n") ;
00528
00529 #if 1
00530 EDIT_blur_volume_3d( nx,ny,nz , 1.0,1.0,1.0 ,
00531 MRI_float , qar ,
00532 VL_twoblur,VL_twoblur,VL_twoblur ) ;
00533 #else
00534 MRI_5blur_inplace_3D( qim ) ;
00535 #endif
00536
00537 clip = 0.025 * mri_max(qim) ;
00538 for( ii=0 ; ii < nxyz ; ii++ )
00539 if( qar[ii] < clip ) qar[ii] = 0.0 ;
00540
00541 mri_3dalign_force_edging( 1 ) ;
00542 albase = mri_3dalign_setup( tp_base , qim ) ;
00543 mri_3dalign_force_edging( 0 ) ;
00544 mri_free(qim) ;
00545 }
00546
00547
00548
00549 if( albase == NULL ){
00550 fprintf(stderr,
00551 "** Can't initialize first pass alignment algorithm\n");
00552 exit(1);
00553 }
00554
00555 dx_1 = (float *) malloc( sizeof(float) * imcount ) ;
00556 dy_1 = (float *) malloc( sizeof(float) * imcount ) ;
00557 dz_1 = (float *) malloc( sizeof(float) * imcount ) ;
00558 roll_1 = (float *) malloc( sizeof(float) * imcount ) ;
00559 pitch_1 = (float *) malloc( sizeof(float) * imcount ) ;
00560 yaw_1 = (float *) malloc( sizeof(float) * imcount ) ;
00561
00562
00563
00564
00565 for( kim=0 ; kim < imcount ; kim++ ){
00566
00567 qim = DSET_BRICK( VL_dset , kim ) ;
00568 fim = mri_to_float( qim ) ;
00569 fim->dx = fabs( DSET_DX(VL_dset) ) ;
00570 fim->dy = fabs( DSET_DY(VL_dset) ) ;
00571 fim->dz = fabs( DSET_DZ(VL_dset) ) ;
00572
00573 if( kim != VL_nbase ){
00574
00575 EDIT_blur_volume_3d( nx,ny,nz , 1.0,1.0,1.0 ,
00576 MRI_float , MRI_FLOAT_PTR(fim) ,
00577 VL_twoblur,VL_twoblur,VL_twoblur ) ;
00578
00579 if( kim > 0 || sx == 66666 ){
00580 if( VL_coarse_del > 0 && VL_coarse_num > 0 ){
00581 if( VL_verbose )
00582 fprintf(stderr,"++ Getting best coarse shift [%d]:",kim) ;
00583 get_best_shift( nx,ny,nz ,
00584 MRI_FLOAT_PTR(tp_base),MRI_FLOAT_PTR(fim) ,
00585 &sx,&sy,&sz ) ;
00586 if( VL_verbose )
00587 fprintf(stderr," %d %d %d\n",sx,sy,sz) ;
00588 } else {
00589 sx = sy = sz = 0 ;
00590 }
00591 }
00592
00593 mri_3dalign_initvals( 0.0 , 0.0 , 0.0 ,
00594 sx*fim->dx , sy*fim->dy , sz*fim->dz ) ;
00595
00596 (void) mri_3dalign_one( albase , fim ,
00597 roll_1+kim , pitch_1+kim , yaw_1+kim ,
00598 dx_1 +kim , dy_1 +kim , dz_1 +kim ) ;
00599
00600 roll_1[kim] *= (180.0/PI) ;
00601 pitch_1[kim] *= (180.0/PI) ;
00602 yaw_1[kim] *= (180.0/PI) ;
00603
00604 } else {
00605 roll_1[kim] =
00606 pitch_1[kim] =
00607 yaw_1[kim] =
00608 dx_1[kim] =
00609 dy_1[kim] =
00610 dz_1[kim] = 0.0 ;
00611 }
00612
00613 mri_free(fim) ;
00614 }
00615
00616 mri_3dalign_cleanup( albase ) ; mri_free( tp_base ) ;
00617
00618 if( VL_verbose ){
00619 cputim = COX_cpu_time() - cputim ;
00620 fprintf(stderr,"++ CPU time for first pass=%.3g s\n" , cputim) ;
00621 }
00622
00623 }
00624
00625
00626
00627
00628 if( VL_verbose ) fprintf(stderr,"++ Initializing alignment base\n") ;
00629
00630 mri_3dalign_params( VL_maxite , VL_dxy , VL_dph , VL_del ,
00631 abs(ax1)-1 , abs(ax2)-1 , abs(ax3)-1 , -1 ) ;
00632
00633
00634
00635
00636
00637 mri_3dalign_method( VL_resam , (VL_verbose>1) , (matvec != 0) , VL_clipit );
00638
00639 if( VL_final < 0 ) VL_final = VL_resam ;
00640 mri_3dalign_final_regmode( VL_final ) ;
00641
00642
00643
00644 if( VL_wtinp ){
00645 VL_imwt = mri_to_float( DSET_BRICK(VL_dset,0) ) ;
00646 mri_3dalign_wproccing( 1 ) ;
00647 }
00648
00649 albase = mri_3dalign_setup( VL_imbase , VL_imwt ) ;
00650 if( albase == NULL ){
00651 fprintf(stderr,"** Can't initialize base image for alignment\n"); exit(1);
00652 }
00653 if( VL_imwt != NULL ) mri_free( VL_imwt ) ;
00654
00655
00656
00657 if( VL_verbose ){
00658 fprintf(stderr,"++ Starting final pass on %d sub-bricks: ",imcount);
00659 cputim = COX_cpu_time();
00660 }
00661
00662 dxbar = dybar = dzbar = rollbar = yawbar = pitchbar = 0.0 ;
00663
00664 for( kim=0 ; kim < imcount ; kim++ ){
00665
00666 if( VL_verbose ) fprintf(stderr,"%d",kim) ;
00667
00668 qim = DSET_BRICK( VL_dset , kim ) ;
00669 fim = mri_to_float( qim ) ;
00670 fim->dx = fabs( DSET_DX(VL_dset) ) ;
00671 fim->dy = fabs( DSET_DY(VL_dset) ) ;
00672 fim->dz = fabs( DSET_DZ(VL_dset) ) ;
00673
00674
00675
00676 if( kim != VL_nbase ){
00677
00678 if( VL_twopass )
00679 mri_3dalign_initvals( roll_1[kim] , pitch_1[kim] , yaw_1[kim] ,
00680 dx_1[kim] , dy_1[kim] , dz_1[kim] ) ;
00681
00682 tim = mri_3dalign_one( albase , fim ,
00683 roll+kim , pitch+kim , yaw+kim ,
00684 &ddx , &ddy , &ddz ) ;
00685
00686 } else {
00687
00688 if( !matvec ) tim = mri_to_float( VL_imbase ) ;
00689 else tim = NULL ;
00690
00691 roll[kim] = pitch[kim] = yaw[kim] = ddx = ddy = ddz = 0.0 ;
00692
00693 }
00694
00695
00696
00697
00698
00699 if( VL_verbose ) fprintf(stderr,".");
00700
00701
00702
00703 roll[kim] *= (180.0/PI) ; if( hax1 < 0 ) roll[kim] = -roll[kim] ;
00704 pitch[kim] *= (180.0/PI) ; if( hax2 < 0 ) pitch[kim] = -pitch[kim] ;
00705 yaw[kim] *= (180.0/PI) ; if( hax3 < 0 ) yaw[kim] = -yaw[kim] ;
00706
00707 switch( new_dset->daxes->xxorient ){
00708 case ORI_R2L_TYPE: dy[kim] = ddx ; break ;
00709 case ORI_L2R_TYPE: dy[kim] = -ddx ; break ;
00710 case ORI_P2A_TYPE: dz[kim] = -ddx ; break ;
00711 case ORI_A2P_TYPE: dz[kim] = ddx ; break ;
00712 case ORI_I2S_TYPE: dx[kim] = ddx ; break ;
00713 case ORI_S2I_TYPE: dx[kim] = -ddx ; break ;
00714 }
00715
00716 switch( new_dset->daxes->yyorient ){
00717 case ORI_R2L_TYPE: dy[kim] = ddy ; break ;
00718 case ORI_L2R_TYPE: dy[kim] = -ddy ; break ;
00719 case ORI_P2A_TYPE: dz[kim] = -ddy ; break ;
00720 case ORI_A2P_TYPE: dz[kim] = ddy ; break ;
00721 case ORI_I2S_TYPE: dx[kim] = ddy ; break ;
00722 case ORI_S2I_TYPE: dx[kim] = -ddy ; break ;
00723 }
00724
00725 switch( new_dset->daxes->zzorient ){
00726 case ORI_R2L_TYPE: dy[kim] = ddz ; break ;
00727 case ORI_L2R_TYPE: dy[kim] = -ddz ; break ;
00728 case ORI_P2A_TYPE: dz[kim] = -ddz ; break ;
00729 case ORI_A2P_TYPE: dz[kim] = ddz ; break ;
00730 case ORI_I2S_TYPE: dx[kim] = ddz ; break ;
00731 case ORI_S2I_TYPE: dx[kim] = -ddz ; break ;
00732 }
00733
00734
00735
00736
00737 if( matvec ){
00738 THD_dvecmat vm1 , vm2 , vmtot ;
00739 char sbuf[128] ;
00740
00741 vm2.mm = rmat ; vm2.vv = tvec ;
00742
00743 sprintf(sbuf,"-rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP" ,
00744 roll[kim],pitch[kim],yaw[kim], dx[kim],dy[kim],dz[kim] ) ;
00745 vm1 = THD_rotcom_to_matvec( new_dset , sbuf ) ;
00746
00747 vmtot = MUL_DVECMAT(vm2,vm1 ) ;
00748
00749
00750
00751 if( npadd ){
00752 MRI_IMAGE * qim = mri_zeropad_3D( 0,0,0,0,npad_neg,npad_pos , fim ) ;
00753 if( qim == NULL ){
00754 fprintf(stderr,"** Can't zeropad at kim=%d -- FATAL ERROR!\n",kim);
00755 exit(1) ;
00756 }
00757 mri_free(fim) ; fim = qim ;
00758 }
00759
00760 THD_rota_method( VL_final ) ;
00761 tim = THD_rota3D_matvec( fim , vmtot.mm,vmtot.vv ) ;
00762
00763 if( VL_clipit &&
00764 (VL_final == MRI_QUINTIC || VL_final==MRI_CUBIC ||
00765 VL_final == MRI_HEPTIC || VL_final==MRI_FOURIER )){
00766
00767 register int ii ;
00768 register float ftop , fbot , * tar ;
00769
00770 ftop = mri_max( fim ); fbot = mri_min( fim );
00771 tar = MRI_FLOAT_PTR(tim) ;
00772 for( ii=0 ; ii < tim->nvox ; ii++ ){
00773 if( tar[ii] < fbot ) tar[ii] = fbot ;
00774 else if( tar[ii] > ftop ) tar[ii] = ftop ;
00775 }
00776 }
00777 }
00778
00779
00780
00781 if( kim == 0 ){
00782 dxtop = dxbot = dx[kim] ;
00783 dytop = dybot = dy[kim] ;
00784 dztop = dzbot = dz[kim] ;
00785 rolltop = rollbot = roll[kim] ;
00786 pitchtop = pitchbot = pitch[kim] ;
00787 yawtop = yawbot = yaw[kim] ;
00788 } else {
00789 dxtop = MAX(dxtop ,dx[kim] ); dxbot = MIN(dxbot ,dx[kim] );
00790 dytop = MAX(dytop ,dy[kim] ); dybot = MIN(dybot ,dy[kim] );
00791 dztop = MAX(dztop ,dz[kim] ); dzbot = MIN(dzbot ,dz[kim] );
00792 rolltop = MAX(rolltop ,roll[kim] ); rollbot = MIN(rollbot ,roll[kim] );
00793 pitchtop = MAX(pitchtop,pitch[kim]); pitchbot = MIN(pitchbot,pitch[kim]);
00794 yawtop = MAX(yawtop ,yaw[kim] ); yawbot = MIN(yawbot ,yaw[kim] );
00795 }
00796
00797 dxbar += dx[kim] ; dybar += dy[kim] ; dzbar += dz[kim] ;
00798 rollbar += roll[kim] ; pitchbar += pitch[kim] ; yawbar += yaw[kim] ;
00799
00800 if( !matvec ){
00801 sum = 0.0 ;
00802 tar = MRI_FLOAT_PTR(tim) ;
00803 for( ii=0 ; ii < tim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
00804 rmsnew[kim] = sqrt( sum / tim->nvox ) ;
00805
00806 sum = 0.0 ;
00807 tar = MRI_FLOAT_PTR(fim) ;
00808 for( ii=0 ; ii < fim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
00809 rmsold[kim] = sqrt( sum / fim->nvox ) ;
00810 } else {
00811 rmsold[kim] = rmsnew[kim] = 0.0 ;
00812 }
00813
00814 mri_free(fim) ;
00815
00816
00817
00818
00819
00820 switch( qim->kind ){
00821
00822 case MRI_float:
00823 EDIT_substitute_brick( new_dset, kim, MRI_float, MRI_FLOAT_PTR(tim) );
00824 mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
00825 break ;
00826
00827 case MRI_short:
00828 fim = mri_to_short(1.0,tim) ; mri_free( tim ) ;
00829 EDIT_substitute_brick( new_dset, kim, MRI_short, MRI_SHORT_PTR(fim) );
00830 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
00831 break ;
00832
00833 case MRI_byte:
00834 fim = mri_to_byte(tim) ; mri_free( tim ) ;
00835 EDIT_substitute_brick( new_dset, kim, MRI_byte, MRI_BYTE_PTR(fim) ) ;
00836 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
00837 break ;
00838
00839
00840
00841 default:
00842 fprintf(stderr,"\n*** Can't register bricks of type %s\n",
00843 MRI_TYPE_name[qim->kind] ) ;
00844 exit(1) ;
00845 }
00846
00847 DSET_unload_one( VL_dset , kim ) ;
00848
00849 if( VL_verbose ) fprintf(stderr,".") ;
00850
00851 }
00852
00853
00854
00855 mri_3dalign_cleanup( albase ) ;
00856 DSET_unload( VL_dset ) ;
00857 mri_free( VL_imbase ) ;
00858 if( VL_twopass ){
00859 free(dx_1);free(dy_1);free(dz_1);free(roll_1);free(pitch_1);free(yaw_1);
00860 }
00861
00862
00863
00864 dxbar /= imcount ; dybar /= imcount ; dzbar /= imcount ;
00865 rollbar /= imcount ; pitchbar /= imcount ; yawbar /= imcount ;
00866
00867 if( VL_verbose ){
00868 cputim = COX_cpu_time() - cputim ;
00869 fprintf(stderr,"\n++ CPU time for realignment=%.3g s" , cputim) ;
00870 if( imcount > 1 ) fprintf(stderr," [=%.3g s/sub-brick]" , cputim/imcount) ;
00871 fprintf(stderr,"\n") ;
00872
00873 fprintf(stderr,"++ Min : roll=%+.3f pitch=%+.3f yaw=%+.3f"
00874 " dS=%+.3f dL=%+.3f dP=%+.3f\n" ,
00875 rollbot , pitchbot , yawbot , dxbot , dybot , dzbot ) ;
00876
00877 fprintf(stderr,"++ Mean: roll=%+.3f pitch=%+.3f yaw=%+.3f"
00878 " dS=%+.3f dL=%+.3f dP=%+.3f\n" ,
00879 rollbar , pitchbar , yawbar , dxbar , dybar , dzbar ) ;
00880
00881 fprintf(stderr,"++ Max : roll=%+.3f pitch=%+.3f yaw=%+.3f"
00882 " dS=%+.3f dL=%+.3f dP=%+.3f\n" ,
00883 rolltop , pitchtop , yawtop , dxtop , dytop , dztop ) ;
00884 }
00885
00886
00887
00888 if( imcount == 1 && VL_rotpar_dset == NULL ){
00889 char *str = NULL ;
00890 str = THD_zzprintf( str , "3dvolreg did: %s" , modes[VL_final] ) ;
00891 if( VL_clipit ) str = THD_zzprintf( str , " -clipit" ) ;
00892 else str = THD_zzprintf( str , " -noclip" ) ;
00893 if( VL_zpad ) str = THD_zzprintf( str , " -zpad %d" , VL_zpad ) ;
00894 str = THD_zzprintf(str,
00895 " -rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP\n",
00896 roll[0],pitch[0],yaw[0], dx[0],dy[0],dz[0] ) ;
00897 tross_Append_History( new_dset , str ) ;
00898 free(str) ;
00899 }
00900
00901
00902
00903
00904 { char sbuf[128] , anam[32] ;
00905 THD_fvec3 cv ;
00906 THD_dmat33 rmat ;
00907 float matar[12] ;
00908
00909
00910
00911 if( VL_rotpar_dset != NULL ){
00912 THD_set_string_atr(new_dset->dblk,"VOLREG_ROTPARENT_IDCODE",
00913 VL_rotpar_dset->idcode.str );
00914 THD_set_string_atr(new_dset->dblk,"VOLREG_ROTPARENT_NAME" ,
00915 DSET_HEADNAME(VL_rotpar_dset));
00916 }
00917
00918 if( VL_gridpar_dset != NULL ){
00919 THD_set_string_atr(new_dset->dblk,"VOLREG_GRIDPARENT_IDCODE",
00920 VL_gridpar_dset->idcode.str );
00921 THD_set_string_atr(new_dset->dblk,"VOLREG_GRIDPARENT_NAME" ,
00922 DSET_HEADNAME(VL_gridpar_dset));
00923 }
00924
00925
00926
00927 THD_set_string_atr( new_dset->dblk, "VOLREG_INPUT_IDCODE",
00928 VL_dset->idcode.str ) ;
00929 THD_set_string_atr( new_dset->dblk, "VOLREG_INPUT_NAME" ,
00930 DSET_HEADNAME(VL_dset) ) ;
00931
00932 cv = THD_dataset_center( new_dset ) ;
00933 THD_set_float_atr( new_dset->dblk , "VOLREG_CENTER_OLD" , 3 , cv.xyz ) ;
00934
00935
00936
00937 if( VL_bset == NULL ) VL_bset = VL_dset ;
00938
00939 THD_set_string_atr( new_dset->dblk , "VOLREG_BASE_IDCODE" ,
00940 VL_bset->idcode.str ) ;
00941 THD_set_string_atr( new_dset->dblk , "VOLREG_BASE_NAME" ,
00942 DSET_HEADNAME(VL_bset) ) ;
00943
00944 cv = THD_dataset_center( VL_bset ) ;
00945 THD_set_float_atr( new_dset->dblk , "VOLREG_CENTER_BASE" , 3 , cv.xyz ) ;
00946
00947
00948
00949 THD_set_int_atr( new_dset->dblk , "VOLREG_ROTCOM_NUM" , 1 , &imcount ) ;
00950
00951
00952
00953 for( kim=0 ; kim < imcount ; kim++ ){
00954 sprintf(anam,"VOLREG_ROTCOM_%06d",kim) ;
00955 sprintf(sbuf,"-rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP" ,
00956 roll[kim],pitch[kim],yaw[kim], dx[kim],dy[kim],dz[kim] ) ;
00957 THD_set_string_atr( new_dset->dblk , anam , sbuf ) ;
00958
00959
00960
00961
00962 rmat = rot_to_matrix( 2 , -(PI/180.0)*roll[kim] ,
00963 0 , -(PI/180.0)*pitch[kim] ,
00964 1 , -(PI/180.0)*yaw[kim] ) ;
00965
00966
00967
00968
00969
00970
00971
00972 UNLOAD_DMAT(rmat,matar[0],matar[1],matar[2],
00973 matar[4],matar[5],matar[6],
00974 matar[8],matar[9],matar[10] ) ;
00975 matar[3] = dy[kim] ; matar[7] = dz[kim] ; matar[11] = dx[kim] ;
00976 sprintf(anam,"VOLREG_MATVEC_%06d",kim) ;
00977 THD_set_float_atr( new_dset->dblk , anam , 12 , matar ) ;
00978 }
00979 }
00980
00981
00982
00983 if( VL_twodup && !VL_intern && VL_rotpar_dset == NULL ){
00984 new_dset->daxes->xxorg = VL_bxorg ;
00985 new_dset->daxes->yyorg = VL_byorg ;
00986 new_dset->daxes->zzorg = VL_bzorg ;
00987
00988 if( new_dset->taxis != NULL && new_dset->taxis->nsl > 0 ){
00989 new_dset->taxis->zorg_sl = new_dset->daxes->zzorg ;
00990 }
00991 }
00992
00993
00994
00995 DSET_write(new_dset) ;
00996 if( VL_verbose )
00997 fprintf(stderr,"++ Wrote dataset to disk in %s",DSET_BRIKNAME(new_dset));
00998 if( VL_verbose ) fprintf(stderr,"\n") ;
00999
01000
01001
01002 if( VL_dfile[0] != '\0' ){
01003 FILE *fp ;
01004
01005 if( THD_is_file(VL_dfile) )
01006 fprintf(stderr,"** Warning: overwriting file %s\n",VL_dfile) ;
01007
01008 fp = fopen( VL_dfile , "w" ) ;
01009 for( kim=0 ; kim < imcount ; kim++ )
01010 fprintf(fp , "%4d %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %11.4g %11.4g\n" ,
01011 kim , roll[kim], pitch[kim], yaw[kim],
01012 dx[kim], dy[kim], dz[kim],
01013 rmsold[kim] , rmsnew[kim] ) ;
01014 fclose(fp) ;
01015 }
01016
01017 if( VL_1Dfile[0] != '\0' ){
01018 FILE *fp ;
01019
01020 if( THD_is_file(VL_1Dfile) )
01021 fprintf(stderr,"** Warning: overwriting file %s\n",VL_1Dfile) ;
01022
01023 fp = fopen( VL_1Dfile , "w" ) ;
01024 for( kim=0 ; kim < imcount ; kim++ )
01025 fprintf(fp , "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n" ,
01026 roll[kim], pitch[kim], yaw[kim],
01027 dx[kim] , dy[kim] , dz[kim] ) ;
01028 fclose(fp) ;
01029 }
01030
01031 if( VL_rotcom ){
01032
01033 printf("\n3drotate fragment%s:\n\n", (imcount > 1)? "s" : "" ) ;
01034 for( kim=0 ; kim < imcount ; kim++ ){
01035 printf("3drotate %s" , modes[VL_final] ) ;
01036 if( VL_clipit ) printf(" -clipit" ) ;
01037 else printf(" -noclip" ) ;
01038 if( VL_zpad ) printf(" -zpad %d" , VL_zpad ) ;
01039 printf(" -rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP\n" ,
01040 roll[kim],pitch[kim],yaw[kim], dx[kim],dy[kim],dz[kim] ) ;
01041 }
01042 printf("\n") ;
01043 }
01044
01045 exit(0) ;
01046 }
01047
01048
01049
01050 void VL_syntax(void)
01051 {
01052 printf(
01053 "Usage: 3dvolreg [options] dataset\n"
01054 "Registers each 3D sub-brick from the input dataset to the base brick.\n"
01055 "'dataset' may contain a sub-brick selector list.\n"
01056 "\n"
01057 "OPTIONS:\n"
01058 " -verbose Print progress reports. Use twice for LOTS of output.\n"
01059 " -Fourier Perform the alignments using Fourier interpolation.\n"
01060 " -heptic Use heptic polynomial interpolation.\n"
01061 " -quintic Use quintic polynomical interpolation.\n"
01062 " -cubic Use cubic polynomial interpolation.\n"
01063 " Default = Fourier [slowest and most accurate interpolator]\n"
01064 " -clipit Clips the values in each output sub-brick to be in the same\n"
01065 " range as the corresponding input volume.\n"
01066 " The interpolation schemes can produce values outside\n"
01067 " the input range, which is sometimes annoying.\n"
01068 " [16 Apr 2002: -clipit is now the default]\n"
01069 " -noclip Turns off -clipit\n"
01070 " -zpad n Zeropad around the edges by 'n' voxels during rotations\n"
01071 " (these edge values will be stripped off in the output)\n"
01072 " N.B.: Unlike to3d, in this program '-zpad' adds zeros in\n"
01073 " all directions.\n"
01074 " N.B.: The environment variable AFNI_ROTA_ZPAD can be used\n"
01075 " to set a nonzero default value for this parameter.\n"
01076 " -prefix fname Use 'fname' for the output dataset prefix.\n"
01077 " The program tries not to overwrite an existing dataset.\n"
01078 " Default = 'volreg'.\n"
01079 "\n"
01080 " -base n Sets the base brick to be the 'n'th sub-brick\n"
01081 " from the input dataset (indexing starts at 0).\n"
01082 " Default = 0 (first sub-brick).\n"
01083 " -base 'bset[n]' Sets the base brick to be the 'n'th sub-brick\n"
01084 " from the dataset specified by 'bset', as in\n"
01085 " -base 'elvis+orig[4]'\n"
01086 " The quotes are needed because the '[]' characters\n"
01087 " are special to the shell.\n"
01088 "\n"
01089 " -dfile dname Save the motion parameters in file 'dname'.\n"
01090 " The output is in 9 ASCII formatted columns:\n"
01091 "\n"
01092 " n roll pitch yaw dS dL dP rmsold rmsnew\n"
01093 "\n"
01094 " where: n = sub-brick index\n"
01095 " roll = rotation about the I-S axis }\n"
01096 " pitch = rotation about the R-L axis } degrees CCW\n"
01097 " yaw = rotation about the A-P axis }\n"
01098 " dS = displacement in the Superior direction }\n"
01099 " dL = displacement in the Left direction } mm\n"
01100 " dP = displacement in the Posterior direction }\n"
01101 " rmsold = RMS difference between input brick and base brick\n"
01102 " rmsnew = RMS difference between output brick and base brick\n"
01103 " N.B.: If the '-dfile' option is not given, the parameters aren't saved.\n"
01104 " N.B.: The motion parameters are those needed to bring the sub-brick\n"
01105 " back into alignment with the base. In 3drotate, it is as if\n"
01106 " the following options were applied to each input sub-brick:\n"
01107 " -rotate <roll>I <pitch>R <yaw>A -ashift <dS>S <dL>L <dP>P\n"
01108 "\n"
01109 " -1Dfile ename Save the motion parameters ONLY in file 'ename'.\n"
01110 " The output is in 6 ASCII formatted columns:\n"
01111 "\n"
01112 " roll pitch yaw dS dL dP\n"
01113 "\n"
01114 " This file can be used in FIM as an 'ort', to detrend\n"
01115 " the data against correlation with the movements.\n"
01116 " This type of analysis can be useful in removing\n"
01117 " errors made in the interpolation.\n"
01118 "\n"
01119 " -rotcom Write the fragmentary 3drotate commands needed to\n"
01120 " perform the realignments to stdout; for example:\n"
01121 " 3drotate -rotate 7.2I 3.2R -5.7A -ashift 2.7S -3.8L 4.9P\n"
01122 " The purpose of this is to make it easier to shift other\n"
01123 " datasets using exactly the same parameters.\n"
01124 "\n"
01125 " -tshift ii If the input dataset is 3D+time and has slice-dependent\n"
01126 " time-offsets (cf. the output of 3dinfo -v), then this\n"
01127 " option tells 3dvolreg to time shift it to the average\n"
01128 " slice time-offset prior to doing the spatial registration.\n"
01129 " The integer 'ii' is the number of time points at the\n"
01130 " beginning to ignore in the time shifting. The results\n"
01131 " should like running program 3dTshift first, then running\n"
01132 " 3dvolreg -- this is primarily a convenience option.\n"
01133 " N.B.: If the base brick is taken from this dataset, as in\n"
01134 " '-base 4', then it will be the time shifted brick.\n"
01135 " If for some bizarre reason this is undesirable, you\n"
01136 " could use '-base this+orig[4]' instead.\n"
01137 "\n"
01138 " -rotparent rset\n"
01139 " Specifies that AFTER the registration algorithm finds the best\n"
01140 " transformation for each sub-brick of the input, an additional\n"
01141 " rotation+translation should be performed before computing the\n"
01142 " final output dataset; this extra transformation is taken from\n"
01143 " the first 3dvolreg transformation found in dataset 'rset'.\n"
01144 " -gridparent gset\n"
01145 " Specifies that the output dataset of 3dvolreg should be shifted to\n"
01146 " match the grid of dataset 'gset'. Can only be used with -rotparent.\n"
01147 " This dataset should be one this is properly aligned with 'rset' when\n"
01148 " overlaid in AFNI.\n"
01149 " * If 'gset' has a different number of slices than the input dataset,\n"
01150 " then the output dataset will be zero-padded in the slice direction\n"
01151 " to match 'gset'.\n"
01152 " * These options are intended to be used to align datasets between sessions:\n"
01153 " S1 = SPGR from session 1 E1 = EPI from session 1\n"
01154 " S2 = SPGR from session 2 E2 = EPI from session 2\n"
01155 " 3dvolreg -twopass -twodup -base S1+orig -prefix S2reg S2+orig\n"
01156 " 3dvolreg -rotparent S2reg+orig -gridparent E1+orig -prefix E2reg \\\n"
01157 " -base 4 E2+orig\n"
01158 " Each sub-brick in E2 is registered to sub-brick E2+orig[4], then the\n"
01159 " rotation from S2 to S2reg is also applied, which shifting+padding\n"
01160 " applied to properly overlap with E1.\n"
01161 " * A similar effect could be done by using commands\n"
01162 " 3dvolreg -twopass -twodup -base S1+orig -prefix S2reg S2+orig\n"
01163 " 3dvolreg -prefix E2tmp -base 4 E2+orig\n"
01164 " 3drotate -rotparent S2reg+orig -gridparent E1+orig -prefix E2reg E2tmp+orig\n"
01165 " The principal difference is that the latter method results in E2\n"
01166 " being interpolated twice to make E2reg: once in the 3dvolreg run to\n"
01167 " produce E2tmp, then again when E2tmp is rotated to make E2reg. Using\n"
01168 " 3dvolreg with the -rotparent and -gridparent options simply skips the\n"
01169 " intermediate interpolation.\n"
01170 "\n"
01171 " *** Please read file README.registration for more ***\n"
01172 " *** information on the use of 3dvolreg and 3drotate ***\n"
01173 "\n"
01174 " Algorithm: Iterated linearized weighted least squares to make each\n"
01175 " sub-brick as like as possible to the base brick.\n"
01176 " This method is useful for finding SMALL MOTIONS ONLY.\n"
01177 " See program 3drotate for the volume shift/rotate algorithm.\n"
01178 " The following options can be used to control the iterations:\n"
01179 " -maxite m = Allow up to 'm' iterations for convergence\n"
01180 " [default = %d].\n"
01181 " -x_thresh x = Iterations converge when maximum movement\n"
01182 " is less than 'x' voxels [default=%f],\n"
01183 " -rot_thresh r = And when maximum rotation is less than\n"
01184 " 'r' degrees [default=%f].\n"
01185 " -delta d = Distance, in voxel size, used to compute\n"
01186 " image derivatives using finite differences\n"
01187 " [default=%f].\n"
01188 " -final mode = Do the final interpolation using the method\n"
01189 " defined by 'mode', which is one of the\n"
01190 " strings 'NN', 'cubic', 'quintic', 'heptic',\n"
01191 " or 'Fourier'\n"
01192 " [default=mode used to estimate parameters].\n"
01193 " -weight 'wset[n]' = Set the weighting applied to each voxel\n"
01194 " proportional to the brick specified here\n"
01195 " [default=smoothed base brick].\n"
01196 " N.B.: if no weight is given, and -twopass is\n"
01197 " engaged, then the first pass weight is the\n"
01198 " blurred sum of the base brick and the first\n"
01199 " data brick to be registered.\n"
01200 " -edging ee = Set the size of the region around the edges of\n"
01201 " the base volume where the default weight will\n"
01202 " be set to zero. If 'ee' is a plain number,\n"
01203 " then it is a voxel count, giving the thickness\n"
01204 " along each face of the 3D brick. If 'ee' is\n"
01205 " of the form '5%%', then it is a fraction of\n"
01206 " of each brick size. For example, '5%%' of\n"
01207 " a 256x256x124 volume means that 13 voxels\n"
01208 " on each side of the xy-axes will get zero\n"
01209 " weight, and 6 along the z-axis. If this\n"
01210 " option is not used, then 'ee' is read from\n"
01211 " the environment variable AFNI_VOLREG_EDGING.\n"
01212 " If that variable is not set, then 5%% is used.\n"
01213 " N.B.: This option has NO effect if the -weight\n"
01214 " option is used.\n"
01215 " N.B.: The largest %% value allowed is 25%%.\n"
01216 " -twopass = Do two passes of the registration algorithm:\n"
01217 " (1) with smoothed base and data bricks, with\n"
01218 " linear interpolation, to get a crude\n"
01219 " alignment, then\n"
01220 " (2) with the input base and data bricks, to\n"
01221 " get a fine alignment.\n"
01222 " This method is useful when aligning high-\n"
01223 " resolution datasets that may need to be\n"
01224 " moved more than a few voxels to be aligned.\n"
01225 " -twoblur bb = 'bb' is the blurring factor for pass 1 of\n"
01226 " the -twopass registration. This should be\n"
01227 " a number >= 2.0 (which is the default).\n"
01228 " Larger values would be reasonable if pass 1\n"
01229 " has to move the input dataset a long ways.\n"
01230 " Use '-verbose -verbose' to check on the\n"
01231 " iterative progress of the passes.\n"
01232 " N.B.: when using -twopass, and you expect the\n"
01233 " data bricks to move a long ways, you might\n"
01234 " want to use '-heptic' rather than\n"
01235 " the default '-Fourier', since you can get\n"
01236 " wraparound from Fourier interpolation.\n"
01237 " -twodup = If this option is set, along with -twopass,\n"
01238 " then the output dataset will have its\n"
01239 " xyz-axes origins reset to those of the\n"
01240 " base dataset. This is equivalent to using\n"
01241 " '3drefit -duporigin' on the output dataset.\n"
01242 " -sinit = When using -twopass registration on volumes\n"
01243 " whose magnitude differs significantly, the\n"
01244 " least squares fitting procedure is started\n"
01245 " by doing a zero-th pass estimate of the\n"
01246 " scale difference between the bricks.\n"
01247 " Use this option to turn this feature OFF.\n"
01248 " -coarse del num = When doing the first pass, the first step is\n"
01249 " to do a number of coarse shifts in order to\n"
01250 " find a starting point for the iterations.\n"
01251 " 'del' is the size of these steps, in voxels;\n"
01252 " 'num' is the number of these steps along\n"
01253 " each direction (+x,-x,+y,-y,+z,-z). The\n"
01254 " default values are del=10 and num=2. If\n"
01255 " you don't want this step performed, set\n"
01256 " num=0. Note that the amount of computation\n"
01257 " grows as num**3, so don't increase num\n"
01258 " past 4, or the program will run forever!\n"
01259 " N.B.: The 'del' parameter cannot be larger than\n"
01260 " 10%% of the smallest dimension of the input\n"
01261 " dataset.\n"
01262 #if 0
01263 " -wtrim = Attempt to trim the intermediate volumes to\n"
01264 " a smaller region (determined by the weight\n"
01265 " volume). The purpose of this is to save\n"
01266 " memory. The use of '-verbose -verbose'\n"
01267 " will report on the trimming usage.\n"
01268 " N.B.: At some point in the future, -wtrim will\n"
01269 " become the default.\n"
01270 #endif
01271 " -wtinp = Use sub-brick[0] of the input dataset as the\n"
01272 " weight brick in the final registration pass.\n"
01273 "\n"
01274 " N.B.: * This program can consume VERY large quantities of memory.\n"
01275 " (Rule of thumb: 40 bytes per input voxel.)\n"
01276 " Use of '-verbose -verbose' will show the amount of workspace,\n"
01277 " and the steps used in each iteration.\n"
01278 " * ALWAYS check the results visually to make sure that the program\n"
01279 " wasn't trapped in a 'false optimum'.\n"
01280 " * The default rotation threshold is reasonable for 64x64 images.\n"
01281 " You may want to decrease it proportionally for larger datasets.\n"
01282 " * -twopass resets the -maxite parameter to 66; if you want to use\n"
01283 " a different value, use -maxite AFTER the -twopass option.\n"
01284 " * The -twopass option can be slow; several CPU minutes for a\n"
01285 " 256x256x124 volume is a typical run time.\n"
01286 " * After registering high-resolution anatomicals, you may need to\n"
01287 " set their origins in 3D space to match. This can be done using\n"
01288 " the '-duporigin' option to program 3drefit, or by using the\n"
01289 " '-twodup' option to this program.\n"
01290
01291 , VL_maxite , VL_dxy , VL_dph , VL_del ) ;
01292
01293 return ;
01294 }
01295
01296
01297
01298 void VL_command_line(void)
01299 {
01300 int ii , nxbase , nybase , nerr , basecode ;
01301 MRI_IMAGE * tim ;
01302 MRI_IMARR * tarr ;
01303 float bdx,bdy,bdz ;
01304
01305
01306
01307 while( Iarg < Argc && Argv[Iarg][0] == '-' ){
01308
01309
01310
01311 if( strcmp(Argv[Iarg],"-sinit") == 0 ){
01312 VL_sinit = 0 ; Iarg++ ; continue ;
01313 }
01314
01315
01316
01317 if( strncmp(Argv[Iarg],"-params",4) == 0 ){
01318 VL_maxite = strtol( Argv[++Iarg] , NULL , 10 ) ;
01319 VL_dxy = strtod( Argv[++Iarg] , NULL ) ;
01320 VL_dph = strtod( Argv[++Iarg] , NULL ) ;
01321 VL_del = strtod( Argv[++Iarg] , NULL ) ;
01322 Iarg++ ; continue ;
01323 }
01324
01325 if( strncmp(Argv[Iarg],"-maxite",4) == 0 ){
01326 int mit = strtol( Argv[++Iarg] , NULL , 10 ) ;
01327 if( mit > 0 ) VL_maxite = mit ;
01328 Iarg++ ; continue ;
01329 }
01330
01331 if( strncmp(Argv[Iarg],"-x_thresh",4) == 0 ){
01332 float dxy = strtod( Argv[++Iarg] , NULL ) ;
01333 if( dxy > 0.0 ) VL_dxy = dxy ;
01334 Iarg++ ; continue ;
01335 }
01336
01337 if( strncmp(Argv[Iarg],"-rot_thresh",6) == 0 ){
01338 float dph = strtod( Argv[++Iarg] , NULL ) ;
01339 if( dph > 0.0 ) VL_dph = dph ;
01340 Iarg++ ; continue ;
01341 }
01342
01343 if( strncmp(Argv[Iarg],"-delta",4) == 0 ){
01344 float del = strtod( Argv[++Iarg] , NULL ) ;
01345 if( del > 0.0 ) VL_del = del ;
01346 Iarg++ ; continue ;
01347 }
01348
01349 if( strncmp(Argv[Iarg],"-final",4) == 0 ){
01350 char * str = Argv[++Iarg] ;
01351
01352 if( strcmp(str,"cubic") == 0 ) VL_final = MRI_CUBIC ;
01353 else if( strcmp(str,"quintic") == 0 ) VL_final = MRI_QUINTIC ;
01354 else if( strcmp(str,"heptic") == 0 ) VL_final = MRI_HEPTIC ;
01355 else if( strcmp(str,"Fourier") == 0 ) VL_final = MRI_FOURIER ;
01356 else if( strcmp(str,"NN") == 0 ) VL_final = MRI_NN ;
01357 else {
01358 fprintf(stderr,"** Illegal mode after -final\n"); exit(1);
01359 }
01360 Iarg++ ; continue ;
01361 }
01362
01363 if( strcmp(Argv[Iarg],"-rotcom") == 0 ){
01364 VL_rotcom++ ;
01365 Iarg++ ; continue ;
01366 }
01367
01368 if( strcmp(Argv[Iarg],"-twopass") == 0 ){
01369 VL_twopass++ ;
01370 if( VL_maxite < 10 ) VL_maxite = 66 ;
01371 Iarg++ ; continue ;
01372 }
01373
01374 if( strcmp(Argv[Iarg],"-twoblur") == 0 ){
01375 VL_twoblur = strtod( Argv[++Iarg] , NULL ) ;
01376 if( VL_twoblur < 2.0 ){
01377 fprintf(stderr,"** ERROR: value after -twoblur is < 2.0\n") ;
01378 exit(1) ;
01379 }
01380 Iarg++ ; continue ;
01381 }
01382
01383
01384
01385 if( strncmp(Argv[Iarg],"-verbose",5) == 0 ){
01386 VL_verbose++ ;
01387 Iarg++ ; continue ;
01388 }
01389
01390
01391
01392 if( strncmp(Argv[Iarg],"-clipit",4) == 0 ){
01393 fprintf(stderr,"++ Notice: -clipit is now the default\n") ;
01394 VL_clipit = 1 ;
01395 Iarg++ ; continue ;
01396 }
01397
01398 if( strncmp(Argv[Iarg],"-noclip",4) == 0 ){
01399 VL_clipit = 0 ;
01400 Iarg++ ; continue ;
01401 }
01402
01403
01404
01405 if( strncmp(Argv[Iarg],"-zpad",5) == 0 ){
01406 if( VL_zpad > 0 )
01407 fprintf(stderr,"++ WARNING: second -zpad option!\n") ;
01408 VL_zpad = (int) strtod( Argv[++Iarg] , NULL ) ;
01409 if( VL_zpad < 0 ){
01410 fprintf(stderr,"** ERROR: Can't use -zpad %d\n",VL_zpad) ;
01411 exit(1) ;
01412 }
01413 THD_rota_setpad(VL_zpad,VL_zpad,VL_zpad) ;
01414 Iarg++ ; continue ;
01415 }
01416
01417
01418
01419 if( strncmp(Argv[Iarg],"-Fourier",4) == 0 ||
01420 strncmp(Argv[Iarg],"-fourier",4) == 0 ){
01421
01422 VL_resam = MRI_FOURIER ;
01423 Iarg++ ; continue ;
01424 }
01425
01426
01427
01428 if( strncmp(Argv[Iarg],"-linear",4) == 0 ||
01429 strncmp(Argv[Iarg],"-Linear",4) == 0 ){
01430
01431 VL_resam = MRI_LINEAR ;
01432 Iarg++ ; continue ;
01433 }
01434
01435
01436
01437 if( strncmp(Argv[Iarg],"-cubic",4) == 0 ||
01438 strncmp(Argv[Iarg],"-Cubic",4) == 0 ){
01439
01440 VL_resam = MRI_CUBIC ;
01441 Iarg++ ; continue ;
01442 }
01443
01444
01445
01446 if( strncmp(Argv[Iarg],"-quintic",4) == 0 ||
01447 strncmp(Argv[Iarg],"-Quintic",4) == 0 ){
01448
01449 VL_resam = MRI_QUINTIC ;
01450 Iarg++ ; continue ;
01451 }
01452
01453
01454
01455 if( strncmp(Argv[Iarg],"-heptic",4) == 0 ||
01456 strncmp(Argv[Iarg],"-Heptic",4) == 0 ){
01457
01458 VL_resam = MRI_HEPTIC ;
01459 Iarg++ ; continue ;
01460 }
01461
01462
01463
01464
01465 if( strncmp(Argv[Iarg],"-prefix",4) == 0 ){
01466 strcpy( VL_prefix , Argv[++Iarg] ) ;
01467 Iarg++ ; continue ;
01468 }
01469
01470
01471
01472 if( strncmp(Argv[Iarg],"-dfile",4) == 0 ){
01473 strcpy( VL_dfile , Argv[++Iarg] ) ;
01474 Iarg++ ; continue ;
01475 }
01476
01477
01478
01479 if( strncmp(Argv[Iarg],"-1Dfile",4) == 0 ){
01480 strcpy( VL_1Dfile , Argv[++Iarg] ) ;
01481 Iarg++ ; continue ;
01482 }
01483
01484
01485
01486 if( strncmp(Argv[Iarg],"-base",4) == 0 ){
01487 int bb,ii ; char * cpt ;
01488
01489 if( VL_imbase != NULL || VL_nbase > 0 ){
01490 fprintf(stderr,"** Can't have two -base arguments\n") ; exit(1) ;
01491 }
01492
01493
01494
01495 bb = strtol( Argv[++Iarg] , &cpt , 10 ) ;
01496 if( bb < 0 ){
01497 fprintf(stderr,"** Illegal number after -base\n"); exit(1);
01498 }
01499
01500 if( *cpt == '\0' ){
01501
01502 VL_nbase = bb ;
01503 VL_imbase = NULL ;
01504 VL_intern = 1 ;
01505
01506 } else {
01507
01508
01509
01510
01511 VL_bset = THD_open_dataset( Argv[Iarg] ) ;
01512 if( VL_bset == NULL ){
01513 fprintf(stderr,"** Couldn't open -base dataset %s\n",Argv[Iarg]) ;
01514 exit(1) ;
01515 }
01516 if( VL_verbose )
01517 fprintf(stderr,
01518 "++ Reading in base dataset %s\n",DSET_BRIKNAME(VL_bset)) ;
01519 DSET_load(VL_bset) ;
01520 if( !DSET_LOADED(VL_bset) ){
01521 fprintf(stderr,"** Couldn't read -base dataset %s\n",
01522 DSET_BRIKNAME(VL_bset)) ;
01523 exit(1) ;
01524 }
01525 if( DSET_NVALS(VL_bset) > 1 )
01526 fprintf(stderr,
01527 "++ WARNING: -base dataset %s has more than 1 sub-brick\n",
01528 Argv[Iarg]) ;
01529
01530 VL_intern = 0 ;
01531
01532 bdx = fabs(DSET_DX(VL_bset)) ;
01533 bdy = fabs(DSET_DY(VL_bset)) ;
01534 bdz = fabs(DSET_DZ(VL_bset)) ;
01535
01536
01537
01538
01539
01540
01541
01542 VL_imbase = mri_to_float( DSET_BRICK(VL_bset,0) ) ;
01543
01544 VL_bxorg = VL_bset->daxes->xxorg ;
01545 VL_byorg = VL_bset->daxes->yyorg ;
01546 VL_bzorg = VL_bset->daxes->zzorg ;
01547
01548 DSET_unload( VL_bset ) ;
01549 }
01550 Iarg++ ; continue ;
01551 }
01552
01553
01554
01555 if( strcmp(Argv[Iarg],"-coarse") == 0 ){
01556 VL_coarse_del = strtol(Argv[++Iarg],NULL,10) ;
01557 VL_coarse_num = strtol(Argv[++Iarg],NULL,10) ;
01558 Iarg++ ; continue ;
01559 }
01560
01561
01562
01563 if( strcmp(Argv[Iarg],"-wtrim") == 0 ){
01564 #if 0
01565 mri_3dalign_wtrimming(1) ;
01566 #else
01567 fprintf(stderr,"++ Notice: -wtrim is now always enabled\n");
01568 #endif
01569 Iarg++ ; continue ;
01570 }
01571
01572
01573
01574 if( strcmp(Argv[Iarg],"-wtinp") == 0 ){
01575 VL_wtinp = 1 ;
01576 Iarg++ ; continue ;
01577 }
01578
01579
01580
01581 if( strcmp(Argv[Iarg],"-edging") == 0 ){
01582 float ee ; char * eq ;
01583 ee = strtod(Argv[++Iarg],&eq) ;
01584 if( ee < 0 ){
01585 fprintf(stderr,"** Illegal value after -edging\n"); exit(1);
01586 }
01587 if( *eq == '%' ){
01588 if( ee > 25.0 ){
01589 fprintf(stderr,"** Illegal percentage after -edging\n"); exit(1);
01590 }
01591 VL_edperc = 1 ; VL_edging = ee ;
01592 } else {
01593 VL_edperc = 0 ; VL_edging = ee ;
01594 }
01595 Iarg++ ; continue ;
01596 }
01597
01598
01599
01600 if( strncmp(Argv[Iarg],"-weight",4) == 0 ){
01601 int bb,ii ; char * cpt ;
01602 THD_3dim_dataset * wset ;
01603 char dname[256] ;
01604
01605 if( VL_imwt != NULL ){
01606 fprintf(stderr,"** Can't have two -weight options\n") ; exit(1) ;
01607 }
01608
01609
01610
01611 cpt = strstr( Argv[++Iarg] , "[" ) ;
01612 if( cpt == NULL || cpt == Argv[Iarg] ){
01613 fprintf(stderr,"** Illegal weight dataset after -weight\n"); exit(1);
01614 }
01615 ii = cpt - Argv[Iarg] ;
01616 memcpy(dname,Argv[Iarg],ii) ; dname[ii] = '\0' ;
01617 bb = -1 ; sscanf( cpt+1 , "%d" , &bb ) ;
01618 if( bb < 0 ){
01619 fprintf(stderr,"** Illegal sub-brick selector after -weight\n"); exit(1);
01620 }
01621 wset = THD_open_one_dataset( dname ) ;
01622 if( wset == NULL ){
01623 fprintf(stderr,"** Can't open weight dataset %s\n",dname); exit(1);
01624 }
01625 if( bb >= DSET_NVALS(wset) ){
01626 fprintf(stderr,"** Illegal sub-brick selector for dataset %s\n",dname);
01627 exit(1) ;
01628 }
01629 if( VL_verbose )
01630 fprintf(stderr,"++ Reading in weight dataset %s\n",DSET_BRIKNAME(wset)) ;
01631 DSET_load(wset) ;
01632 VL_imwt = mri_to_float( DSET_BRICK(wset,bb) ) ;
01633 DSET_delete( wset ) ;
01634 Iarg++ ; continue ;
01635 }
01636
01637
01638
01639 if( strcmp(Argv[Iarg],"-twodup") == 0 ){
01640 VL_twodup++ ; Iarg++ ; continue ;
01641 }
01642
01643
01644
01645 if( strcmp(Argv[Iarg],"-tshift") == 0 ){
01646 VL_tshift = 1 ;
01647 VL_tshift_ignore = (int) strtod(Argv[++Iarg],NULL) ;
01648 if( VL_tshift_ignore < 0 ) ERREX("-tshift parameter is negative!") ;
01649 Iarg++ ; continue ;
01650 }
01651
01652
01653
01654 if( strncmp(Argv[Iarg],"-rotpar",7) == 0 ){
01655 ATR_float * atr ;
01656
01657 if( VL_rotpar_dset != NULL )
01658 ERREX("Can't use -2 -rotparent options!") ;
01659
01660 VL_rotpar_dset = THD_open_one_dataset( Argv[++Iarg] ) ;
01661 if( VL_rotpar_dset == NULL )
01662 ERREX("Can't open -rotparent dataset!\n") ;
01663
01664 atr = THD_find_float_atr( VL_rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ;
01665 if( atr == NULL || atr->nfl < 12 )
01666 ERREX("-rotparent dataset doesn't have VOLREG attributes!?") ;
01667
01668 Iarg++ ; continue ;
01669 }
01670
01671 if( strncmp(Argv[Iarg],"-gridpar",7) == 0 ){
01672
01673 if( VL_gridpar_dset != NULL )
01674 ERREX("Can't use -2 -gridparent options!") ;
01675
01676 VL_gridpar_dset = THD_open_one_dataset( Argv[++Iarg] ) ;
01677 if( VL_gridpar_dset == NULL )
01678 ERREX("Can't open -gridparent dataset!\n") ;
01679
01680 Iarg++ ; continue ;
01681 }
01682
01683
01684
01685 fprintf(stderr,"** Unknown option: %s\a\n",Argv[Iarg]) ; exit(1) ;
01686 }
01687
01688
01689
01690
01691 if( VL_twodup ){
01692 if( !VL_twopass || VL_intern ) VL_twodup == 0 ;
01693 }
01694
01695
01696
01697 if( VL_gridpar_dset != NULL && VL_rotpar_dset == NULL ){
01698 fprintf(stderr,
01699 "++ WARNING: -gridparent means nothing without -rotparent!\n");
01700 DSET_delete( VL_gridpar_dset ) ;
01701 VL_gridpar_dset = NULL ;
01702 }
01703
01704 if( VL_rotpar_dset != NULL && VL_twopass ){
01705 fprintf(stderr,
01706 "++ WARNING: Combining -rotparent and -twopass isn't recommended!\n");
01707 }
01708
01709 if( VL_gridpar_dset != NULL && VL_twodup ){
01710 fprintf(stderr,"++ WARNING: -gridparent overrides -twodup!\n") ;
01711 VL_twodup = 0 ;
01712 }
01713
01714
01715
01716 if( VL_edperc < 0 ){
01717 char *ef=my_getenv("AFNI_VOLREG_EDGING") , *eq ;
01718 float ee ;
01719
01720 if( ef == NULL ){
01721 VL_edperc = 1 ; VL_edging = 5.0 ;
01722 } else {
01723 ee = strtod(ef,&eq) ;
01724 if( ee < 0 ){
01725 fprintf(stderr,"** Illegal value in AFNI_VOLREG_EDGING\n"); exit(1);
01726 }
01727 if( *eq == '%' ){
01728 if( ee > 25.0 ){
01729 fprintf(stderr,"** Illegal percentage in AFNI_VOLREG_EDGING\n"); exit(1);
01730 }
01731 VL_edperc = 1 ; VL_edging = ee ;
01732 } else {
01733 VL_edperc = 0 ; VL_edging = ee ;
01734 }
01735 }
01736 }
01737
01738
01739
01740 if( Iarg > Argc ){
01741 fprintf(stderr,"** Too few arguments!? Last=%s\n",Argv[Argc-1]) ; exit(1) ;
01742 } else if( Iarg < Argc-1 ){
01743 fprintf(stderr,"** Too many arguments?! Dataset=%s?\n",Argv[Iarg]) ; exit(1) ;
01744 }
01745
01746 VL_dset = THD_open_dataset( Argv[Iarg] ) ;
01747
01748
01749
01750 if( VL_dset == NULL ){
01751 fprintf(stderr,"** Can't open dataset %s\n",Argv[Iarg]) ; exit(1) ;
01752 }
01753
01754 if( VL_tshift ){
01755 if( DSET_NVALS(VL_dset) < 2 ){
01756 fprintf(stderr,"++ WARNING: -tshift used on a 1 brick dataset!\n") ;
01757 VL_tshift = 0 ;
01758 } else if( VL_dset->taxis == NULL || VL_dset->taxis->nsl < DSET_NZ(VL_dset) ){
01759 fprintf(stderr,"++ WARNING: -tshift used on a dataset with no time-offsets!\n") ;
01760 VL_tshift = 0 ;
01761 } else if( VL_tshift_ignore > DSET_NVALS(VL_dset)-5 ){
01762 fprintf(stderr,"++ WARNING: -tshift ignore is too large for this dataset!\n") ;
01763 VL_tshift = 0 ;
01764 } else if( VL_imbase == NULL && VL_nbase < VL_tshift_ignore ){
01765 fprintf(stderr,"++ WARNING: base brick is prior to -tshift ignore point.\n") ;
01766 }
01767 }
01768
01769 if( VL_imbase == NULL && VL_nbase >= DSET_NVALS(VL_dset) ){
01770 fprintf(stderr,"** Dataset %s doesn't have base brick index = %d\n",
01771 Argv[Iarg] , VL_nbase ) ;
01772 exit(1) ;
01773 }
01774
01775
01776 #if 1
01777 if( VL_bset != NULL ){
01778 int mm = THD_dataset_mismatch( VL_dset , VL_bset ) , nn=0 ;
01779
01780 if( mm & MISMATCH_DIMEN ){
01781 fprintf(stderr,
01782 "** Input %s and base %s don't have same dimensions!\n"
01783 " Input: nx=%d ny=%d nz=%d\n"
01784 " Base: nx=%d ny=%d nz=%d\n",
01785 DSET_HEADNAME(VL_dset) , DSET_HEADNAME(VL_bset) ,
01786 DSET_NX(VL_dset) , DSET_NY(VL_dset) , DSET_NZ(VL_dset) ,
01787 DSET_NX(VL_bset) , DSET_NY(VL_bset) , DSET_NZ(VL_bset) ) ;
01788 nn++ ;
01789 }
01790
01791 if( mm & MISMATCH_DELTA ){
01792 fprintf(stderr,
01793 "** Input %s and base %s don't have same grid spacing!\n"
01794 " Input: dx=%6.3f dy=%6.3f dz=%6.3f\n"
01795 " Base: dx=%6.3f dy=%6.3f dz=%6.3f\n",
01796 DSET_HEADNAME(VL_dset) , DSET_HEADNAME(VL_bset) ,
01797 DSET_DX(VL_dset) , DSET_DY(VL_dset) , DSET_DZ(VL_dset) ,
01798 DSET_DX(VL_bset) , DSET_DY(VL_bset) , DSET_DZ(VL_bset) ) ;
01799 nn++ ;
01800 }
01801
01802 if( mm & MISMATCH_ORIENT ){
01803 fprintf(stderr,
01804 "** Input %s and base %s don't have same orientation!\n"
01805 " Input: %s %s %s\n"
01806 " Base: %s %s %s \n" ,
01807 DSET_HEADNAME(VL_dset) , DSET_HEADNAME(VL_bset) ,
01808 ORIENT_shortstr[VL_dset->daxes->xxorient] ,
01809 ORIENT_shortstr[VL_dset->daxes->yyorient] ,
01810 ORIENT_shortstr[VL_dset->daxes->zzorient] ,
01811 ORIENT_shortstr[VL_bset->daxes->xxorient] ,
01812 ORIENT_shortstr[VL_bset->daxes->yyorient] ,
01813 ORIENT_shortstr[VL_bset->daxes->zzorient] ) ;
01814 nn++ ;
01815 }
01816
01817 if( nn > 0 ){
01818 fprintf(stderr,
01819 "** FATAL ERROR: perhaps you could make your datasets match?\n") ;
01820 exit(1) ;
01821 }
01822 }
01823 #else
01824 if( VL_imbase != NULL && ( VL_imbase->nx != DSET_NX(VL_dset) ||
01825 VL_imbase->ny != DSET_NY(VL_dset) ||
01826 VL_imbase->nz != DSET_NZ(VL_dset) ) ){
01827
01828 fprintf(stderr,"** Dataset %s doesn't conform to dimensions of base brick\n",
01829 Argv[Iarg] ) ;
01830 fprintf(stderr," Base has nx = %d ny = %d nz = %d\n",
01831 VL_imbase->nx, VL_imbase->ny, VL_imbase->nz ) ;
01832 fprintf(stderr," Dataset has nx = %d ny = %d nz = %d\n",
01833 DSET_NX(VL_dset) , DSET_NY(VL_dset) , DSET_NZ(VL_dset) ) ;
01834 exit(1) ;
01835 }
01836
01837 if( VL_imbase != NULL &&
01838 ( fabs(DSET_DX(VL_dset)) != bdx ||
01839 fabs(DSET_DY(VL_dset)) != bdy ||
01840 fabs(DSET_DZ(VL_dset)) != bdz ) ){
01841
01842 fprintf(stderr,"++ WARNING:\n"
01843 "++ Dataset %s and base have different grid spacings:\n"
01844 "++ Dataset: dx=%9.3f dy=%9.3f dz=%9.3f\n"
01845 "++ Base: dx=%9.3f dy=%9.3f dz=%9.3f\n" ,
01846 Argv[Iarg] ,
01847 fabs(DSET_DX(VL_dset)),fabs(DSET_DY(VL_dset)),fabs(DSET_DZ(VL_dset)),
01848 bdx,bdy,bdz ) ;
01849 }
01850 #endif
01851
01852 if( VL_imwt != NULL && ( VL_imwt->nx != DSET_NX(VL_dset) ||
01853 VL_imwt->ny != DSET_NY(VL_dset) ||
01854 VL_imwt->nz != DSET_NZ(VL_dset) ) ){
01855
01856 fprintf(stderr,"** Dataset %s doesn't conform to dimensions of weight brick\n",
01857 Argv[Iarg] ) ;
01858 fprintf(stderr," Weight has nx = %d ny = %d nz = %d\n",
01859 VL_imwt->nx, VL_imwt->ny, VL_imwt->nz ) ;
01860 fprintf(stderr," Dataset has nx = %d ny = %d nz = %d\n",
01861 DSET_NX(VL_dset) , DSET_NY(VL_dset) , DSET_NZ(VL_dset) ) ;
01862 exit(1) ;
01863 }
01864
01865 if( VL_intern && DSET_NVALS(VL_dset) == 1 ){
01866 fprintf(stderr,"** You can't register a 1 brick dataset to itself!\n") ;
01867 exit(1) ;
01868 }
01869
01870
01871
01872 if( VL_twopass && VL_coarse_del > 0 && VL_coarse_num > 0 ){
01873 int mm ;
01874 mm = MIN( DSET_NX(VL_dset) , DSET_NY(VL_dset) ) ;
01875 mm = MIN( DSET_NZ(VL_dset) , mm ) ;
01876 mm = (int)(0.1*mm + 0.499) ;
01877 if( VL_coarse_del > mm ){
01878 fprintf(stderr,"++ Coarse del was %d, replaced with %d\n",VL_coarse_del,mm) ;
01879 VL_coarse_del = mm ;
01880 }
01881 }
01882
01883
01884
01885 if( VL_imwt != NULL && VL_wtinp ){
01886 fprintf(stderr,"++ Input weight file overrides -wtinp option!\n") ;
01887 VL_wtinp = 0 ;
01888 }
01889
01890
01891
01892 return ;
01893 }
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905 #define B(i,j,k) b[(i)+(j)*nx+(k)*nxy]
01906 #define V(i,j,k) v[(i)+(j)*nx+(k)*nxy]
01907
01908 float voldif( int nx, int ny, int nz, float *b,
01909 int dx, int dy, int dz, float *v, int edge )
01910 {
01911 int nxy=nx*ny, nxtop=nx-edge, nytop=ny-edge, nztop=nz-edge , ii,jj,kk ;
01912 float bbsum=0.0 , vvsum=0.0 , bvsum=0.0 , bb,vv ;
01913
01914 for( kk=edge ; kk < nztop ; kk++ ){
01915 for( jj=edge ; jj < nytop ; jj++ ){
01916 for( ii=edge ; ii < nxtop ; ii++ ){
01917 bb = B(ii,jj,kk) ; vv = V(ii-dx,jj-dy,kk-dz) ;
01918 bbsum += bb*bb ; vvsum += vv*vv ; bvsum += bb*vv ;
01919 }
01920 }
01921 }
01922
01923 if( vvsum > 0.0 ) bbsum -= bvsum*bvsum/vvsum ;
01924 return bbsum ;
01925 }
01926
01927
01928
01929
01930
01931
01932 void get_best_shift( int nx, int ny, int nz,
01933 float *b, float *v ,
01934 int *dxp , int *dyp , int *dzp )
01935 {
01936 int bdx=0 , bdy=0 , bdz=0 , dx,dy,dz , nxyz=nx*ny*nz ;
01937 float bsum , sum ;
01938
01939 int shift = VL_coarse_del, numsh = VL_coarse_num,
01940 shtop = shift*numsh , edge = shtop+shift , sqtop = shtop*shtop ;
01941
01942 bsum = 0.0 ;
01943 for( dx=0 ; dx < nxyz ; dx++ ) bsum += b[dx]*b[dx] ;
01944
01945 for( dz=-shtop ; dz <= shtop ; dz+=shift ){
01946 for( dy=-shtop ; dy <= shtop ; dy+=shift ){
01947 for( dx=-shtop ; dx <= shtop ; dx+=shift ){
01948 if( dx*dx+dy*dy+dz*dz > sqtop ) continue ;
01949 sum = voldif( nx,ny,nz , b , dx,dy,dz , v , edge ) ;
01950 if( sum < bsum ){
01951 #if 0
01952 fprintf(stderr," get_best_shift: bsum=%g dx=%d dy=%d dz=%d\n",sum,dx,dy,dz) ;
01953 #endif
01954 bsum = sum ; bdx = dx ; bdy = dy ; bdz = dz ;
01955 }
01956 }
01957 }
01958 }
01959
01960 *dxp = bdx ; *dyp = bdy ; *dzp = bdz ; return ;
01961 }
