/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2000, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ #include "afni.h" #ifndef ALLOW_PLUGINS # error "Plugins not properly set up -- see machdep.h" #endif /*********************************************************************** Plugin to do what imreg.c does ************************************************************************/ /*--------------------- string to 'help' the user --------------------*/ static char helpstring[] = "Purpose: 2D (slice-wise) registration of 3D+time dataset\n" "[See also programs 'imreg' and '2dImReg']\n" "\n" "Input items to this plugin are:\n" " Datasets: Input = 3D+time dataset to process\n" " Output = Prefix for new dataset\n" " Parameters: Base = Time index for base image\n" "----\n" "This option is only for experimenting with the parameters\n" "of the final (fine) step of the registration algorithm:\n" " Fine Fit: Blur = FWHM of blurring prior to registration\n" " Dxy = Convergence tolerance for translations\n" " Dphi = Convergence tolerance for rotations\n" "Author -- RW Cox" ; /*----------------- prototypes for internal routines -----------------*/ char * IMREG_main( PLUGIN_interface * ) ; /* the entry point */ /*---------------------------- global data ---------------------------*/ static PLUGIN_interface * global_plint = NULL ; /*********************************************************************** Set up the interface to the user: 1) Create a new interface using "PLUTO_new_interface"; 2) For each line of inputs, create the line with "PLUTO_add_option" (this line of inputs can be optional or mandatory); 3) For each item on the line, create the item with "PLUTO_add_dataset" for a dataset chooser, "PLUTO_add_string" for a string chooser, "PLUTO_add_number" for a number chooser. ************************************************************************/ DEFINE_PLUGIN_PROTOTYPE PLUGIN_interface * PLUGIN_init( int ncall ) { PLUGIN_interface * plint ; /* will be the output of this routine */ if( ncall > 0 ) return NULL ; /* only one interface */ CHECK_IF_ALLOWED("2DREGISTRATION","2D Registration") ; /* 30 Sep 2016 */ /*---------------- set titles and call point ----------------*/ plint = PLUTO_new_interface( "2D Registration" , "2D Registration of a 3D+time Dataset" , helpstring , PLUGIN_CALL_VIA_MENU , IMREG_main ) ; PLUTO_add_hint( plint , "2D Registration of a 3D+time Dataset" ) ; global_plint = plint ; /* make global copy */ PLUTO_set_sequence( plint , "A:newdset:reg" ) ; /*--------- 1st line ---------*/ PLUTO_add_option( plint , "Datasets" , /* label at left of input line */ "DAtasets" , /* tag to return to plugin */ TRUE /* is this mandatory? */ ) ; PLUTO_add_dataset( plint , "Input" , /* label next to button */ ANAT_ALL_MASK , /* take any anat datasets */ FUNC_FIM_MASK , /* only allow fim funcs */ DIMEN_4D_MASK | /* need 3D+time datasets */ BRICK_ALLREAL_MASK /* need real-valued datasets */ ) ; PLUTO_add_string( plint , "Output" , /* label next to textfield */ 0,NULL , /* no fixed strings to choose among */ 19 /* 19 spaces for typing in value */ ) ; /*---------- 2nd line --------*/ PLUTO_add_option( plint , "Parameters" , /* label at left of input line */ "Parameters" , /* tag to return to plugin */ TRUE /* is this mandatory? */ ) ; PLUTO_add_number( plint , "Base" , /* label next to chooser */ 0 , /* smallest possible value */ 98 , /* largest possible value */ 0 , /* decimal shift (none in this case) */ 3 , /* default value */ FALSE /* allow user to edit value? */ ) ; /*---------- 3rd line --------*/ PLUTO_add_option( plint , "Fine Fit" , "Fine Fit" , FALSE ) ; PLUTO_add_number( plint , "Blur" , 0 , 40 , 1 , 10 , FALSE ) ; PLUTO_add_number( plint , "Dxy" , 1 , 20 , 2 , 7 , FALSE ) ; PLUTO_add_number( plint , "Dphi" , 1 , 50 , 2 , 21 , FALSE ) ; /*--------- done with interface setup ---------*/ return plint ; } /*************************************************************************** Main routine for this plugin (will be called from AFNI). If the return string is not NULL, some error transpired, and AFNI will popup the return string in a message box. ****************************************************************************/ #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0) #define FREE_WORKSPACE \ do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \ FREEUP(bout) ; FREEUP(sout) ; FREEUP(fout) ; \ FREEUP(dxar) ; FREEUP(dyar) ; FREEUP(phiar); \ } while(0) ; char * IMREG_main( PLUGIN_interface * plint ) { MCW_idcode * idc ; /* input dataset idcode */ THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */ char * new_prefix , * str ; /* strings from user */ int base , ntime , datum , nx,ny,nz , ii,kk , npix ; float dx,dy,dz ; MRI_IMARR * ims_in , * ims_out ; MRI_IMAGE * im , * imbase ; byte ** bptr = NULL , ** bout = NULL ; short ** sptr = NULL , ** sout = NULL ; float ** fptr = NULL , ** fout = NULL ; float * dxar = NULL , * dyar = NULL , * phiar = NULL ; /*--------------------------------------------------------------------*/ /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/ /*--------- go to first input line ---------*/ PLUTO_next_option(plint) ; idc = PLUTO_get_idcode(plint) ; /* get dataset item */ old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */ if( old_dset == NULL ) return "*************************\n" "Cannot find Input Dataset\n" "*************************" ; ntime = DSET_NUM_TIMES(old_dset) ; if( ntime < 2 ) return "*****************************\n" "Dataset has only 1 time point\n" "*****************************" ; ii = DSET_NVALS_PER_TIME(old_dset) ; if( ii > 1 ) return "************************************\n" "Dataset has > 1 value per time point\n" "************************************" ; nx = old_dset->daxes->nxx ; dx = old_dset->daxes->xxdel ; ny = old_dset->daxes->nyy ; dy = old_dset->daxes->yydel ; npix = nx*ny ; nz = old_dset->daxes->nzz ; dz = old_dset->daxes->zzdel ; if( nx != ny || fabs(dx) != fabs(dy) ){ #ifdef IMREG_DEBUG fprintf(stderr,"\nIMREG: nx=%d ny=%d nz=%d dx=%f dy=%f dz=%f\n", nx,ny,nz,dx,dy,dz ) ; #endif return "***********************************\n" "Dataset does not have square slices\n" "***********************************" ; } new_prefix = PLUTO_get_string(plint) ; /* get string item (the output prefix) */ if( ! PLUTO_prefix_ok(new_prefix) ) /* check if it is OK */ return "************************\n" "Output Prefix is illegal\n" "************************" ; /*--------- go to next input line ---------*/ PLUTO_next_option(plint) ; base = PLUTO_get_number(plint) ; if( base >= ntime ) return "********************\n" "Base value too large\n" "********************" ; /*--------- see if the 3rd option line is present --------*/ str = PLUTO_get_optiontag( plint ) ; if( str != NULL ){ float fsig , fdxy , fdph ; fsig = PLUTO_get_number(plint) * 0.42466090 ; fdxy = PLUTO_get_number(plint) ; fdph = PLUTO_get_number(plint) ; mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ; /* fprintf(stderr,"Set fine params = %f %f %f\n",fsig,fdxy,fdph) ; */ } /*------------- ready to compute new dataset -----------*/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: loading dataset\n") ; #endif DSET_load( old_dset ) ; /*** 1) Copy the dataset in toto ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: Copying dataset\n") ; #endif new_dset = PLUTO_copy_dset( old_dset , new_prefix ) ; if( new_dset == NULL ) return "****************************\n" "Failed to copy input dataset\n" "****************************" ; /*** 2) Make an array of empty images ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: making empty images\n") ; #endif datum = DSET_BRICK_TYPE(new_dset,0) ; INIT_IMARR(ims_in) ; for( ii=0 ; ii < ntime ; ii++ ){ im = mri_new_vol_empty( nx , ny , 1 , datum ) ; ADDTO_IMARR(ims_in,im) ; } imbase = mri_new_vol_empty( nx , ny , 1 , datum ) ; dxar = (float *) malloc( sizeof(float) * ntime ) ; dyar = (float *) malloc( sizeof(float) * ntime ) ; phiar = (float *) malloc( sizeof(float) * ntime ) ; /*** 3) Get pointers to sub-bricks in old and new datasets ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: getting input brick pointers\n") ; #endif switch( datum ){ /* pointer type depends on input datum type */ case MRI_byte: bptr = (byte **) malloc( sizeof(byte *) * ntime ) ; bout = (byte **) malloc( sizeof(byte *) * ntime ) ; for( ii=0 ; ii < ntime ; ii++ ){ bptr[ii] = (byte *) DSET_ARRAY(old_dset,ii) ; bout[ii] = (byte *) DSET_ARRAY(new_dset,ii) ; } break ; case MRI_short: sptr = (short **) malloc( sizeof(short *) * ntime ) ; sout = (short **) malloc( sizeof(short *) * ntime ) ; for( ii=0 ; ii < ntime ; ii++ ){ sptr[ii] = (short *) DSET_ARRAY(old_dset,ii) ; sout[ii] = (short *) DSET_ARRAY(new_dset,ii) ; } #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: sptr[0] = %p sout[0] = %p\n",sptr[0],sout[0]) ; #endif break ; case MRI_float: fptr = (float **) malloc( sizeof(float *) * ntime ) ; fout = (float **) malloc( sizeof(float *) * ntime ) ; for( ii=0 ; ii < ntime ; ii++ ){ fptr[ii] = (float *) DSET_ARRAY(old_dset,ii) ; fout[ii] = (float *) DSET_ARRAY(new_dset,ii) ; } break ; } /*** 4) Loop over slices ***/ PLUTO_popup_meter(plint) ; for( kk=0 ; kk < nz ; kk++ ){ /*** 4a) Setup ims_in images to point to input slices ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: slice %d -- setup input images\n",kk) ; #endif for( ii=0 ; ii < ntime ; ii++ ){ im = IMARR_SUBIMAGE(ims_in,ii) ; switch( datum ){ case MRI_byte: mri_fix_data_pointer( bptr[ii] + kk*npix, im ) ; break ; case MRI_short: mri_fix_data_pointer( sptr[ii] + kk*npix, im ) ; break ; case MRI_float: mri_fix_data_pointer( fptr[ii] + kk*npix, im ) ; break ; } } /*** 4b) Setup im to point to base image ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: slice %d -- setup base image\n",kk) ; #endif switch( datum ){ case MRI_byte: mri_fix_data_pointer( bptr[base] + kk*npix, imbase ) ; break ; case MRI_short: mri_fix_data_pointer( sptr[base] + kk*npix, imbase ) ; break ; case MRI_float: mri_fix_data_pointer( fptr[base] + kk*npix, imbase ) ; break ; } /*** 4c) Register this slice at all times ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: slice %d -- register\n",kk) ; #endif ims_out = mri_align_dfspace( imbase , NULL , ims_in , ALIGN_REGISTER_CODE , dxar,dyar,phiar ) ; if( ims_out == NULL ) fprintf(stderr,"IMREG: mri_align_dfspace return NULL\n") ; /*** 4d) Put the output back in on top of the input; note that the output is always in MRI_float format ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: slice %d -- put output back into dataset\n",kk) ; #endif for( ii=0 ; ii < ntime ; ii++ ){ switch( datum ){ case MRI_byte: im = mri_to_mri( MRI_byte , IMARR_SUBIMAGE(ims_out,ii) ) ; memcpy( bout[ii] + kk*npix , MRI_BYTE_PTR(im) , sizeof(byte)*npix ) ; mri_free(im) ; break ; case MRI_short: #ifdef IMREG_DEBUG if( ii==0 )fprintf(stderr,"IMREG: conversion to short at ii=%d\n",ii) ; #endif im = mri_to_mri( MRI_short , IMARR_SUBIMAGE(ims_out,ii) ) ; #ifdef IMREG_DEBUG if( ii==0 )fprintf(stderr,"IMREG: copying to %p from %p\n",sout[ii] + kk*npix,MRI_SHORT_PTR(im)) ; #endif memcpy( sout[ii] + kk*npix , MRI_SHORT_PTR(im) , sizeof(short)*npix ) ; #ifdef IMREG_DEBUG if( ii==0 )fprintf(stderr,"IMREG: freeing\n") ; #endif mri_free(im) ; break ; case MRI_float: im = IMARR_SUBIMAGE(ims_out,ii) ; memcpy( fout[ii] + kk*npix , MRI_FLOAT_PTR(im) , sizeof(float)*npix ) ; break ; } } PLUTO_set_meter(plint, (100*(kk+1))/nz ) ; /*** 4e) Destroy the output images ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: destroying aligned output\n") ; #endif DESTROY_IMARR( ims_out ) ; } /*** 5) Destroy the empty images and other workspaces ***/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: destroy workspaces\n") ; #endif mri_clear_data_pointer(imbase) ; mri_free(imbase) ; for( ii=0 ; ii < ntime ; ii++ ){ im = IMARR_SUBIMAGE(ims_in,ii) ; mri_clear_data_pointer(im) ; } DESTROY_IMARR(ims_in) ; FREE_WORKSPACE ; /*------------- let AFNI know about the new dataset ------------*/ #ifdef IMREG_DEBUG fprintf(stderr,"IMREG: send result to AFNI\n") ; #endif PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ; return NULL ; /* null string returned means all was OK */ }