00001
00002
00003
00004
00005
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
00016
00017
00018
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
00038
00039 char * IMREG_main( PLUGIN_interface * ) ;
00040
00041
00042
00043 static PLUGIN_interface * global_plint = NULL ;
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 DEFINE_PLUGIN_PROTOTYPE
00060
00061 PLUGIN_interface * PLUGIN_init( int ncall )
00062 {
00063 PLUGIN_interface * plint ;
00064
00065 if( ncall > 0 ) return NULL ;
00066
00067
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 ;
00077
00078 PLUTO_set_sequence( plint , "A:newdset:reg" ) ;
00079
00080
00081
00082 PLUTO_add_option( plint ,
00083 "Datasets" ,
00084 "DAtasets" ,
00085 TRUE
00086 ) ;
00087
00088 PLUTO_add_dataset( plint ,
00089 "Input" ,
00090 ANAT_ALL_MASK ,
00091 FUNC_FIM_MASK ,
00092 DIMEN_4D_MASK |
00093 BRICK_ALLREAL_MASK
00094 ) ;
00095
00096 PLUTO_add_string( plint ,
00097 "Output" ,
00098 0,NULL ,
00099 19
00100 ) ;
00101
00102
00103
00104 PLUTO_add_option( plint ,
00105 "Parameters" ,
00106 "Parameters" ,
00107 TRUE
00108 ) ;
00109
00110 PLUTO_add_number( plint ,
00111 "Base" ,
00112 0 ,
00113 98 ,
00114 0 ,
00115 3 ,
00116 FALSE
00117 ) ;
00118
00119
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
00132
00133 return plint ;
00134 }
00135
00136
00137
00138
00139
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 ;
00153 THD_3dim_dataset * old_dset , * new_dset ;
00154 char * new_prefix , * str ;
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
00168
00169
00170
00171 PLUTO_next_option(plint) ;
00172
00173 idc = PLUTO_get_idcode(plint) ;
00174 old_dset = PLUTO_find_dset(idc) ;
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) ;
00209 if( ! PLUTO_prefix_ok(new_prefix) )
00210 return "************************\n"
00211 "Output Prefix is illegal\n"
00212 "************************" ;
00213
00214
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
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
00234 }
00235
00236
00237
00238 #ifdef IMREG_DEBUG
00239 fprintf(stderr,"IMREG: loading dataset\n") ;
00240 #endif
00241
00242 DSET_load( old_dset ) ;
00243
00244
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
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
00277
00278 #ifdef IMREG_DEBUG
00279 fprintf(stderr,"IMREG: getting input brick pointers\n") ;
00280 #endif
00281
00282 switch( datum ){
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
00317
00318 PLUTO_popup_meter(plint) ;
00319
00320 for( kk=0 ; kk < nz ; kk++ ){
00321
00322
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
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
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
00362
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
00406
00407 #ifdef IMREG_DEBUG
00408 fprintf(stderr,"IMREG: destroying aligned output\n") ;
00409 #endif
00410
00411 DESTROY_IMARR( ims_out ) ;
00412 }
00413
00414
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
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 ;
00437 }