00001
00002
00003
00004
00005
00006 #include "mrilib.h"
00007
00008
00009
00010 #define MAXPAR 99
00011
00012 typedef struct { int np ; float vp ; } fixed_param ;
00013
00014 static int nparfix ;
00015 static fixed_param parfix[MAXPAR] ;
00016
00017 static float parvec[MAXPAR] ;
00018
00019 static void (*warp_parset)(void) = NULL ;
00020
00021 void load_parvec( int np, float *pv ){
00022 memcpy( parvec , pv , sizeof(float)*np ) ;
00023 if( warp_parset != NULL ) warp_parset() ;
00024 }
00025
00026
00027
00028 static void (*warp_for)( float,float,float , float *,float *,float *) ;
00029 static void (*warp_inv)( float,float,float , float *,float *,float *) ;
00030
00031 static THD_vecmat ijk_to_xyz , xyz_to_ijk ;
00032
00033
00034
00035 void ijk_warp_for( float ii , float jj , float kk ,
00036 float *pp , float *qq , float *rr )
00037 {
00038 THD_fvec3 xxx , yyy ;
00039
00040 LOAD_FVEC3( xxx , ii,jj,kk ) ;
00041 yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
00042 warp_for( yyy.xyz[0] , yyy.xyz[1] , yyy.xyz[2] ,
00043 &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
00044 yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
00045 *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
00046 }
00047
00048
00049
00050 void ijk_warp_inv( float ii , float jj , float kk ,
00051 float *pp , float *qq , float *rr )
00052 {
00053 THD_fvec3 xxx , yyy ;
00054
00055 LOAD_FVEC3( xxx , ii,jj,kk ) ;
00056 yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
00057 warp_inv( yyy.xyz[0] , yyy.xyz[1] , yyy.xyz[2] ,
00058 &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
00059 yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
00060 *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
00061 }
00062
00063
00064
00065 #define WARPDRIVE_SHIFT 1
00066 #define WARPDRIVE_ROTATE 2
00067 #define WARPDRIVE_SCALE 3
00068 #define WARPDRIVE_AFFINE 4
00069 #define WARPDRIVE_BILINEAR 5
00070
00071 #define WARPDRIVE_IS_AFFINE(wc) \
00072 ( (wc) >= WARPDRIVE_SHIFT && (wc) <= WARPDRIVE_AFFINE )
00073
00074
00075
00076
00077 static float xsh , ysh , zsh ;
00078
00079 void parset_shift(void)
00080 {
00081 xsh = parvec[0] ; ysh = parvec[1] ; zsh = parvec[2] ;
00082 }
00083
00084 void warper_shift_for( float aa , float bb , float cc ,
00085 float *p , float *q , float *r )
00086 {
00087 *p = aa+xsh ; *q = bb+ysh ; *r = cc+zsh ;
00088 }
00089
00090 void warper_shift_inv( float aa , float bb , float cc ,
00091 float *p , float *q , float *r )
00092 {
00093 *p = aa-xsh ; *q = bb-ysh ; *r = cc-zsh ;
00094 }
00095
00096
00097
00098
00099 static THD_vecmat mv_for , mv_inv ;
00100
00101 void warper_affine_for( float aa , float bb , float cc ,
00102 float *p , float *q , float *r )
00103 {
00104 THD_fvec3 v , w ;
00105 LOAD_FVEC3( v , aa,bb,cc ) ;
00106 w = VECMAT_VEC( mv_for , v ) ;
00107 *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
00108 }
00109
00110 void warper_affine_inv( float aa , float bb , float cc ,
00111 float *p , float *q , float *r )
00112 {
00113 THD_fvec3 v , w ;
00114 LOAD_FVEC3( v , aa,bb,cc ) ;
00115 w = VECMAT_VEC( mv_inv , v ) ;
00116 *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
00117 }
00118
00119
00120
00121
00122
00123
00124 static THD_mat33 rot_matrix( int ax1, double th1,
00125 int ax2, double th2, int ax3, double th3 )
00126 {
00127 THD_mat33 q , p ;
00128 LOAD_ROT_MAT( q , th1 , ax1 ) ;
00129 LOAD_ROT_MAT( p , th2 , ax2 ) ; q = MAT_MUL( p , q ) ;
00130 LOAD_ROT_MAT( p , th3 , ax3 ) ; q = MAT_MUL( p , q ) ;
00131 return q ;
00132 }
00133
00134 #define D2R (PI/180.0)
00135
00136 #define MATORDER_SDU 1
00137 #define MATORDER_SUD 2
00138 #define MATORDER_DSU 3
00139 #define MATORDER_DUS 4
00140 #define MATORDER_USD 5
00141 #define MATORDER_UDS 6
00142
00143 static int matorder = MATORDER_SDU ;
00144 static int dcode = DELTA_AFTER ;
00145
00146 #define SMAT_UPPER 1
00147 #define SMAT_LOWER 2
00148
00149 static int smat = SMAT_LOWER ;
00150
00151 void parset_affine(void)
00152 {
00153 THD_mat33 ss,dd,uu,aa,bb ;
00154 THD_fvec3 vv ;
00155
00156 #if 0
00157 { int ii;fprintf(stderr,"\nparset:");
00158 for(ii=0;ii<12;ii++)fprintf(stderr," %g",parvec[ii]); fprintf(stderr,"\n");}
00159 #endif
00160
00161
00162
00163 uu = rot_matrix( 2, D2R*parvec[3] , 0, D2R*parvec[4] , 1, D2R*parvec[5] ) ;
00164
00165
00166
00167 LOAD_DIAG_MAT( dd , parvec[6] , parvec[7] , parvec[8] ) ;
00168
00169
00170
00171 switch( smat ){
00172 default:
00173 case SMAT_LOWER:
00174 LOAD_MAT( ss , 1.0 , 0.0 , 0.0 ,
00175 parvec[9] , 1.0 , 0.0 ,
00176 parvec[10] , parvec[11] , 1.0 ) ;
00177 break ;
00178
00179 case SMAT_UPPER:
00180 LOAD_MAT( ss , 1.0 , parvec[9] , parvec[10] ,
00181 0.0 , 1.0 , parvec[11] ,
00182 0.0 , 0.0 , 1.0 ) ;
00183 break ;
00184 }
00185
00186
00187
00188 switch( matorder ){
00189 case MATORDER_SDU: aa = MAT_MUL(ss,dd) ; bb = uu ; break ;
00190 case MATORDER_SUD: aa = MAT_MUL(ss,uu) ; bb = dd ; break ;
00191 case MATORDER_DSU: aa = MAT_MUL(dd,ss) ; bb = uu ; break ;
00192 case MATORDER_DUS: aa = MAT_MUL(dd,uu) ; bb = ss ; break ;
00193 case MATORDER_USD: aa = MAT_MUL(uu,ss) ; bb = dd ; break ;
00194 case MATORDER_UDS: aa = MAT_MUL(uu,dd) ; bb = ss ; break ;
00195 }
00196 mv_for.mm = MAT_MUL(aa,bb) ;
00197
00198 LOAD_FVEC3( vv , parvec[0] , parvec[1] , parvec[2] ) ;
00199
00200 switch( dcode ){
00201 case DELTA_AFTER: mv_for.vv = vv ; break ;
00202 case DELTA_BEFORE: mv_for.vv = MATVEC( mv_for.mm , vv ) ; break ;
00203 }
00204
00205 mv_inv = INV_VECMAT( mv_for ) ;
00206
00207 #if 0
00208 DUMP_VECMAT("mv_for",mv_for) ; DUMP_VECMAT("mv_inv",mv_inv) ;
00209 #endif
00210 }
00211
00212
00213
00214
00215 static float dd_fac = 1.0f ;
00216 static float dd_for[3][3][3] , dd_inv[3][3][3] ;
00217
00218 void warper_bilinear_for( float aa , float bb , float cc ,
00219 float *p , float *q , float *r )
00220 {
00221 THD_fvec3 v,w ; THD_mat33 dd,ee ;
00222
00223 LOAD_FVEC3( v , aa,bb,cc ) ;
00224 w = VECMAT_VEC( mv_for , v ) ;
00225
00226 dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
00227 dd.mat[0][1] = dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
00228 dd.mat[0][2] = dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
00229 dd.mat[1][0] = dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
00230 dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
00231 dd.mat[1][2] = dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
00232 dd.mat[2][0] = dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
00233 dd.mat[2][1] = dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
00234 dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
00235
00236 ee = MAT_INV(dd) ;
00237 v = MATVEC(ee,w) ;
00238
00239 *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
00240 }
00241
00242 void warper_bilinear_inv( float aa , float bb , float cc ,
00243 float *p , float *q , float *r )
00244 {
00245 THD_fvec3 v,w ; THD_mat33 dd,ee ;
00246
00247 LOAD_FVEC3( v , aa,bb,cc ) ;
00248 w = VECMAT_VEC( mv_inv , v ) ;
00249
00250 dd.mat[0][0] = 1.0f + dd_inv[0][0][0]*aa + dd_inv[0][0][1]*bb + dd_inv[0][0][2]*cc ;
00251 dd.mat[0][1] = dd_inv[0][1][0]*aa + dd_inv[0][1][1]*bb + dd_inv[0][1][2]*cc ;
00252 dd.mat[0][2] = dd_inv[0][2][0]*aa + dd_inv[0][2][1]*bb + dd_inv[0][2][2]*cc ;
00253 dd.mat[1][0] = dd_inv[1][0][0]*aa + dd_inv[1][0][1]*bb + dd_inv[1][0][2]*cc ;
00254 dd.mat[1][1] = 1.0f + dd_inv[1][1][0]*aa + dd_inv[1][1][1]*bb + dd_inv[1][1][2]*cc ;
00255 dd.mat[1][2] = dd_inv[1][2][0]*aa + dd_inv[1][2][1]*bb + dd_inv[1][2][2]*cc ;
00256 dd.mat[2][0] = dd_inv[2][0][0]*aa + dd_inv[2][0][1]*bb + dd_inv[2][0][2]*cc ;
00257 dd.mat[2][1] = dd_inv[2][1][0]*aa + dd_inv[2][1][1]*bb + dd_inv[2][1][2]*cc ;
00258 dd.mat[2][2] = 1.0f + dd_inv[2][2][0]*aa + dd_inv[2][2][1]*bb + dd_inv[2][2][2]*cc ;
00259
00260 ee = MAT_INV(dd) ;
00261 v = MATVEC(ee,w) ;
00262
00263 *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
00264 }
00265
00266 float warper_bilinear_det( float ii , float jj , float kk )
00267 {
00268 THD_fvec3 x,v,w ; THD_mat33 dd,ee ; int i,j ; float edet,adet,ddet , aa,bb,cc ;
00269 static int first=1 ;
00270
00271 LOAD_FVEC3( x , ii,jj,kk ) ; v = VECMAT_VEC( ijk_to_xyz , x ) ;
00272 UNLOAD_FVEC3( v , aa,bb,cc ) ; w = VECMAT_VEC( mv_for , v ) ;
00273
00274 dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
00275 dd.mat[0][1] = dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
00276 dd.mat[0][2] = dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
00277 dd.mat[1][0] = dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
00278 dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
00279 dd.mat[1][2] = dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
00280 dd.mat[2][0] = dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
00281 dd.mat[2][1] = dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
00282 dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
00283
00284 ddet = MAT_DET(dd) ;
00285 if( first && fabs(ddet) < 0.001f ){
00286 fprintf(stderr,"******* ddet=%g ii,jj,kk=%g %g %g aa,bb,cc=%g %g %g\n",
00287 ddet,ii,jj,kk, aa,bb,cc ) ;
00288 DUMP_MAT33("dd",dd) ;
00289 DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
00290 DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
00291 first = 0 ;
00292 }
00293
00294 ee = MAT_INV(dd) ; edet = MAT_DET(ee) ; v = MATVEC(ee,w) ;
00295
00296 for( i=0 ; i < 3 ; i++ )
00297 for( j=0 ; j < 3 ; j++ )
00298 dd.mat[i][j] = mv_for.mm.mat[i][j] - dd_for[i][0][j]*v.xyz[0]
00299 - dd_for[i][1][j]*v.xyz[1]
00300 - dd_for[i][2][j]*v.xyz[2] ;
00301
00302 adet = MAT_DET(dd) ;
00303 return (adet*edet) ;
00304 }
00305
00306 void parset_bilinear(void)
00307 {
00308 THD_mat33 ai ; THD_fvec3 df,di ; int i,j,k ;
00309
00310 parset_affine() ;
00311
00312
00313
00314 dd_for[0][0][0] = parvec[12]; dd_for[0][0][1] = parvec[13]; dd_for[0][0][2] = parvec[14];
00315 dd_for[0][1][0] = parvec[15]; dd_for[0][1][1] = parvec[16]; dd_for[0][1][2] = parvec[17];
00316 dd_for[0][2][0] = parvec[18]; dd_for[0][2][1] = parvec[19]; dd_for[0][2][2] = parvec[20];
00317 dd_for[1][0][0] = parvec[21]; dd_for[1][0][1] = parvec[22]; dd_for[1][0][2] = parvec[23];
00318 dd_for[1][1][0] = parvec[24]; dd_for[1][1][1] = parvec[25]; dd_for[1][1][2] = parvec[26];
00319 dd_for[1][2][0] = parvec[27]; dd_for[1][2][1] = parvec[28]; dd_for[1][2][2] = parvec[29];
00320 dd_for[2][0][0] = parvec[30]; dd_for[2][0][1] = parvec[31]; dd_for[2][0][2] = parvec[32];
00321 dd_for[2][1][0] = parvec[33]; dd_for[2][1][1] = parvec[34]; dd_for[2][1][2] = parvec[35];
00322 dd_for[2][2][0] = parvec[36]; dd_for[2][2][1] = parvec[37]; dd_for[2][2][2] = parvec[38];
00323
00324 for( i=0 ; i < 3 ; i++ )
00325 for( j=0 ; j < 3 ; j++ )
00326 for( k=0 ; k < 3 ; k++ ) dd_for[i][j][k] *= dd_fac ;
00327
00328
00329
00330 ai = mv_inv.mm ;
00331 for( k=0 ; k < 3 ; k++ ){
00332 for( j=0 ; j < 3 ; j++ ){
00333 LOAD_FVEC3( df , -dd_for[0][k][j] , -dd_for[1][k][j] , -dd_for[2][k][j] ) ;
00334 di = MATVEC( ai , df ) ;
00335 UNLOAD_FVEC3( di , dd_inv[0][j][k] , dd_inv[1][j][k] , dd_inv[2][j][k] ) ;
00336 }
00337 }
00338
00339 #if 0
00340 fprintf(stderr,"++++++++++++ parset_bilinear: dd_fac = %g\n",dd_fac) ;
00341 for( k=0 ; k < 3 ; k++ ){
00342 fprintf(stderr," +-------+ dd_for[][][%d] | dd_inv[][][%d]\n"
00343 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n"
00344 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n"
00345 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n" ,
00346 k , k ,
00347 dd_for[0][0][k] , dd_for[0][1][k] , dd_for[0][2][k] ,
00348 dd_inv[0][0][k] , dd_inv[0][1][k] , dd_inv[0][2][k] ,
00349 dd_for[1][0][k] , dd_for[1][1][k] , dd_for[1][2][k] ,
00350 dd_inv[1][0][k] , dd_inv[1][1][k] , dd_inv[1][2][k] ,
00351 dd_for[2][0][k] , dd_for[2][1][k] , dd_for[2][2][k] ,
00352 dd_inv[2][0][k] , dd_inv[2][1][k] , dd_inv[2][2][k] ) ;
00353 }
00354 #endif
00355 }
00356
00357
00358
00359 int main( int argc , char * argv[] )
00360 {
00361 THD_3dim_dataset *inset=NULL, *outset=NULL, *baset=NULL, *wtset=NULL ;
00362 MRI_warp3D_align_basis abas ;
00363 char *prefix="warpdriven" ;
00364 int warpdrive_code=-1 , nerr , nx,ny,nz , nopt=1 ;
00365 float dx,dy,dz , vp ;
00366 int kim , nvals , kpar , np , nfree ;
00367 MRI_IMAGE *qim , *tim , *fim ;
00368 float clip_baset=0.0f , clip_inset=0.0f ;
00369 char *W_1Dfile=NULL ;
00370 float **parsave=NULL ;
00371 int output_float=0 ;
00372 char *base_idc=NULL , *wt_idc=NULL ;
00373 int ctstart = NI_clock_time() ;
00374
00375
00376
00377 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00378 printf("\n"
00379 "Usage: 3dWarpDrive [options] dataset\n"
00380 "Warp a dataset to match another one (the base).\n"
00381 "\n"
00382 "This program is a generalization of 3dvolreg. It tries to find\n"
00383 "a spatial transformation that warps a given dataset to match an\n"
00384 "input dataset (given by the -base option). It will be slow.\n"
00385 "\n"
00386 "--------------------------\n"
00387 "Transform Defining Options: [exactly one of these must be used]\n"
00388 "--------------------------\n"
00389 " -shift_only = 3 parameters (shifts)\n"
00390 " -shift_rotate = 6 parameters (shifts + angles)\n"
00391 " -shift_rotate_scale = 9 parameters (shifts + angles + scale factors)\n"
00392 " -affine_general = 12 parameters (3 shifts + 3x3 matrix)\n"
00393 " -bilinear_general = 39 parameters (3 + 3x3 + 3x3x3)\n"
00394 "\n"
00395 " N.B.: At this time, the image intensity is NOT \n"
00396 " adjusted for the Jacobian of the transformation.\n"
00397 " N.B.: -bilinear_general is not yet implemented.\n"
00398 "\n"
00399 "-------------\n"
00400 "Other Options:\n"
00401 "-------------\n"
00402 " -linear }\n"
00403 " -cubic } = Chooses spatial interpolation method.\n"
00404 " -NN } = [default = linear; inaccurate but fast]\n"
00405 " -quintic } [for accuracy, try '-cubic -final quintic']\n"
00406 "\n"
00407 " -base bbb = Load dataset 'bbb' as the base to which the\n"
00408 " input dataset will be matched.\n"
00409 " [This is a mandatory option]\n"
00410 "\n"
00411 " -verb = Print out lots of information along the way.\n"
00412 " -prefix ppp = Sets the prefix of the output dataset.\n"
00413 " -input ddd = You can put the input dataset anywhere in the\n"
00414 " command line option list by using the '-input'\n"
00415 " option, instead of always putting it last.\n"
00416 "\n"
00417 "-----------------\n"
00418 "Technical Options:\n"
00419 "-----------------\n"
00420 " -maxite m = Allow up to 'm' iterations for convergence.\n"
00421 " -delta d = Distance, in voxel size, used to compute\n"
00422 " image derivatives using finite differences.\n"
00423 " [Default=1.0]\n"
00424 " -weight wset = Set the weighting applied to each voxel\n"
00425 " proportional to the brick specified here.\n"
00426 " [Default=computed by program from base]\n"
00427 " -thresh t = Set the convergence parameter to be RMS 't' voxels\n"
00428 " movement between iterations. [Default=0.03]\n"
00429 " -twopass = Do the parameter estimation in two passes,\n"
00430 " coarse-but-fast first, then fine-but-slow second\n"
00431 " (much like the same option in program 3dvolreg).\n"
00432 " This is useful if large-ish warping is needed to\n"
00433 " align the volumes.\n"
00434 " -final 'mode' = Set the final warp to be interpolated using 'mode'\n"
00435 " instead of the spatial interpolation method used\n"
00436 " to find the warp parameters.\n"
00437 " -parfix n v = Fix the n'th parameter of the warp model to\n"
00438 " the value 'v'. More than one -parfix option\n"
00439 " can be used, to fix multiple parameters.\n"
00440 " -1Dfile ename = Write out the warping parameters to the file\n"
00441 " named 'ename'. Each sub-brick of the input\n"
00442 " dataset gets one line in this file. Each\n"
00443 " parameter in the model gets one column.\n"
00444 " -float = Write output dataset in float format, even if\n"
00445 " input dataset is short or byte.\n"
00446 "\n"
00447 "----------------------\n"
00448 "AFFINE TRANSFORMATIONS:\n"
00449 "----------------------\n"
00450 "The options below control how the affine tranformations\n"
00451 "(-shift_rotate, -shift_rotate_scale, -affine_general)\n"
00452 "are structured in terms of 3x3 matrices:\n"
00453 "\n"
00454 " -SDU or -SUD }= Set the order of the matrix multiplication\n"
00455 " -DSU or -DUS }= for the affine transformations:\n"
00456 " -USD or -UDS }= S = triangular shear (params #10-12)\n"
00457 " D = diagonal scaling matrix (params #7-9)\n"
00458 " U = rotation matrix (params #4-6)\n"
00459 " Default order is '-SDU', which means that\n"
00460 " the U matrix is applied first, then the\n"
00461 " D matrix, then the S matrix.\n"
00462 "\n"
00463 " -Supper }= Set the S matrix to be upper or lower\n"
00464 " -Slower }= triangular [Default=lower triangular]\n"
00465 "\n"
00466 " -ashift OR }= Apply the shift parameters (#1-3) after OR\n"
00467 " -bshift }= before the matrix transformation. [Default=after]\n"
00468 "\n"
00469 "The matrices are specified in DICOM-ordered (x=-R+L,y=-A+P,z=-I+S)\n"
00470 "coordinates as:\n"
00471 "\n"
00472 " [U] = [Rotate_y(param#6)] [Rotate_x(param#5)] [Rotate_z(param #4)]\n"
00473 " (angles are in degrees)\n"
00474 "\n"
00475 " [D] = diag( param#7 , param#8 , param#9 )\n"
00476 "\n"
00477 " [ 1 0 0 ] [ 1 param#10 param#11 ]\n"
00478 " [S] = [ param#10 1 0 ] OR [ 0 1 param#12 ]\n"
00479 " [ param#11 param#12 1 ] [ 0 0 1 ]\n"
00480 "\n"
00481 " For example, the default (-SDU/-ashift/-Slower) has the warp\n"
00482 " specified as [x]_warped = [S] [D] [U] [x]_in + [shift].\n"
00483 " The shift vector comprises parameters #1, #2, and #3.\n"
00484 "\n"
00485 " The goal of the program is to find the warp parameters such that\n"
00486 " I([x]_warped) = s * J([x]_in)\n"
00487 " as closely as possible in a weighted least squares sense, where\n"
00488 " 's' is a scaling factor (an extra, invisible, parameter), J(x)\n"
00489 " is the base image, I(x) is the input image, and the weight image\n"
00490 " is a blurred copy of J(x).\n"
00491 "\n"
00492 " Using '-parfix', you can specify that some of these parameters\n"
00493 " are fixed. For example, '-shift_rotate_scale' is equivalent\n"
00494 " '-affine_general -parfix 10 0 -parfix 11 0 -parfix 12 0'.\n"
00495 " Don't attempt to use the '-parfix' option unless you understand\n"
00496 " this example!\n"
00497 "\n"
00498 "-------------------------\n"
00499 " RWCox - November 2004\n"
00500 "-------------------------\n"
00501 ) ;
00502 exit(0) ;
00503 }
00504
00505
00506
00507 mainENTRY("3dWarpDrive main"); machdep(); AFNI_logger("3dWarpDrive",argc,argv);
00508 PRINT_VERSION("3dWarpDrive") ;
00509
00510
00511
00512 abas.nparam = 0 ;
00513 abas.param = NULL ;
00514 abas.scale_init = 1.0f ;
00515 abas.delfac = 1.0f ;
00516 abas.tolfac = 0.03f ;
00517 abas.twoblur = 0.0f ;
00518 abas.regmode = MRI_LINEAR ;
00519 abas.regfinal = -1 ;
00520 abas.verb = 0 ;
00521 abas.max_iter = 0 ;
00522 abas.wtproc = 1 ;
00523 abas.imbase = NULL ;
00524 abas.imwt = NULL ;
00525 abas.vwfor = ijk_warp_for ;
00526 abas.vwinv = ijk_warp_inv ;
00527 abas.vwset = load_parvec ;
00528 abas.vwdet = NULL ;
00529
00530 abas.xedge = abas.yedge = abas.zedge = -1 ;
00531 abas.imww = abas.imap = abas.imps = abas.imsk = NULL ;
00532 abas.imps_blur = NULL ;
00533
00534 nparfix = 0 ;
00535
00536
00537
00538 while( nopt < argc && argv[nopt][0] == '-' ){
00539
00540
00541
00542 if( strcmp(argv[nopt],"-float") == 0 ){
00543 output_float = 1 ; nopt++ ; continue ;
00544 }
00545
00546
00547
00548 if( strcmp(argv[nopt],"-1Dfile") == 0 ){
00549 if( ++nopt >= argc )
00550 ERROR_exit("Need 1 parameter afer -1Dfile!\n");
00551 W_1Dfile = strdup( argv[nopt] ) ;
00552 if( !THD_filename_ok(W_1Dfile) )
00553 ERROR_exit("Name after -1Dfile has bad characters!\n") ;
00554 nopt++ ; continue ;
00555 }
00556
00557
00558
00559 if( strcmp(argv[nopt],"-twopass") == 0 ){
00560 float bbb = AFNI_numenv("AFNI_WARPDRIVE_TWOBLUR") ;
00561 abas.twoblur = (bbb==0.0f) ? 3.0f : bbb ;
00562 nopt++ ; continue ;
00563 }
00564
00565
00566
00567 if( strcmp(argv[nopt],"-SDU") == 0 ){
00568 matorder = MATORDER_SDU ; nopt++ ; continue ;
00569 }
00570 if( strcmp(argv[nopt],"-SUD") == 0 ){
00571 matorder = MATORDER_SUD ; nopt++ ; continue ;
00572 }
00573 if( strcmp(argv[nopt],"-DSU") == 0 ){
00574 matorder = MATORDER_DSU ; nopt++ ; continue ;
00575 }
00576 if( strcmp(argv[nopt],"-DUS") == 0 ){
00577 matorder = MATORDER_DUS ; nopt++ ; continue ;
00578 }
00579 if( strcmp(argv[nopt],"-USD") == 0 ){
00580 matorder = MATORDER_USD ; nopt++ ; continue ;
00581 }
00582 if( strcmp(argv[nopt],"-UDS") == 0 ){
00583 matorder = MATORDER_UDS ; nopt++ ; continue ;
00584 }
00585 if( strcmp(argv[nopt],"-ashift") == 0 ){
00586 dcode = DELTA_AFTER ; nopt++ ; continue ;
00587 }
00588 if( strcmp(argv[nopt],"-bshift") == 0 ){
00589 dcode = DELTA_BEFORE ; nopt++ ; continue ;
00590 }
00591 if( strcmp(argv[nopt],"-Slower") == 0 ){
00592 smat = SMAT_LOWER ; nopt++ ; continue ;
00593 }
00594 if( strcmp(argv[nopt],"-Supper") == 0 ){
00595 smat = SMAT_UPPER ; nopt++ ; continue ;
00596 }
00597
00598
00599
00600 if( strcmp(argv[nopt],"-final") == 0 ){
00601 char *str ;
00602
00603 if( ++nopt >= argc )
00604 ERROR_exit("Need 1 parameter afer -final!\n") ;
00605 str = argv[nopt] ; if( *str == '-' ) str++ ;
00606
00607 if( strcmp(str,"cubic") == 0 ) abas.regfinal = MRI_CUBIC ;
00608 else if( strcmp(str,"quintic") == 0 ) abas.regfinal = MRI_QUINTIC ;
00609 else if( strcmp(str,"linear") == 0 ) abas.regfinal = MRI_LINEAR ;
00610 else if( strcmp(str,"NN") == 0 ) abas.regfinal = MRI_NN ;
00611 else
00612 ERROR_exit("Illegal mode after -final\n");
00613 nopt++ ; continue ;
00614 }
00615
00616
00617
00618 if( strcmp(argv[nopt],"-parfix") == 0 ){
00619 int np , ip ; float vp ;
00620 if( ++nopt >= argc-1 )
00621 ERROR_exit("Need 2 parameters afer -parfix!\n");
00622 np = strtol( argv[nopt] , NULL , 10 ) ; nopt++ ;
00623 vp = strtod( argv[nopt] , NULL ) ;
00624 if( np <= 0 || np > MAXPAR )
00625 ERROR_exit("Param #%d after -parfix is illegal!\n",np) ;
00626 for( ip=0 ; ip < nparfix ; ip++ ){
00627 if( parfix[ip].np == np ){
00628 WARNING_message("Ignoring later -parfix option for param #%d\n",ip);
00629 break ;
00630 }
00631 }
00632 if( ip == nparfix ) nparfix++ ;
00633 parfix[ip].np = np ; parfix[ip].vp = vp ;
00634 nopt++ ; continue ;
00635 }
00636
00637
00638
00639 if( strcmp(argv[nopt],"-shift_only") == 0 ){
00640 warpdrive_code = WARPDRIVE_SHIFT ; nopt++ ; continue ;
00641 }
00642
00643 if( strcmp(argv[nopt],"-shift_rotate") == 0 ){
00644 warpdrive_code = WARPDRIVE_ROTATE ; nopt++ ; continue ;
00645 }
00646
00647 if( strcmp(argv[nopt],"-shift_rotate_scale") == 0 ){
00648 warpdrive_code = WARPDRIVE_SCALE ; nopt++ ; continue ;
00649 }
00650
00651 if( strcmp(argv[nopt],"-affine_general") == 0 ){
00652 warpdrive_code = WARPDRIVE_AFFINE ; nopt++ ; continue ;
00653 }
00654
00655 if( strcmp(argv[nopt],"-bilinear_general") == 0 ){
00656 #if 0
00657 ERROR_exit("3dWarpDrive -bilinear_general NOT IMPLEMENTED!\n");
00658 #else
00659 WARNING_message("-bilinear_general transformations are experimental!") ;
00660 #endif
00661 warpdrive_code = WARPDRIVE_BILINEAR ; nopt++ ; continue ;
00662 }
00663
00664
00665
00666 if( strcmp(argv[nopt],"-NN") == 0 ){
00667 abas.regmode = MRI_NN ; nopt++ ; continue ;
00668 }
00669 if( strcmp(argv[nopt],"-linear") == 0 ){
00670 abas.regmode = MRI_LINEAR ; nopt++ ; continue ;
00671 }
00672 if( strcmp(argv[nopt],"-cubic") == 0 ){
00673 abas.regmode = MRI_CUBIC ; nopt++ ; continue ;
00674 }
00675 if( strcmp(argv[nopt],"-quintic") == 0 ){
00676 abas.regmode = MRI_QUINTIC ; nopt++ ; continue ;
00677 }
00678
00679
00680
00681 if( strcmp(argv[nopt],"-prefix") == 0 ){
00682 if( ++nopt >= argc )
00683 ERROR_exit("Need an argument after -prefix!\n");
00684 if( !THD_filename_ok(argv[nopt]) )
00685 ERROR_exit("-prefix argument is invalid!\n");
00686 prefix = argv[nopt] ; nopt++ ; continue ;
00687 }
00688
00689
00690
00691 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00692 abas.verb++ ; nopt++ ; continue ;
00693 }
00694
00695
00696
00697 if( strcmp(argv[nopt],"-base") == 0 ){
00698 if( ++nopt >= argc )
00699 ERROR_exit("Need an argument after -base!\n");
00700 baset = THD_open_dataset( argv[nopt] ) ;
00701 if( baset == NULL )
00702 ERROR_exit("Can't open -base dataset %s\n",argv[nopt]);
00703 if( DSET_NVALS(baset) > 1 )
00704 WARNING_message(
00705 "-base dataset %s has %d sub-bricks; will only use #0\n",
00706 argv[nopt],DSET_NVALS(baset) ) ;
00707 nopt++ ; continue ;
00708 }
00709
00710
00711
00712 if( strcmp(argv[nopt],"-weight") == 0 ){
00713 if( ++nopt >= argc )
00714 ERROR_exit("Need an argument after -weight!\n");
00715 wtset = THD_open_dataset( argv[nopt] ) ;
00716 if( wtset == NULL )
00717 ERROR_exit("Can't open -weight dataset %s\n",argv[nopt]);
00718 if( DSET_NVALS(wtset) > 1 )
00719 WARNING_message(
00720 "-weight dataset %s has %d sub-bricks; will only use #0\n",
00721 argv[nopt],DSET_NVALS(wtset) ) ;
00722 nopt++ ; continue ;
00723 }
00724
00725
00726
00727 if( strcmp(argv[nopt],"-input") == 0 ){
00728 if( ++nopt >= argc )
00729 ERROR_exit("Need an argument after -input!\n");
00730 inset = THD_open_dataset( argv[nopt] ) ;
00731 if( inset == NULL )
00732 ERROR_exit("Can't open -input dataset %s\n",argv[nopt]);
00733 nopt++ ; continue ;
00734 }
00735
00736
00737
00738 if( strcmp(argv[nopt],"-maxite") == 0 ){
00739 int ival ;
00740 if( ++nopt >= argc )
00741 ERROR_exit("Need an argument after -maxite!\n");
00742 ival = strtol( argv[nopt] , NULL , 10 ) ;
00743 if( ival > 1 ) abas.max_iter = ival ;
00744 nopt++ ; continue ;
00745 }
00746
00747
00748
00749 if( strcmp(argv[nopt],"-delta") == 0 ){
00750 float val ;
00751 if( ++nopt >= argc )
00752 ERROR_exit("Need an argument after -delta!\n");
00753 val = strtod( argv[nopt] , NULL ) ;
00754 if( val > 0.0499 && val < 49.99 ) abas.delfac = val ;
00755 nopt++ ; continue ;
00756 }
00757
00758
00759
00760 if( strcmp(argv[nopt],"-thresh") == 0 ){
00761 float val ;
00762 if( ++nopt >= argc )
00763 ERROR_exit("Need an argument after -thresh!\n");
00764 val = strtod( argv[nopt] , NULL ) ;
00765 if( val > 0.001 && val < 3.01 ) abas.tolfac = val ;
00766 nopt++ ; continue ;
00767 }
00768
00769
00770
00771 ERROR_exit("Unknown option %s\n",argv[nopt]) ;
00772
00773 }
00774
00775
00776
00777 if( inset == NULL && nopt != argc-1 )
00778 ERROR_exit("Command line should have exactly 1 dataset!\n"
00779 "** Whereas there seems to be %d of them!\n",
00780 argc-nopt ) ;
00781
00782
00783
00784 if( inset == NULL ){
00785 inset = THD_open_dataset( argv[nopt] ) ;
00786 if( !ISVALID_DSET(inset) )
00787 ERROR_exit("Can't open dataset %s\n",argv[nopt]);
00788 }
00789
00790 nx = DSET_NX(inset) ; ny = DSET_NY(inset) ; nz = DSET_NZ(inset) ;
00791 dx = DSET_DX(inset) ; dy = DSET_DY(inset) ; dz = DSET_DZ(inset) ;
00792
00793 if( abas.verb ) INFO_message("Checking inputs\n") ;
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 #define ADDPAR(nm,bb,tt,id,dd,ll) \
00806 do{ int p=abas.nparam ; \
00807 abas.param = (MRI_warp3D_param_def *) realloc( \
00808 (void *)abas.param , \
00809 sizeof(MRI_warp3D_param_def)*(p+1) ) ; \
00810 abas.param[p].min = (bb) ; abas.param[p].max = (tt) ; \
00811 abas.param[p].delta = (dd) ; abas.param[p].toler = (ll) ; \
00812 abas.param[p].ident = abas.param[p].val_init = (id) ; \
00813 strcpy( abas.param[p].name , (nm) ) ; \
00814 abas.param[p].fixed = 0 ; \
00815 abas.nparam = p+1 ; \
00816 } while(0)
00817
00818 nerr = 0 ;
00819 if( warpdrive_code <= 0 ){
00820 ERROR_message("Need a transform-specifying option!\n");
00821 nerr++ ;
00822 } else if( WARPDRIVE_IS_AFFINE(warpdrive_code) ) {
00823
00824 char *lab09, *lab10, *lab11 ;
00825
00826
00827
00828 ADDPAR( "x-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00829 ADDPAR( "y-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00830 ADDPAR( "z-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00831
00832 ADDPAR( "z-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00833 ADDPAR( "x-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00834 ADDPAR( "y-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00835
00836 ADDPAR( "x-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00837 ADDPAR( "y-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00838 ADDPAR( "z-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00839
00840 switch( smat ){
00841 default:
00842 case SMAT_LOWER:
00843 lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
00844 break ;
00845
00846 case SMAT_UPPER:
00847 lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
00848 break ;
00849 }
00850 ADDPAR( lab09 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00851 ADDPAR( lab10 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00852 ADDPAR( lab11 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00853
00854
00855
00856 for( kpar=0 ; kpar < 12 ; kpar++ )
00857 parvec[kpar] = abas.param[kpar].ident ;
00858
00859
00860
00861 if( warpdrive_code == WARPDRIVE_SHIFT ){
00862 warp_parset = parset_shift ;
00863 warp_for = warper_shift_for ;
00864 warp_inv = warper_shift_inv ;
00865 } else {
00866 warp_parset = parset_affine ;
00867 warp_for = warper_affine_for ;
00868 warp_inv = warper_affine_inv ;
00869 }
00870
00871
00872
00873 switch( warpdrive_code ){
00874 case WARPDRIVE_SHIFT: abas.nparam = 3 ; break ;
00875 case WARPDRIVE_ROTATE: abas.nparam = 6 ; break ;
00876 case WARPDRIVE_SCALE: abas.nparam = 9 ; break ;
00877 case WARPDRIVE_AFFINE: abas.nparam = 12 ; break ;
00878 }
00879
00880 } else if( warpdrive_code == WARPDRIVE_BILINEAR ){
00881
00882 char *lab09, *lab10, *lab11 , labxx[16] ;
00883 float xr,yr,zr,rr ;
00884
00885
00886
00887 ADDPAR( "x-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00888 ADDPAR( "y-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00889 ADDPAR( "z-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00890
00891 ADDPAR( "z-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00892 ADDPAR( "x-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00893 ADDPAR( "y-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00894
00895 ADDPAR( "x-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00896 ADDPAR( "y-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00897 ADDPAR( "z-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00898
00899 switch( smat ){
00900 default:
00901 case SMAT_LOWER:
00902 lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
00903 break ;
00904
00905 case SMAT_UPPER:
00906 lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
00907 break ;
00908 }
00909 ADDPAR( lab09 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00910 ADDPAR( lab10 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00911 ADDPAR( lab11 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00912
00913 xr = 0.5f*fabs(dx)*nx ; yr = 0.5f*fabs(dy)*ny ; zr = 0.5f*fabs(dz)*nz ;
00914 rr = MAX(xr,yr) ; rr = MAX(rr,zr) ; dd_fac = 1.0f / rr ;
00915 for( kpar=12 ; kpar < 39 ; kpar++ ){
00916 sprintf(labxx,"blin_%02d",kpar+1) ;
00917 ADDPAR( labxx , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00918 }
00919
00920
00921
00922 for( kpar=0 ; kpar < 39 ; kpar++ )
00923 parvec[kpar] = abas.param[kpar].ident ;
00924
00925 abas.vwdet = warper_bilinear_det ;
00926
00927
00928
00929 warp_parset = parset_bilinear ;
00930 warp_for = warper_bilinear_for ;
00931 warp_inv = warper_bilinear_inv ;
00932
00933
00934
00935 abas.nparam = 39 ;
00936
00937 } else {
00938 ERROR_message("Unimplemented transform model!\n") ;
00939 nerr++ ;
00940 }
00941
00942
00943
00944 nfree = abas.nparam ;
00945 for( kpar=0 ; kpar < nparfix ; kpar++ ){
00946 np = parfix[kpar].np - 1 ; vp = parfix[kpar].vp ;
00947 if( np >= 0 && np < abas.nparam ){
00948 if( vp >= abas.param[np].min && vp <= abas.param[np].max ){
00949 abas.param[np].fixed = 1 ;
00950 abas.param[np].val_fixed = vp ;
00951 nfree -- ;
00952 } else {
00953 ERROR_message("-parfix for param #%d has illegal value!\n",np+1) ;
00954 nerr++ ;
00955 }
00956 } else {
00957 WARNING_message("-parfix for param #%d is out of range 1..%d\n",
00958 np+1 , abas.nparam+1 ) ;
00959 }
00960 }
00961 if( nfree <= 0 ){
00962 ERROR_message("No free parameters in transform model!\n") ;
00963 nerr++ ;
00964 }
00965
00966
00967
00968 if( abas.max_iter <= 0 ) abas.max_iter = 11*nfree+5 ;
00969
00970
00971
00972 if( baset == NULL ){
00973 ERROR_message("Need to specify a base dataset!\n") ;
00974 nerr++ ;
00975 }
00976
00977 if( nerr ) exit(1) ;
00978
00979 if( abas.verb ) INFO_message("Creating empty output dataset\n") ;
00980
00981 outset = EDIT_empty_copy( inset ) ;
00982
00983 EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;
00984
00985 if( output_float ){
00986 EDIT_dset_items( outset , ADN_datum_all , MRI_float , ADN_none ) ;
00987 for( kim=0 ; kim < nvals ; kim++ )
00988 EDIT_BRICK_FACTOR( outset , kim , 0.0 ) ;
00989 }
00990
00991 if( THD_is_file( DSET_HEADNAME(outset) ) )
00992 ERROR_exit("Output file %s already exists!\n",DSET_HEADNAME(outset) ) ;
00993
00994 tross_Copy_History( inset , outset ) ;
00995 tross_Make_History( "3dWarpDrive" , argc,argv , outset ) ;
00996
00997 outset->daxes->xxorg = inset->daxes->xxorg ;
00998 outset->daxes->yyorg = inset->daxes->yyorg ;
00999 outset->daxes->zzorg = inset->daxes->zzorg ;
01000
01001
01002
01003 if( DSET_NX(baset) != nx || DSET_NY(baset) != ny || DSET_NZ(baset) != nz ){
01004 ERROR_message("base and input grid dimensions don't match!\n") ;
01005 ERROR_message("base is %d X %d X %d voxels\n",
01006 DSET_NX(baset),DSET_NY(baset),DSET_NZ(baset) ) ;
01007 ERROR_exit ("input is %d X %d X %d voxels\n",nx,ny,nz) ;
01008 }
01009
01010
01011
01012 if( FLDIF(DSET_DX(baset),dx) ||
01013 FLDIF(DSET_DY(baset),dy) || FLDIF(DSET_DZ(baset),dz) ){
01014 WARNING_message("base and input grid spacings don't match!\n") ;
01015 WARNING_message("base grid = %.5f X %.5f X %.5f mm\n",
01016 DSET_DX(baset),DSET_DY(baset),DSET_DZ(baset) ) ;
01017 WARNING_message("input grid = %.5f X %.5f X %.5f mm\n",
01018 DSET_DX(inset),DSET_DY(inset),DSET_DZ(inset) ) ;
01019 }
01020
01021 if( FLDIF(DSET_XORG(baset),DSET_XORG(inset)) ||
01022 FLDIF(DSET_YORG(baset),DSET_YORG(inset)) ||
01023 FLDIF(DSET_ZORG(baset),DSET_ZORG(inset)) ){
01024 WARNING_message("base and input grid offsets don't match!\n") ;
01025 WARNING_message("base offsets = %.5f X %.5f X %.5f mm\n",
01026 DSET_XORG(baset),DSET_YORG(baset),DSET_ZORG(baset) ) ;
01027 WARNING_message("input offsets = %.5f X %.5f X %.5f mm\n",
01028 DSET_XORG(inset),DSET_YORG(inset),DSET_ZORG(inset) ) ;
01029 }
01030
01031 if( baset->daxes->xxorient != inset->daxes->xxorient ||
01032 baset->daxes->yyorient != inset->daxes->yyorient ||
01033 baset->daxes->zzorient != inset->daxes->zzorient ){
01034 WARNING_message("base and input orientations don't match!\n") ;
01035 WARNING_message("base = %s X %s X %s\n",
01036 ORIENT_shortstr[baset->daxes->xxorient] ,
01037 ORIENT_shortstr[baset->daxes->yyorient] ,
01038 ORIENT_shortstr[baset->daxes->zzorient] ) ;
01039 WARNING_message("input = %s X %s X %s\n",
01040 ORIENT_shortstr[inset->daxes->xxorient] ,
01041 ORIENT_shortstr[inset->daxes->yyorient] ,
01042 ORIENT_shortstr[inset->daxes->zzorient] ) ;
01043 }
01044
01045
01046
01047 if( wtset != NULL &&
01048 (DSET_NX(wtset) != nx || DSET_NY(wtset) != ny || DSET_NZ(wtset) != nz) ){
01049 ERROR_message("weight and input grid dimensions don't match!\n") ;
01050 ERROR_message("weight is %d X %d X %d voxels\n",
01051 DSET_NX(wtset),DSET_NY(wtset),DSET_NZ(wtset) ) ;
01052 ERROR_exit ("input is %d X %d X %d voxels\n",nx,ny,nz) ;
01053 }
01054
01055
01056
01057 if( abas.verb ) INFO_message("Loading datasets\n") ;
01058
01059 DSET_load(inset) ;
01060 if( !DSET_LOADED(inset) ){
01061 ERROR_exit("Can't load input dataset into memory!\n") ;
01062 } else {
01063 nvals = DSET_NVALS(inset) ;
01064 if( nvals == 1 ){
01065 clip_inset = THD_cliplevel( DSET_BRICK(inset,0) , 0.0 ) ;
01066 if( DSET_BRICK_FACTOR(inset,0) > 0.0f )
01067 clip_inset *= DSET_BRICK_FACTOR(inset,0) ;
01068 } else {
01069 qim = THD_median_brick( inset ) ;
01070 clip_inset = THD_cliplevel( qim , 0.0 ) ;
01071 mri_free(qim) ;
01072 }
01073 }
01074
01075 DSET_load(baset) ;
01076 if( !DSET_LOADED(baset) ){
01077 ERROR_exit("Can't load base dataset into memory!\n") ;
01078 } else {
01079 clip_baset = THD_cliplevel( DSET_BRICK(baset,0) , 0.0 ) ;
01080 abas.imbase = mri_to_float( DSET_BRICK(baset,0) ) ;
01081 base_idc = strdup(baset->idcode.str) ;
01082 DSET_delete(baset) ; baset = NULL ;
01083 }
01084
01085 if( wtset != NULL ){
01086 DSET_load(wtset) ;
01087 if( !DSET_LOADED(wtset) ){
01088 ERROR_exit("Can't load weight dataset into memory!\n") ;
01089 } else {
01090 abas.imwt = mri_to_float( DSET_BRICK(wtset,0) ) ;
01091 wt_idc = strdup(wtset->idcode.str) ;
01092 DSET_delete(wtset) ; wtset = NULL ;
01093 }
01094 }
01095
01096
01097
01098 { THD_vecmat ijk_to_inset_xyz , xyz_to_dicom ;
01099
01100 LOAD_DIAG_MAT( ijk_to_inset_xyz.mm ,
01101 inset->daxes->xxdel ,
01102 inset->daxes->yydel , inset->daxes->zzdel );
01103
01104 if( warpdrive_code == WARPDRIVE_BILINEAR ){
01105
01106
01107 LOAD_FVEC3 ( ijk_to_inset_xyz.vv ,
01108 -0.5*(nx-1)*inset->daxes->xxdel ,
01109 -0.5*(ny-1)*inset->daxes->yydel ,
01110 -0.5*(nz-1)*inset->daxes->zzdel ) ;
01111
01112 } else {
01113
01114
01115 LOAD_FVEC3 ( ijk_to_inset_xyz.vv ,
01116 DSET_XORG(inset) , DSET_YORG(inset), DSET_ZORG(inset) ) ;
01117 }
01118
01119 xyz_to_dicom.mm = inset->daxes->to_dicomm ;
01120 LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
01121
01122 ijk_to_xyz = MUL_VECMAT( xyz_to_dicom , ijk_to_inset_xyz ) ;
01123 xyz_to_ijk = INV_VECMAT( ijk_to_xyz ) ;
01124
01125 #if 0
01126 if( abas.verb ){
01127 DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
01128 DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
01129 }
01130 #endif
01131 }
01132
01133
01134
01135 if( clip_baset > 0.0f && clip_inset > 0.0f ){
01136 float fac = clip_inset / clip_baset ;
01137 if( fac > 0.005 && fac < 2000.0 ){
01138 abas.scale_init = fac ;
01139 if( abas.verb ) INFO_message("Scale factor set to %.2f/%.2f=%.2g\n",
01140 clip_baset , clip_inset , fac ) ;
01141 if (fac > 200) {
01142 fprintf(stderr,"++ Warning: Scale init = %g is large. Check output.\n", fac) ;
01143 }
01144 } else {
01145 fprintf(stderr,"-- Large scale difference between datasets.\n"
01146 " Scale init = %g\n"
01147 " 3dWarpDrive might not converge.",fac) ;
01148 }
01149 }
01150
01151
01152
01153 if( abas.verb ) INFO_message("Beginning alignment setup\n") ;
01154
01155
01156
01157 parsave = (float **)malloc( sizeof(float *) * abas.nparam ) ;
01158 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01159 parsave[kpar] = (float *)calloc( sizeof(float) , nvals ) ;
01160
01161 mri_warp3D_align_setup( &abas ) ;
01162
01163 if( abas.verb ) INFO_message("Beginning alignment loop\n") ;
01164
01165
01166
01167 for( kim=0 ; kim < nvals ; kim++ ){
01168
01169 for( kpar=0 ; kpar < abas.nparam ; kpar++ ){
01170 if( abas.param[kpar].fixed )
01171 abas.param[kpar].val_init = abas.param[kpar].val_fixed ;
01172 else
01173 abas.param[kpar].val_init = abas.param[kpar].ident ;
01174 }
01175
01176
01177
01178
01179 qim = mri_scale_to_float( DSET_BRICK_FACTOR(inset,kim) ,
01180 DSET_BRICK(inset,kim) ) ;
01181 tim = mri_warp3d_align_one( &abas , qim ) ;
01182 mri_free( qim ) ; DSET_unload_one( inset , kim ) ;
01183
01184 if( abas.verb ){ DUMP_VECMAT( "end mv_for" , mv_for ) ; }
01185
01186
01187
01188 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01189 parsave[kpar][kim] = abas.param[kpar].val_out ;
01190
01191
01192
01193 switch( DSET_BRICK_TYPE(outset,kim) ){
01194
01195 default:
01196 ERROR_message("Can't store bricks of type %s\n",
01197 MRI_TYPE_name[DSET_BRICK_TYPE(inset,kim)] ) ;
01198
01199
01200 case MRI_float:
01201 EDIT_substitute_brick( outset, kim, MRI_float, MRI_FLOAT_PTR(tim) );
01202 mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
01203 break ;
01204
01205 case MRI_short:{
01206 double fac = DSET_BRICK_FACTOR(inset,kim) ;
01207 fac = (fac > 0.0) ? 1.0/fac : 1.0 ;
01208 fim = mri_to_short(fac,tim) ; mri_free( tim ) ;
01209 EDIT_substitute_brick( outset, kim, MRI_short, MRI_SHORT_PTR(fim) );
01210 fac = (fac != 1.0) ? 1.0/fac : 0.0 ;
01211 EDIT_BRICK_FACTOR( outset , kim , fac ) ;
01212 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
01213 }
01214 break ;
01215
01216 case MRI_byte:
01217 vp = mri_min(tim) ;
01218 if( vp < 0.0f ){
01219 WARNING_message(
01220 "output sub-brick #%d is byte, but has negative values\n",
01221 kim ) ;
01222 }
01223 fim = mri_to_byte(tim) ; mri_free( tim ) ;
01224 EDIT_substitute_brick( outset, kim, MRI_byte, MRI_BYTE_PTR(fim) ) ;
01225 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
01226 break ;
01227 }
01228 }
01229 DSET_unload( inset ) ;
01230
01231
01232
01233 mri_warp3D_align_cleanup( &abas ) ;
01234
01235
01236
01237
01238 THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_IDCODE" ,
01239 inset->idcode.str ) ;
01240 THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_NAME" ,
01241 DSET_HEADNAME(inset) ) ;
01242 if( base_idc != NULL )
01243 THD_set_string_atr( outset->dblk , "WARPDRIVE_BASE_IDCODE" , base_idc ) ;
01244 if( wt_idc != NULL )
01245 THD_set_string_atr( outset->dblk , "WARPDRIVE_WEIGHT_IDCODE" , wt_idc ) ;
01246
01247 if( WARPDRIVE_IS_AFFINE(warpdrive_code) ){
01248 float matar[12] ; char anam[64] ;
01249
01250 for( kpar=0 ; kpar < 12 ; kpar++ ) parvec[kpar] = 0.0 ;
01251
01252 for( kim=0 ; kim < nvals ; kim++ ){
01253 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01254 parvec[kpar] = parsave[kpar][kim] ;
01255 parset_affine() ;
01256
01257 UNLOAD_MAT(mv_for.mm,matar[0],matar[1],matar[2],
01258 matar[4],matar[5],matar[6],
01259 matar[8],matar[9],matar[10] ) ;
01260 UNLOAD_FVEC3(mv_for.vv,matar[3],matar[7],matar[11]) ;
01261 sprintf(anam,"WARPDRIVE_MATVEC_FOR_%06d",kim) ;
01262 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
01263
01264 UNLOAD_MAT(mv_inv.mm,matar[0],matar[1],matar[2],
01265 matar[4],matar[5],matar[6],
01266 matar[8],matar[9],matar[10] ) ;
01267 UNLOAD_FVEC3(mv_inv.vv,matar[3],matar[7],matar[11]) ;
01268 sprintf(anam,"WARPDRIVE_MATVEC_INV_%06d",kim) ;
01269 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
01270 }
01271 }
01272
01273
01274
01275 DSET_write( outset ) ; DSET_unload( outset ) ;
01276 if( abas.verb ) WROTE_DSET(outset) ;
01277
01278 if( W_1Dfile != NULL ){
01279 FILE *fp ;
01280 if( abas.verb ) INFO_message("Writing 1Dfile: %s\n",W_1Dfile) ;
01281 if( THD_is_file(W_1Dfile) )
01282 WARNING_message("Overwriting file %s\n",W_1Dfile) ;
01283
01284 fp = fopen( W_1Dfile , "w" ) ;
01285 if( fp != NULL ){
01286
01287 fprintf(fp,"#") ;
01288 for( kim=0 ; kim < argc ; kim++ ) fprintf(fp," %s",argv[kim]) ;
01289 fprintf(fp,"\n") ;
01290
01291 fprintf(fp,"#") ;
01292 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01293 fprintf(fp," %-13.13s",abas.param[kpar].name) ;
01294 fprintf(fp,"\n") ;
01295
01296 fprintf(fp,"#") ;
01297 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01298 fprintf(fp," -------------") ;
01299 fprintf(fp,"\n") ;
01300
01301 for( kim=0 ; kim < nvals ; kim++ ){
01302 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01303 fprintf(fp," %13.6g",parsave[kpar][kim]) ;
01304 fprintf(fp,"\n") ;
01305 }
01306 fclose(fp) ;
01307 }
01308
01309 }
01310
01311 if( abas.verb ){
01312 double tt = (NI_clock_time()-ctstart)*0.001 ;
01313 INFO_message("Total elapsed time = %.2f s\n",tt) ;
01314 }
01315
01316 exit(0) ;
01317 }
