00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "afni.h"
00008 
00009 #ifndef ALLOW_PLUGINS
00010 #  error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012 
00013 /***********************************************************************
00014   Plugin to do what 3dvolreg.c does.  Most of the code is
00015   adapted from plug_imreg.c, 3dvolreg.c, and plug_realtime.c.
00016 ************************************************************************/
00017 
00018 /*--------------------- string to 'help' the user --------------------*/
00019 
00020 static char helpstring[] =
00021   "Purpose: 3D (volumetric) registration of 3D+time dataset\n"
00022   "[See also program '3dvolreg']\n"
00023   "\n"
00024   "Input items to this plugin are:\n"
00025   "   Datasets:   Input  = 3D+time dataset to process\n"
00026   "               Output = Prefix for new dataset\n"
00027   "\n"
00028   "   Parameters: Base Brick = Time index for base image\n"
00029   "               Resampling = Resampling method\n"
00030   "\n"
00031   "   Basis:      Dataset = Dataset from which to take base brick\n"
00032   "                         [if not selected, uses input dataset]\n"
00033   "\n"
00034   "   Outputs:    dfile = If not blank, name of file in which to\n"
00035   "                       save estimated movement parameters plus\n"
00036   "                       before-and-after residuals;\n"
00037   "                       see '3dvolreg -help' for details\n"
00038   "              1Dfile = If not blank, name of file in which to save\n"
00039   "                       ONLY the estimated movement parameters;\n"
00040   "                       this file will also be added to the AFNI\n"
00041   "                       timeseries list, so that it can be used\n"
00042   "                       as an 'ort' with FIM.\n"
00043   "               Graph = If 'Yes', will plot the estimated movement\n"
00044   "                       parameters when the registration is done.\n"
00045   "\n"
00046   "WARNING: This plugin can be slow, particularly when registering\n"
00047   "         large datasets using Fourier resampling.  The algorithm\n"
00048   "         used will only work to correct for small movements.\n"
00049   "\n"
00050   "Author -- RW Cox, November 1998"
00051 ;
00052 
00053 /*----------------- prototypes for internal routines -----------------*/
00054 
00055 char * VOLREG_main( PLUGIN_interface * ) ;  /* the entry point */
00056 
00057 /*---------------------------- global data ---------------------------*/
00058 
00059 static PLUGIN_interface * global_plint = NULL ;
00060 
00061 #define NRESAM 4
00062 static char * REG_resam_strings[NRESAM] = {
00063    "Cubic" , "Quintic" , "Heptic" , "Fourier"
00064  } ;
00065 
00066 static int REG_resam_ints[NRESAM] = {
00067    MRI_CUBIC , MRI_QUINTIC , MRI_HEPTIC , MRI_FOURIER
00068 } ;
00069 
00070 #define NGRAPH 2
00071 static char * GRAPH_strings[NGRAPH] = { "No" , "Yes" } ;
00072 
00073 static int         VL_nbase  = 3 ;
00074 static int         VL_resam  = 3 ;
00075 static int         VL_clipit = 1 ;
00076 static int         VL_intern = 1 ;
00077 static int         VL_graph  = 0 ;
00078 static MRI_IMAGE * VL_imbase = NULL ;
00079 
00080 static char VL_dfile[256]  = "\0" ;
00081 static char VL_1Dfile[256] = "\0" ;              /* 14 Apr 2000 */
00082 static char VL_prefix[THD_MAX_PREFIX] = "\0" ;
00083 
00084 static THD_3dim_dataset * VL_dset = NULL ;
00085 
00086 static int VL_maxite = 9 ;
00087 static float VL_dxy  = 0.05 ;  /* voxels */
00088 static float VL_dph  = 0.07 ;  /* degrees */
00089 static float VL_del  = 0.70 ;  /* voxels */
00090 
00091 /***********************************************************************
00092    Set up the interface to the user:
00093     1) Create a new interface using "PLUTO_new_interface";
00094 
00095     2) For each line of inputs, create the line with "PLUTO_add_option"
00096          (this line of inputs can be optional or mandatory);
00097 
00098     3) For each item on the line, create the item with
00099         "PLUTO_add_dataset" for a dataset chooser,
00100         "PLUTO_add_string"  for a string chooser,
00101         "PLUTO_add_number"  for a number chooser.
00102 ************************************************************************/
00103 
00104 
00105 DEFINE_PLUGIN_PROTOTYPE
00106 
00107 PLUGIN_interface * PLUGIN_init( int ncall )
00108 {
00109    PLUGIN_interface * plint ;     /* will be the output of this routine */
00110 
00111    if( ncall > 0 ) return NULL ;  /* only one interface */
00112 
00113    /*---------------- set titles and call point ----------------*/
00114 
00115    plint = PLUTO_new_interface( "3D Registration" ,
00116                                 "3D Registration of a 3D+time Dataset" ,
00117                                 helpstring ,
00118                                 PLUGIN_CALL_VIA_MENU , VOLREG_main  ) ;
00119 
00120    PLUTO_add_hint( plint , "3D Registration of a 3D+time Dataset" ) ;
00121 
00122    global_plint = plint ;  /* make global copy */
00123 
00124    PLUTO_set_sequence( plint , "A:newdset:reg" ) ;
00125 
00126    /*--------- 1st line ---------*/
00127 
00128    PLUTO_add_option( plint ,
00129                      "Datasets" ,  /* label at left of input line */
00130                      "DAtasets" ,  /* tag to return to plugin */
00131                      TRUE          /* is this mandatory? */
00132                    ) ;
00133 
00134    PLUTO_add_dataset(  plint ,
00135                        "Input" ,          /* label next to button   */
00136                        ANAT_ALL_MASK ,    /* take any anat datasets */
00137                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00138                        DIMEN_ALL_MASK | BRICK_ALLREAL_MASK
00139                     ) ;
00140 
00141    PLUTO_add_string( plint ,
00142                      "Output" ,  /* label next to textfield */
00143                      0,NULL ,    /* no fixed strings to choose among */
00144                      19          /* 19 spaces for typing in value */
00145                    ) ;
00146 
00147    /*---------- 2nd line --------*/
00148 
00149    PLUTO_add_option( plint ,
00150                      "Parameters" ,  /* label at left of input line */
00151                      "Parameters" ,  /* tag to return to plugin */
00152                      TRUE            /* is this mandatory? */
00153                    ) ;
00154 
00155    PLUTO_add_number( plint ,
00156                      "Base Brick" ,  /* label next to chooser */
00157                      0 ,             /* smallest possible value */
00158                      98 ,            /* largest possible value */
00159                      0 ,             /* decimal shift (none in this case) */
00160                      VL_nbase ,      /* default value */
00161                      FALSE           /* allow user to edit value? */
00162                    ) ;
00163 
00164    PLUTO_add_string( plint, "Resampling", NRESAM, REG_resam_strings, VL_resam ) ;
00165 
00166    /*---------- 3rd line --------*/
00167 
00168    PLUTO_add_option( plint , "Basis" , "Basis" , FALSE ) ;
00169 
00170    PLUTO_add_dataset(  plint ,
00171                        "Dataset" ,        /* label next to button   */
00172                        ANAT_ALL_MASK ,    /* take any anat datasets */
00173                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00174                        DIMEN_ALL_MASK | BRICK_ALLREAL_MASK
00175                     ) ;
00176 
00177    /*---------- 4th line ---------*/
00178 
00179    PLUTO_add_option( plint , "Outputs" , "Outputs" , FALSE ) ;
00180 
00181    PLUTO_add_string( plint ,
00182                      "dfile" ,   /* label next to textfield */
00183                      0,NULL ,    /* no fixed strings to choose among */
00184                      19          /* 19 spaces for typing in value */
00185                    ) ;
00186 
00187    /* 14 Apr 2000: add 1Dfile */
00188 
00189    PLUTO_add_string( plint ,
00190                      "1Dfile" ,  /* label next to textfield */
00191                      0,NULL ,    /* no fixed strings to choose among */
00192                      19          /* 19 spaces for typing in value */
00193                    ) ;
00194 
00195    PLUTO_add_string( plint , "Graph" , NGRAPH , GRAPH_strings , VL_graph ) ;
00196 
00197    /*--------- done with interface setup ---------*/
00198 
00199    return plint ;
00200 }
00201 
00202 /***************************************************************************
00203   Main routine for this plugin (will be called from AFNI).
00204   If the return string is not NULL, some error transpired, and
00205   AFNI will popup the return string in a message box.
00206 ****************************************************************************/
00207 
00208 #undef FREEUP
00209 #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
00210 
00211 char * VOLREG_main( PLUGIN_interface * plint )
00212 {
00213    MCW_idcode * idc ;
00214    char * str ;
00215    MRI_3dalign_basis * albase ;
00216    THD_3dim_dataset * new_dset ;
00217    MRI_IMAGE * qim , * tim , * fim ;
00218    float *dx, *dy, *dz, *roll, *yaw, *pitch, *rmsnew, *rmsold, *imb, *tar ;
00219    float ddx,ddy,ddz , sum ;
00220    float dxtop,dytop,dztop , rolltop,yawtop,pitchtop ;
00221    float dxbot,dybot,dzbot , rollbot,yawbot,pitchbot ;
00222    float dxbar,dybar,dzbar , rollbar,yawbar,pitchbar ;
00223    int kim,ii , imcount , iha , ax1,ax2,ax3 , hax1,hax2,hax3 ;
00224 
00225    double cputim = COX_cpu_time();
00226 
00227    /*--------------------------------------------------------------------*/
00228    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00229 
00230    /*--------- go to first input line ---------*/
00231 
00232    PLUTO_next_option(plint) ;
00233 
00234    idc     = PLUTO_get_idcode(plint) ;   /* get dataset item */
00235    VL_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00236    if( VL_dset == NULL )
00237       return "*************************\n"
00238              "Cannot find Input Dataset\n"
00239              "*************************"  ;
00240 
00241    ii = DSET_NVALS_PER_TIME(VL_dset) ;
00242    if( ii > 1 )
00243       return "************************************\n"
00244              "Dataset has > 1 value per time point\n"
00245              "************************************"  ;
00246 
00247    str = PLUTO_get_string(plint) ;   /* get the output prefix */
00248    if( ! PLUTO_prefix_ok(str) )      /* check if it is OK */
00249       return "************************\n"
00250              "Output Prefix is illegal\n"
00251              "************************"  ;
00252    strcpy( VL_prefix , str ) ;
00253 
00254    /*--------- go to next input line ---------*/
00255 
00256    PLUTO_next_option(plint) ;
00257 
00258    VL_nbase = PLUTO_get_number(plint) ;
00259 
00260    str      = PLUTO_get_string(plint) ;
00261    VL_resam = PLUTO_string_index( str , NRESAM , REG_resam_strings ) ;
00262    VL_resam = REG_resam_ints[VL_resam] ;
00263 
00264    /*--------- see if the other (non-mandatory) option line are present --------*/
00265 
00266    VL_imbase   = NULL ;
00267    VL_dfile[0] = '\0' ;
00268    VL_1Dfile[0]= '\0' ;  /* 14 Apr 2000 */
00269    VL_graph    = 0 ;
00270 
00271    while( 1 ){
00272       str = PLUTO_get_optiontag( plint ) ; if( str == NULL ) break ;
00273 
00274       if( strcmp(str,"Basis") == 0 ){  /* Alternate basis dataset */
00275          THD_3dim_dataset * bset ;
00276 
00277          idc  = PLUTO_get_idcode(plint) ;   /* get dataset item */
00278          bset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00279          if( bset == NULL )
00280             return "*************************\n"
00281                    "Cannot find Basis Dataset\n"
00282                    "*************************"  ;
00283 
00284           if( VL_nbase >= DSET_NVALS(bset) || VL_nbase < 0 )
00285              return "*****************************************\n"
00286                     "Base index out of range for Basis dataset\n"
00287                     "*****************************************"  ;
00288 
00289           if( DSET_NX(bset) != DSET_NX(VL_dset) ||
00290               DSET_NY(bset) != DSET_NY(VL_dset) || DSET_NZ(bset) != DSET_NZ(VL_dset) )
00291              return "***************************\n"
00292                     "Brick size mismatch between\n"
00293                     " Input and Basis datasets! \n"
00294                     "***************************"  ;
00295 
00296           DSET_load(bset) ;
00297           if( ! DSET_LOADED(bset) ) return "************************\n"
00298                                            "Can't load Basis dataset\n"
00299                                            "************************"  ;
00300 
00301           VL_imbase = mri_to_float( DSET_BRICK(bset,VL_nbase) ) ;  /* copy this */
00302           VL_intern = 0 ;
00303           DSET_unload(bset) ;
00304 
00305       } else if( strcmp(str,"Outputs") == 0 ){  /* Optional outputs */
00306 
00307          str = PLUTO_get_string(plint) ;
00308          if( str != NULL && str[0] != '\0' ){
00309             if( THD_filename_ok(str) ) strcpy( VL_dfile , str ) ;
00310             else                       return "*************************\n"
00311                                               "dfile name not acceptable\n"
00312                                               "*************************"  ;
00313          }
00314 
00315          /* 14 Apr 2000: get 1Dfile */
00316 
00317          str = PLUTO_get_string(plint) ;
00318          if( str != NULL && str[0] != '\0' ){
00319             if( THD_filename_ok(str) ) strcpy( VL_1Dfile , str ) ;
00320             else                       return "**************************\n"
00321                                               "1Dfile name not acceptable\n"
00322                                               "**************************"  ;
00323          }
00324 
00325          str      = PLUTO_get_string(plint) ;
00326          VL_graph = PLUTO_string_index( str , NGRAPH , GRAPH_strings ) ;
00327       }
00328 
00329    }
00330 
00331    /*-- must manufacture base image --*/
00332 
00333    if( VL_imbase == NULL ){
00334 
00335       if( DSET_NVALS(VL_dset) < 2 )
00336          return "******************************************\n"
00337                 "Can't register a 1 brick dataset to itself\n"
00338                 "******************************************"  ;
00339 
00340       if( VL_nbase >= DSET_NVALS(VL_dset) || VL_nbase < 0 )
00341           return "*****************************************\n"
00342                  "Base index out of range for Input dataset\n"
00343                  "*****************************************"  ;
00344 
00345       DSET_load(VL_dset) ;
00346       if( ! DSET_LOADED(VL_dset) ) return "************************\n"
00347                                           "Can't load Input dataset\n"
00348                                           "************************"  ;
00349 
00350       VL_imbase = mri_to_float( DSET_BRICK(VL_dset,VL_nbase) ) ;  /* copy this */
00351       VL_intern = 1 ;
00352    }
00353 
00354    VL_imbase->dx = fabs( DSET_DX(VL_dset) ) ;  /* must set the voxel dimensions */
00355    VL_imbase->dy = fabs( DSET_DY(VL_dset) ) ;
00356    VL_imbase->dz = fabs( DSET_DZ(VL_dset) ) ;
00357    imb = MRI_FLOAT_PTR( VL_imbase ) ;          /* need this to compute rms */
00358 
00359    /*------------- set up to compute new dataset -----------*/
00360 
00361    DSET_load( VL_dset ) ;
00362    if( ! DSET_LOADED(VL_dset) ){
00363       mri_free( VL_imbase ) ;
00364       return "************************\n"
00365              "Can't load Input dataset\n"
00366              "************************"  ;
00367    }
00368 
00369    imcount = DSET_NVALS( VL_dset ) ;
00370    dx      = (float *) malloc( sizeof(float) * imcount ) ;
00371    dy      = (float *) malloc( sizeof(float) * imcount ) ;
00372    dz      = (float *) malloc( sizeof(float) * imcount ) ;
00373    roll    = (float *) malloc( sizeof(float) * imcount ) ;
00374    pitch   = (float *) malloc( sizeof(float) * imcount ) ;
00375    yaw     = (float *) malloc( sizeof(float) * imcount ) ;
00376    rmsnew  = (float *) malloc( sizeof(float) * imcount ) ;
00377    rmsold  = (float *) malloc( sizeof(float) * imcount ) ;
00378 
00379    iha = THD_handedness( VL_dset )   ;                     /* LH or RH? */
00380    ax1 = THD_axcode( VL_dset , 'I' ) ; hax1 = ax1 * iha ;  /* roll */
00381    ax2 = THD_axcode( VL_dset , 'R' ) ; hax2 = ax2 * iha ;  /* pitch */
00382    ax3 = THD_axcode( VL_dset , 'A' ) ; hax3 = ax3 * iha ;  /* yaw */
00383 
00384    mri_3dalign_params( VL_maxite , VL_dxy , VL_dph , VL_del ,
00385                        abs(ax1)-1 , abs(ax2)-1 , abs(ax3)-1 , -1 ) ;
00386 
00387    mri_3dalign_method( VL_resam , 0 , 0 , VL_clipit ) ;
00388 
00389    new_dset = EDIT_empty_copy( VL_dset ) ;
00390    EDIT_dset_items( new_dset , ADN_prefix , VL_prefix , ADN_none ) ;
00391 
00392    { char * his = PLUTO_commandstring(plint) ;
00393      tross_Copy_History( VL_dset , new_dset ) ;
00394      tross_Append_History( new_dset , his ) ; free(his) ;
00395    }
00396 
00397    albase = mri_3dalign_setup( VL_imbase , NULL ) ;
00398 
00399    /*-- loop over sub-bricks and register them --*/
00400 
00401    PLUTO_popup_meter(plint) ;
00402 
00403    dxbar = dybar = dzbar = rollbar = yawbar = pitchbar = 0.0 ;
00404 
00405    for( kim=0 ; kim < imcount ; kim++ ){
00406 
00407       qim     = DSET_BRICK( VL_dset , kim ) ; /* the sub-brick in question */
00408       fim     = mri_to_float( qim ) ;         /* make a float copy */
00409       fim->dx = fabs( DSET_DX(VL_dset) ) ;    /* must set voxel dimensions */
00410       fim->dy = fabs( DSET_DY(VL_dset) ) ;
00411       fim->dz = fabs( DSET_DZ(VL_dset) ) ;
00412 
00413       /*-- the actual registration --*/
00414 
00415       if( !VL_intern || kim != VL_nbase ){
00416          tim = mri_3dalign_one( albase , fim ,
00417                                 roll+kim , pitch+kim , yaw+kim ,
00418                                 &ddx     , &ddy      , &ddz     ) ;
00419       } else {
00420          tim = mri_to_float( VL_imbase ) ;
00421          roll[kim] = pitch[kim] = yaw[kim] = ddx = ddy = ddz = 0.0 ;
00422       }
00423 
00424       /*-- need to massage output parameters --*/
00425 
00426       roll[kim]  *= (180.0/PI) ; if( hax1 < 0 ) roll[kim]  = -roll[kim] ;
00427       pitch[kim] *= (180.0/PI) ; if( hax2 < 0 ) pitch[kim] = -pitch[kim] ;
00428       yaw[kim]   *= (180.0/PI) ; if( hax3 < 0 ) yaw[kim]   = -yaw[kim] ;
00429 
00430       switch( new_dset->daxes->xxorient ){
00431          case ORI_R2L_TYPE: dy[kim] =  ddx ; break ;
00432          case ORI_L2R_TYPE: dy[kim] = -ddx ; break ;
00433          case ORI_P2A_TYPE: dz[kim] = -ddx ; break ;
00434          case ORI_A2P_TYPE: dz[kim] =  ddx ; break ;
00435          case ORI_I2S_TYPE: dx[kim] =  ddx ; break ;
00436          case ORI_S2I_TYPE: dx[kim] = -ddx ; break ;
00437       }
00438 
00439       switch( new_dset->daxes->yyorient ){
00440          case ORI_R2L_TYPE: dy[kim] =  ddy ; break ;
00441          case ORI_L2R_TYPE: dy[kim] = -ddy ; break ;
00442          case ORI_P2A_TYPE: dz[kim] = -ddy ; break ;
00443          case ORI_A2P_TYPE: dz[kim] =  ddy ; break ;
00444          case ORI_I2S_TYPE: dx[kim] =  ddy ; break ;
00445          case ORI_S2I_TYPE: dx[kim] = -ddy ; break ;
00446       }
00447 
00448       switch( new_dset->daxes->zzorient ){
00449          case ORI_R2L_TYPE: dy[kim] =  ddz ; break ;
00450          case ORI_L2R_TYPE: dy[kim] = -ddz ; break ;
00451          case ORI_P2A_TYPE: dz[kim] = -ddz ; break ;
00452          case ORI_A2P_TYPE: dz[kim] =  ddz ; break ;
00453          case ORI_I2S_TYPE: dx[kim] =  ddz ; break ;
00454          case ORI_S2I_TYPE: dx[kim] = -ddz ; break ;
00455       }
00456 
00457       /*-- collect statistics --*/
00458 
00459       if( kim == 0 ){
00460          dxtop    = dxbot    = dx[kim]    ;
00461          dytop    = dybot    = dy[kim]    ;
00462          dztop    = dzbot    = dz[kim]    ;
00463          rolltop  = rollbot  = roll[kim]  ;
00464          pitchtop = pitchbot = pitch[kim] ;
00465          yawtop   = yawbot   = yaw[kim]   ;
00466       } else {
00467          dxtop    = MAX(dxtop   , dx[kim]   ); dxbot    = MIN(dxbot   , dx[kim]   );
00468          dytop    = MAX(dytop   , dy[kim]   ); dybot    = MIN(dybot   , dy[kim]   );
00469          dztop    = MAX(dztop   , dz[kim]   ); dzbot    = MIN(dzbot   , dz[kim]   );
00470          rolltop  = MAX(rolltop , roll[kim] ); rollbot  = MIN(rollbot , roll[kim] );
00471          pitchtop = MAX(pitchtop, pitch[kim]); pitchbot = MIN(pitchbot, pitch[kim]);
00472          yawtop   = MAX(yawtop  , yaw[kim]  ); yawbot   = MIN(yawbot  , yaw[kim]  );
00473       }
00474 
00475       dxbar   += dx[kim]   ; dybar    += dy[kim]    ; dzbar  += dz[kim]  ;
00476       rollbar += roll[kim] ; pitchbar += pitch[kim] ; yawbar += yaw[kim] ;
00477 
00478       sum = 0.0 ;
00479       tar = MRI_FLOAT_PTR(tim) ;
00480       for( ii=0 ; ii < tim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
00481       rmsnew[kim] = sqrt( sum / tim->nvox ) ;
00482 
00483       sum = 0.0 ;
00484       tar = MRI_FLOAT_PTR(fim) ;
00485       for( ii=0 ; ii < fim->nvox ; ii++ ) sum += SQR( imb[ii] - tar[ii] ) ;
00486       rmsold[kim] = sqrt( sum / fim->nvox ) ;
00487 
00488       mri_free(fim) ;  /* only needed this to compute rmsold */
00489 
00490       /*-- Attach the registered brick to output dataset,
00491            converting it to the correct type, if necessary
00492            (the new registered brick in "tim" is stored as floats). --*/
00493 
00494       switch( qim->kind ){
00495 
00496          case MRI_float:
00497             EDIT_substitute_brick( new_dset , kim , MRI_float , MRI_FLOAT_PTR(tim) ) ;
00498             mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
00499          break ;
00500 
00501          case MRI_short:
00502             fim = mri_to_short(1.0,tim) ; mri_free( tim ) ;
00503             EDIT_substitute_brick( new_dset , kim , MRI_short , MRI_SHORT_PTR(fim) ) ;
00504             mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
00505          break ;
00506 
00507          case MRI_byte:
00508             fim = mri_to_byte(tim) ; mri_free( tim ) ;
00509             EDIT_substitute_brick( new_dset , kim , MRI_byte , MRI_BYTE_PTR(fim) ) ;
00510             mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
00511          break ;
00512 
00513          /*-- should not ever get here, but who knows? --*/
00514 
00515          default:
00516             mri_free( VL_imbase ) ;
00517             FREEUP(dx)     ; FREEUP(dy)     ; FREEUP(dz) ;
00518             FREEUP(roll)   ; FREEUP(yaw)    ; FREEUP(pitch) ;
00519             FREEUP(rmsnew) ; FREEUP(rmsold) ;
00520 
00521             fprintf(stderr,"\n*** Can't register bricks of type %s\n",
00522                     MRI_TYPE_name[qim->kind] ) ;
00523             return("*** Can't register bricks of this type");
00524             /* EXIT(1) ;*/
00525       }
00526 
00527       DSET_unload_one( VL_dset , kim ) ;      /* don't need this anymore */
00528 
00529       PLUTO_set_meter(plint, (100*(kim+1))/imcount ) ;
00530 
00531    }  /* end of loop over sub-bricks */
00532 
00533    /*-- done with registration --*/
00534 
00535    mri_3dalign_cleanup( albase ) ;
00536    mri_free( VL_imbase ) ;
00537    DSET_unload( VL_dset ) ;
00538 
00539    /*-- show some statistics of the results --*/
00540 
00541    dxbar   /= imcount ; dybar    /= imcount ; dzbar  /= imcount ;
00542    rollbar /= imcount ; pitchbar /= imcount ; yawbar /= imcount ;
00543 
00544    { char buf[512] ; int nbuf ;
00545 
00546       cputim = COX_cpu_time() - cputim ;
00547       sprintf(buf,"  \nCPU time for realignment=%.3g s\n\n" , cputim) ;
00548       nbuf = strlen(buf) ;
00549 
00550       sprintf(buf+nbuf,"Min : roll=%+.3f  pitch=%+.3f  yaw=%+.3f"
00551                        "  dS=%+.3f  dL=%+.3f  dP=%+.3f\n\n" ,
00552               rollbot , pitchbot , yawbot , dxbot , dybot , dzbot ) ;
00553       nbuf = strlen(buf) ;
00554 
00555       sprintf(buf+nbuf,"Mean: roll=%+.3f  pitch=%+.3f  yaw=%+.3f"
00556                        "  dS=%+.3f  dL=%+.3f  dP=%+.3f\n\n" ,
00557               rollbar , pitchbar , yawbar , dxbar , dybar , dzbar ) ;
00558       nbuf = strlen(buf) ;
00559 
00560       sprintf(buf+nbuf,"Max : roll=%+.3f  pitch=%+.3f  yaw=%+.3f"
00561                        "  dS=%+.3f  dL=%+.3f  dP=%+.3f\n" ,
00562               rolltop , pitchtop , yawtop , dxtop , dytop , dztop ) ;
00563 
00564       PLUTO_popup_message( plint , buf ) ;
00565    }
00566 
00567    /*-- save to a file --*/
00568 
00569    if( VL_dfile[0] != '\0' ){
00570       FILE * fp ;
00571 
00572       if( THD_is_file(VL_dfile) )
00573          PLUTO_popup_transient( plint , "** Warning:\n"
00574                                         "** Overwriting dfile" ) ;
00575 
00576       fp = fopen( VL_dfile , "w" ) ;
00577       for( kim=0 ; kim < imcount ; kim++ )
00578          fprintf(fp , "%4d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f  %11.4g %11.4g\n" ,
00579                  kim , roll[kim], pitch[kim], yaw[kim],
00580                        dx[kim], dy[kim], dz[kim], rmsold[kim] , rmsnew[kim]  ) ;
00581       fclose(fp) ;
00582    }
00583 
00584    if( VL_1Dfile[0] != '\0' ){  /* 14 Apr 2000 */
00585       FILE * fp ;
00586       char fn[256] ;
00587       MRI_IMAGE * tsim ;
00588       float * tsar ;
00589 
00590       strcpy(fn,VL_1Dfile) ;
00591       if( strstr(VL_1Dfile,"1D") == NULL ) strcat(fn,".1D") ;
00592       
00593       if( THD_is_file(fn) )
00594          PLUTO_popup_transient( plint , "** Warning:\n"
00595                                         "** Overwriting 1Dfile" ) ;
00596 
00597       tsim = mri_new( imcount , 6 , MRI_float ) ;
00598       tsar = MRI_FLOAT_PTR(tsim) ;
00599       fp = fopen( fn , "w" ) ;
00600       for( kim=0 ; kim < imcount ; kim++ ){
00601          fprintf(fp , "%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n" ,
00602                  roll[kim], pitch[kim], yaw[kim],
00603                  dx[kim]  , dy[kim]   , dz[kim]  ) ;
00604          tsar[0*imcount+kim] = roll[kim] ;
00605          tsar[1*imcount+kim] = pitch[kim] ;
00606          tsar[2*imcount+kim] = yaw[kim] ;
00607          tsar[3*imcount+kim] = dx[kim] ;
00608          tsar[4*imcount+kim] = dy[kim] ;
00609          tsar[5*imcount+kim] = dz[kim] ;
00610       }
00611       fclose(fp) ;
00612       PLUTO_register_timeseries( fn , tsim ) ;
00613       mri_free(tsim) ;
00614    }
00615 
00616    /*-- graph --*/
00617 
00618    if( VL_graph && imcount > 1 ){
00619       float * yar[7] , dt ;
00620       int nn = imcount ;
00621       static char * nar[6] = {
00622          "\\Delta I-S [mm]" , "\\Delta R-L [mm]" , "\\Delta A-P [mm]" ,
00623          "Roll [\\degree]" , "Pitch [\\degree]" , "Yaw [\\degree]"  } ;
00624 
00625       yar[0] = (float *) malloc( sizeof(float) * nn ) ;
00626       dt     = DSET_TIMESTEP(VL_dset) ; if( dt <= 0.0 ) dt = 1.0 ;
00627       for( ii=0 ; ii < nn ; ii++ ) yar[0][ii] = ii * dt ;
00628 
00629       yar[1] = dx   ; yar[2] = dy    ; yar[3] = dz  ;
00630       yar[4] = roll ; yar[5] = pitch ; yar[6] = yaw ;
00631 
00632       plot_ts_lab( GLOBAL_library.dc->display ,
00633                    nn , yar[0] , -6 , yar+1 ,
00634                    "time" , NULL , DSET_FILECODE(VL_dset) , nar , NULL ) ;
00635 
00636       free(yar[0]) ;
00637    }
00638 
00639    /* done with these statistics */
00640 
00641    FREEUP(dx)     ; FREEUP(dy)     ; FREEUP(dz) ;
00642    FREEUP(roll)   ; FREEUP(yaw)    ; FREEUP(pitch) ;
00643    FREEUP(rmsnew) ; FREEUP(rmsold) ;
00644 
00645    /*-- save new dataset to disk --*/
00646 
00647    PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00648    return NULL ;
00649 }