00001 
00002 /*****************************************************************************
00003    Major portions of this software are copyrighted by the Medical College
00004    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00005    License, Version 2.  See the file README.Copyright for details.
00006 ******************************************************************************/
00007 
00008 #include "afni.h"
00009 
00010 #ifndef ALLOW_PLUGINS
00011 #  error "Plugins not properly set up -- see machdep.h"
00012 #endif
00013 
00014 /***********************************************************************
00015   Plugin to do what imreg.c does
00016 ************************************************************************/
00017 
00018 /*--------------------- string to 'help' the user --------------------*/
00019 
00020 static char helpstring[] =
00021   "Purpose: 2D (slice-wise) registration of 3D+time dataset\n"
00022   "[See also programs 'imreg' and '2dImReg']\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   "   Parameters: Base   = Time index for base image\n"
00028   "----\n"
00029   "This option is only for experimenting with the parameters\n"
00030   "of the final (fine) step of the registration algorithm:\n"
00031   "   Fine Fit:   Blur   = FWHM of blurring prior to registration\n"
00032   "               Dxy    = Convergence tolerance for translations\n"
00033   "               Dphi   = Convergence tolerance for rotations\n"
00034   "Author -- RW Cox"
00035 ;
00036 
00037 /*----------------- prototypes for internal routines -----------------*/
00038 
00039 char * IMREG_main( PLUGIN_interface * ) ;  /* the entry point */
00040 
00041 /*---------------------------- global data ---------------------------*/
00042 
00043 static PLUGIN_interface * global_plint = NULL ;
00044 
00045 /***********************************************************************
00046    Set up the interface to the user:
00047     1) Create a new interface using "PLUTO_new_interface";
00048 
00049     2) For each line of inputs, create the line with "PLUTO_add_option"
00050          (this line of inputs can be optional or mandatory);
00051 
00052     3) For each item on the line, create the item with
00053         "PLUTO_add_dataset" for a dataset chooser,
00054         "PLUTO_add_string"  for a string chooser,
00055         "PLUTO_add_number"  for a number chooser.
00056 ************************************************************************/
00057 
00058 
00059 DEFINE_PLUGIN_PROTOTYPE
00060 
00061 PLUGIN_interface * PLUGIN_init( int ncall )
00062 {
00063    PLUGIN_interface * plint ;     /* will be the output of this routine */
00064 
00065    if( ncall > 0 ) return NULL ;  /* only one interface */
00066 
00067    /*---------------- set titles and call point ----------------*/
00068 
00069    plint = PLUTO_new_interface( "2D Registration" ,
00070                                 "2D Registration of a 3D+time Dataset" ,
00071                                 helpstring ,
00072                                 PLUGIN_CALL_VIA_MENU , IMREG_main  ) ;
00073 
00074    PLUTO_add_hint( plint , "2D Registration of a 3D+time Dataset" ) ;
00075 
00076    global_plint = plint ;  /* make global copy */
00077 
00078    PLUTO_set_sequence( plint , "A:newdset:reg" ) ;
00079 
00080    /*--------- 1st line ---------*/
00081 
00082    PLUTO_add_option( plint ,
00083                      "Datasets" ,  /* label at left of input line */
00084                      "DAtasets" ,  /* tag to return to plugin */
00085                      TRUE          /* is this mandatory? */
00086                    ) ;
00087 
00088    PLUTO_add_dataset(  plint ,
00089                        "Input" ,          /* label next to button   */
00090                        ANAT_ALL_MASK ,    /* take any anat datasets */
00091                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00092                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00093                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00094                     ) ;
00095 
00096    PLUTO_add_string( plint ,
00097                      "Output" ,  /* label next to textfield */
00098                      0,NULL ,    /* no fixed strings to choose among */
00099                      19          /* 19 spaces for typing in value */
00100                    ) ;
00101 
00102    /*---------- 2nd line --------*/
00103 
00104    PLUTO_add_option( plint ,
00105                      "Parameters" ,  /* label at left of input line */
00106                      "Parameters" ,  /* tag to return to plugin */
00107                      TRUE            /* is this mandatory? */
00108                    ) ;
00109 
00110    PLUTO_add_number( plint ,
00111                      "Base" ,    /* label next to chooser */
00112                      0 ,         /* smallest possible value */
00113                      98 ,        /* largest possible value */
00114                      0 ,         /* decimal shift (none in this case) */
00115                      3 ,         /* default value */
00116                      FALSE       /* allow user to edit value? */
00117                    ) ;
00118 
00119    /*---------- 3rd line --------*/
00120 
00121    PLUTO_add_option( plint ,
00122                      "Fine Fit" ,
00123                      "Fine Fit" ,
00124                      FALSE
00125                    ) ;
00126 
00127    PLUTO_add_number( plint , "Blur" , 0 , 40 , 1 , 10 , FALSE ) ;
00128    PLUTO_add_number( plint , "Dxy"  , 1 , 20 , 2 ,  7 , FALSE ) ;
00129    PLUTO_add_number( plint , "Dphi" , 1 , 50 , 2 , 21 , FALSE ) ;
00130 
00131    /*--------- done with interface setup ---------*/
00132 
00133    return plint ;
00134 }
00135 
00136 /***************************************************************************
00137   Main routine for this plugin (will be called from AFNI).
00138   If the return string is not NULL, some error transpired, and
00139   AFNI will popup the return string in a message box.
00140 ****************************************************************************/
00141 
00142 #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
00143 
00144 #define FREE_WORKSPACE                             \
00145   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
00146       FREEUP(bout) ; FREEUP(sout) ; FREEUP(fout) ; \
00147       FREEUP(dxar) ; FREEUP(dyar) ; FREEUP(phiar); \
00148     } while(0) ;
00149 
00150 char * IMREG_main( PLUGIN_interface * plint )
00151 {
00152    MCW_idcode * idc ;                          /* input dataset idcode */
00153    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00154    char * new_prefix , * str ;                 /* strings from user */
00155    int base , ntime , datum , nx,ny,nz , ii,kk , npix ;
00156    float                      dx,dy,dz ;
00157    MRI_IMARR * ims_in , * ims_out ;
00158    MRI_IMAGE * im , * imbase ;
00159 
00160    byte   ** bptr = NULL , ** bout = NULL ;
00161    short  ** sptr = NULL , ** sout = NULL ;
00162    float  ** fptr = NULL , ** fout = NULL ;
00163 
00164    float * dxar = NULL , * dyar = NULL , * phiar = NULL ;
00165 
00166    /*--------------------------------------------------------------------*/
00167    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00168 
00169    /*--------- go to first input line ---------*/
00170 
00171    PLUTO_next_option(plint) ;
00172 
00173    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00174    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00175    if( old_dset == NULL )
00176       return "*************************\n"
00177              "Cannot find Input Dataset\n"
00178              "*************************"  ;
00179 
00180    ntime = DSET_NUM_TIMES(old_dset) ;
00181    if( ntime < 2 )
00182       return "*****************************\n"
00183              "Dataset has only 1 time point\n"
00184              "*****************************"  ;
00185 
00186    ii = DSET_NVALS_PER_TIME(old_dset) ;
00187    if( ii > 1 )
00188       return "************************************\n"
00189              "Dataset has > 1 value per time point\n"
00190              "************************************"  ;
00191 
00192    nx = old_dset->daxes->nxx ; dx = old_dset->daxes->xxdel ;
00193    ny = old_dset->daxes->nyy ; dy = old_dset->daxes->yydel ; npix = nx*ny ;
00194    nz = old_dset->daxes->nzz ; dz = old_dset->daxes->zzdel ;
00195 
00196    if( nx != ny || fabs(dx) != fabs(dy) ){
00197 
00198 #ifdef IMREG_DEBUG
00199 fprintf(stderr,"\nIMREG: nx=%d ny=%d nz=%d  dx=%f dy=%f dz=%f\n",
00200         nx,ny,nz,dx,dy,dz ) ;
00201 #endif
00202 
00203       return "***********************************\n"
00204              "Dataset does not have square slices\n"
00205              "***********************************"  ;
00206    }
00207 
00208    new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00209    if( ! PLUTO_prefix_ok(new_prefix) )      /* check if it is OK */
00210       return "************************\n"
00211              "Output Prefix is illegal\n"
00212              "************************"  ;
00213 
00214    /*--------- go to next input line ---------*/
00215 
00216    PLUTO_next_option(plint) ;
00217 
00218    base = PLUTO_get_number(plint) ;
00219    if( base >= ntime )
00220       return "********************\n"
00221              "Base value too large\n"
00222              "********************"  ;
00223 
00224    /*--------- see if the 3rd option line is present --------*/
00225 
00226    str = PLUTO_get_optiontag( plint ) ;
00227    if( str != NULL ){
00228       float fsig , fdxy , fdph ;
00229       fsig = PLUTO_get_number(plint) * 0.42466090 ;
00230       fdxy = PLUTO_get_number(plint) ;
00231       fdph = PLUTO_get_number(plint) ;
00232       mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
00233 /* fprintf(stderr,"Set fine params = %f %f %f\n",fsig,fdxy,fdph) ; */
00234    }
00235 
00236    /*------------- ready to compute new dataset -----------*/
00237 
00238 #ifdef IMREG_DEBUG
00239 fprintf(stderr,"IMREG: loading dataset\n") ;
00240 #endif
00241 
00242    DSET_load( old_dset ) ;
00243 
00244    /*** 1) Copy the dataset in toto ***/
00245 
00246 #ifdef IMREG_DEBUG
00247 fprintf(stderr,"IMREG: Copying dataset\n") ;
00248 #endif
00249 
00250    new_dset = PLUTO_copy_dset( old_dset , new_prefix ) ;
00251    if( new_dset == NULL )
00252       return "****************************\n"
00253              "Failed to copy input dataset\n"
00254              "****************************"  ;
00255 
00256    /*** 2) Make an array of empty images ***/
00257 
00258 #ifdef IMREG_DEBUG
00259 fprintf(stderr,"IMREG: making empty images\n") ;
00260 #endif
00261 
00262    datum = DSET_BRICK_TYPE(new_dset,0) ;
00263 
00264    INIT_IMARR(ims_in) ;
00265    for( ii=0 ; ii < ntime ; ii++ ){
00266       im = mri_new_vol_empty( nx , ny , 1 , datum ) ;
00267       ADDTO_IMARR(ims_in,im) ;
00268    }
00269 
00270    imbase = mri_new_vol_empty( nx , ny , 1 , datum ) ;
00271 
00272    dxar  = (float *) malloc( sizeof(float) * ntime ) ;
00273    dyar  = (float *) malloc( sizeof(float) * ntime ) ;
00274    phiar = (float *) malloc( sizeof(float) * ntime ) ;
00275 
00276    /*** 3) Get pointers to sub-bricks in old and new datasets ***/
00277 
00278 #ifdef IMREG_DEBUG
00279 fprintf(stderr,"IMREG: getting input brick pointers\n") ;
00280 #endif
00281 
00282    switch( datum ){  /* pointer type depends on input datum type */
00283       case MRI_byte:
00284          bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
00285          bout = (byte **) malloc( sizeof(byte *) * ntime ) ;
00286          for( ii=0 ; ii < ntime ; ii++ ){
00287             bptr[ii] = (byte *) DSET_ARRAY(old_dset,ii) ;
00288             bout[ii] = (byte *) DSET_ARRAY(new_dset,ii) ;
00289          }
00290       break ;
00291 
00292       case MRI_short:
00293          sptr = (short **) malloc( sizeof(short *) * ntime ) ;
00294          sout = (short **) malloc( sizeof(short *) * ntime ) ;
00295          for( ii=0 ; ii < ntime ; ii++ ){
00296             sptr[ii] = (short *) DSET_ARRAY(old_dset,ii) ;
00297             sout[ii] = (short *) DSET_ARRAY(new_dset,ii) ;
00298          }
00299 
00300 #ifdef IMREG_DEBUG
00301 fprintf(stderr,"IMREG: sptr[0] = %p  sout[0] = %p\n",sptr[0],sout[0]) ;
00302 #endif
00303 
00304       break ;
00305 
00306       case MRI_float:
00307          fptr = (float **) malloc( sizeof(float *) * ntime ) ;
00308          fout = (float **) malloc( sizeof(float *) * ntime ) ;
00309          for( ii=0 ; ii < ntime ; ii++ ){
00310             fptr[ii] = (float *) DSET_ARRAY(old_dset,ii) ;
00311             fout[ii] = (float *) DSET_ARRAY(new_dset,ii) ;
00312          }
00313       break ;
00314    }
00315 
00316    /*** 4) Loop over slices ***/
00317 
00318    PLUTO_popup_meter(plint) ;
00319 
00320    for( kk=0 ; kk < nz ; kk++ ){
00321 
00322       /*** 4a) Setup ims_in images to point to input slices ***/
00323 
00324 #ifdef IMREG_DEBUG
00325 fprintf(stderr,"IMREG: slice %d -- setup input images\n",kk) ;
00326 #endif
00327 
00328       for( ii=0 ; ii < ntime ; ii++ ){
00329          im = IMARR_SUBIMAGE(ims_in,ii) ;
00330          switch( datum ){
00331             case MRI_byte:  mri_fix_data_pointer( bptr[ii] + kk*npix, im ) ; break ;
00332             case MRI_short: mri_fix_data_pointer( sptr[ii] + kk*npix, im ) ; break ;
00333             case MRI_float: mri_fix_data_pointer( fptr[ii] + kk*npix, im ) ; break ;
00334          }
00335       }
00336 
00337       /*** 4b) Setup im to point to base image ***/
00338 
00339 #ifdef IMREG_DEBUG
00340 fprintf(stderr,"IMREG: slice %d -- setup base image\n",kk) ;
00341 #endif
00342 
00343       switch( datum ){
00344          case MRI_byte:  mri_fix_data_pointer( bptr[base] + kk*npix, imbase ) ; break ;
00345          case MRI_short: mri_fix_data_pointer( sptr[base] + kk*npix, imbase ) ; break ;
00346          case MRI_float: mri_fix_data_pointer( fptr[base] + kk*npix, imbase ) ; break ;
00347       }
00348 
00349       /*** 4c) Register this slice at all times ***/
00350 
00351 #ifdef IMREG_DEBUG
00352 fprintf(stderr,"IMREG: slice %d -- register\n",kk) ;
00353 #endif
00354 
00355       ims_out = mri_align_dfspace( imbase , NULL , ims_in ,
00356                                    ALIGN_REGISTER_CODE , dxar,dyar,phiar ) ;
00357 
00358 if( ims_out == NULL )
00359    fprintf(stderr,"IMREG: mri_align_dfspace return NULL\n") ;
00360 
00361       /*** 4d) Put the output back in on top of the input;
00362                note that the output is always in MRI_float format ***/
00363 
00364 #ifdef IMREG_DEBUG
00365 fprintf(stderr,"IMREG: slice %d -- put output back into dataset\n",kk) ;
00366 #endif
00367 
00368       for( ii=0 ; ii < ntime ; ii++ ){
00369          switch( datum ){
00370            case MRI_byte:
00371               im = mri_to_mri( MRI_byte , IMARR_SUBIMAGE(ims_out,ii) ) ;
00372               memcpy( bout[ii] + kk*npix , MRI_BYTE_PTR(im) , sizeof(byte)*npix ) ;
00373               mri_free(im) ;
00374            break ;
00375 
00376            case MRI_short:
00377 #ifdef IMREG_DEBUG
00378 if( ii==0 )fprintf(stderr,"IMREG: conversion to short at ii=%d\n",ii) ;
00379 #endif
00380 
00381               im = mri_to_mri( MRI_short , IMARR_SUBIMAGE(ims_out,ii) ) ;
00382 
00383 #ifdef IMREG_DEBUG
00384 if( ii==0 )fprintf(stderr,"IMREG: copying to %p from %p\n",sout[ii] + kk*npix,MRI_SHORT_PTR(im)) ;
00385 #endif
00386 
00387               memcpy( sout[ii] + kk*npix , MRI_SHORT_PTR(im) , sizeof(short)*npix ) ;
00388 
00389 #ifdef IMREG_DEBUG
00390 if( ii==0 )fprintf(stderr,"IMREG: freeing\n") ;
00391 #endif
00392 
00393               mri_free(im) ;
00394            break ;
00395 
00396            case MRI_float:
00397               im = IMARR_SUBIMAGE(ims_out,ii) ;
00398               memcpy( fout[ii] + kk*npix , MRI_FLOAT_PTR(im) , sizeof(float)*npix ) ;
00399            break ;
00400          }
00401       }
00402 
00403       PLUTO_set_meter(plint, (100*(kk+1))/nz ) ;
00404 
00405       /*** 4e) Destroy the output images ***/
00406 
00407 #ifdef IMREG_DEBUG
00408 fprintf(stderr,"IMREG: destroying aligned output\n") ;
00409 #endif
00410 
00411       DESTROY_IMARR( ims_out ) ;
00412    }
00413 
00414    /*** 5) Destroy the empty images and other workspaces ***/
00415 
00416 #ifdef IMREG_DEBUG
00417 fprintf(stderr,"IMREG: destroy workspaces\n") ;
00418 #endif
00419 
00420    mri_clear_data_pointer(imbase) ; mri_free(imbase) ;
00421    for( ii=0 ; ii < ntime ; ii++ ){
00422       im = IMARR_SUBIMAGE(ims_in,ii) ;
00423       mri_clear_data_pointer(im) ;
00424    }
00425    DESTROY_IMARR(ims_in) ;
00426    FREE_WORKSPACE ;
00427 
00428    /*------------- let AFNI know about the new dataset ------------*/
00429 
00430 #ifdef IMREG_DEBUG
00431 fprintf(stderr,"IMREG: send result to AFNI\n") ;
00432 #endif
00433 
00434    PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00435 
00436    return NULL ;  /* null string returned means all was OK */
00437 }