1 /*------------------------------------------------------------------------
2 *** This program does something, but nobody is quite sure what. ***
3 *** But what it does is very important; nobody doubts that, either. ***
4 --------------------------------------------------------------------------*/
5
6 #include "mrilib.h"
7
8 /*--------------------------------------------------------------------------*/
9 /* 02 Dec 2005: prototype for quick initializer for shift+rotate parameters */
10
11 static float get_best_shiftrot( MRI_warp3D_align_basis *bas ,
12 MRI_IMAGE *base , MRI_IMAGE *vol ,
13 float *roll , float *pitch , float *yaw ,
14 int *dxp , int *dyp , int *dzp ) ;
15
16 static int VL_coarse_del = 12 ; /* variables for the above, */
17 static int VL_coarse_num = 3 ; /* copied from 3dvolreg.c */
18 static int VL_coarse_rot = 0 ;
19 /*--------------------------------------------------------------------------*/
20
21 #define MAXPAR 99
22
23 typedef struct { int np ; float vp ; } fixed_param ;
24
25 static int nparfix ;
26 static fixed_param parfix[MAXPAR] ;
27
28 static float parvec[MAXPAR] ;
29
30 static void (*warp_parset)(void) = NULL ;
31
load_parvec(int np,float * pv)32 void load_parvec( int np, float *pv ){
33 memcpy( parvec , pv , sizeof(float)*np ) ;
34 if( warp_parset != NULL ) warp_parset() ;
35 }
36
37 /*--------------------------------------------------------------------------*/
38
39 static void (*warp_for)( float,float,float , float *,float *,float *) ;
40 static void (*warp_inv)( float,float,float , float *,float *,float *) ;
41
42 static THD_vecmat ijk_to_xyz , xyz_to_ijk ;
43
44 /*--------------------------------------------------------------------------*/
45
ijk_warp_for(float ii,float jj,float kk,float * pp,float * qq,float * rr)46 void ijk_warp_for( float ii , float jj , float kk ,
47 float *pp , float *qq , float *rr )
48 {
49 THD_fvec3 xxx , yyy ;
50
51 LOAD_FVEC3( xxx , ii,jj,kk ) ;
52 yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
53 warp_for( yyy.xyz[0] , yyy.xyz[1] , yyy.xyz[2] ,
54 &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
55 yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
56 *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
57 }
58
59 /*--------------------------------------------------------------------------*/
60
ijk_warp_inv(float ii,float jj,float kk,float * pp,float * qq,float * rr)61 void ijk_warp_inv( float ii , float jj , float kk ,
62 float *pp , float *qq , float *rr )
63 {
64 THD_fvec3 xxx , yyy ;
65
66 LOAD_FVEC3( xxx , ii,jj,kk ) ;
67 yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
68 warp_inv( yyy.xyz[0] , yyy.xyz[1] , yyy.xyz[2] ,
69 &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
70 yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
71 *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
72 }
73
74 /*--------------------------------------------------------------------------*/
75
76 #define WARPDRIVE_SHIFT 1
77 #define WARPDRIVE_ROTATE 2
78 #define WARPDRIVE_SCALE 3
79 #define WARPDRIVE_AFFINE 4
80 #define WARPDRIVE_BILINEAR 5
81
82 #define WARPDRIVE_IS_AFFINE(wc) \
83 ( (wc) >= WARPDRIVE_SHIFT && (wc) <= WARPDRIVE_AFFINE )
84
85 /*--------------------------------------------------------------------------*/
86 /* For shift-only 'warps' */
87
88 static float xsh , ysh , zsh ;
89
parset_shift(void)90 void parset_shift(void)
91 {
92 xsh = parvec[0] ; ysh = parvec[1] ; zsh = parvec[2] ;
93 }
94
warper_shift_for(float aa,float bb,float cc,float * p,float * q,float * r)95 void warper_shift_for( float aa , float bb , float cc ,
96 float *p , float *q , float *r )
97 {
98 *p = aa+xsh ; *q = bb+ysh ; *r = cc+zsh ;
99 }
100
warper_shift_inv(float aa,float bb,float cc,float * p,float * q,float * r)101 void warper_shift_inv( float aa , float bb , float cc ,
102 float *p , float *q , float *r )
103 {
104 *p = aa-xsh ; *q = bb-ysh ; *r = cc-zsh ;
105 }
106
107 /*--------------------------------------------------------------------------*/
108 /* For affine warps */
109
110 static THD_vecmat mv_for , mv_inv ;
111
warper_affine_for(float aa,float bb,float cc,float * p,float * q,float * r)112 void warper_affine_for( float aa , float bb , float cc ,
113 float *p , float *q , float *r )
114 {
115 THD_fvec3 v , w ;
116 LOAD_FVEC3( v , aa,bb,cc ) ;
117 w = VECMAT_VEC( mv_for , v ) ;
118 *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
119 }
120
warper_affine_inv(float aa,float bb,float cc,float * p,float * q,float * r)121 void warper_affine_inv( float aa , float bb , float cc ,
122 float *p , float *q , float *r )
123 {
124 THD_fvec3 v , w ;
125 LOAD_FVEC3( v , aa,bb,cc ) ;
126 w = VECMAT_VEC( mv_inv , v ) ;
127 *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
128 }
129
130 /*--------------------------------------------------------------------------
131 Compute a rotation matrix specified by 3 angles:
132 Q = R3 R2 R1, where Ri is rotation about axis axi by angle thi.
133 ----------------------------------------------------------------------------*/
134
rot_matrix(int ax1,double th1,int ax2,double th2,int ax3,double th3)135 static THD_mat33 rot_matrix( int ax1, double th1,
136 int ax2, double th2, int ax3, double th3 )
137 {
138 THD_mat33 q , p ;
139 LOAD_ROT_MAT( q , th1 , ax1 ) ;
140 LOAD_ROT_MAT( p , th2 , ax2 ) ; q = MAT_MUL( p , q ) ;
141 LOAD_ROT_MAT( p , th3 , ax3 ) ; q = MAT_MUL( p , q ) ;
142 return q ;
143 }
144
145 #define D2R (PI/180.0) /* angles are in degrees */
146
147 static int matorder = MATORDER_SDU ; /* cf. mrilib.h */
148 static int smat = SMAT_LOWER ;
149 static int dcode = DELTA_AFTER ; /* cf. 3ddata.h */
150
parset_affine(void)151 void parset_affine(void)
152 {
153 THD_mat33 ss,dd,uu,aa,bb ;
154 THD_fvec3 vv ;
155 float a,b,c ;
156
157 #if 0
158 { int ii;fprintf(stderr,"\nparset:");
159 for(ii=0;ii<12;ii++)fprintf(stderr," %g",parvec[ii]); fprintf(stderr,"\n");}
160 #endif
161 memset(&bb,0,sizeof(THD_mat33)); memset(&aa,0,sizeof(THD_mat33));
162 /* rotation */
163
164 uu = rot_matrix( 2, D2R*parvec[3] , 0, D2R*parvec[4] , 1, D2R*parvec[5] ) ;
165
166 /* scaling */
167
168 a = parvec[6] ; if( a <= 0.01f ) a = 1.0f ;
169 b = parvec[7] ; if( b <= 0.01f ) b = 1.0f ;
170 c = parvec[8] ; if( c <= 0.01f ) c = 1.0f ;
171 LOAD_DIAG_MAT( dd , a,b,c ) ;
172
173 /* shear */
174
175 a = parvec[ 9] ; if( fabs(a) > 0.3333f ) a = 0.0f ;
176 b = parvec[10] ; if( fabs(b) > 0.3333f ) b = 0.0f ;
177 c = parvec[11] ; if( fabs(c) > 0.3333f ) c = 0.0f ;
178
179 switch( smat ){
180 default:
181 case SMAT_LOWER:
182 LOAD_MAT( ss , 1.0 , 0.0 , 0.0 ,
183 a , 1.0 , 0.0 ,
184 b , c , 1.0 ) ;
185 break ;
186
187 case SMAT_UPPER:
188 LOAD_MAT( ss , 1.0 , a , b ,
189 0.0 , 1.0 , c ,
190 0.0 , 0.0 , 1.0 ) ;
191 break ;
192 }
193
194 /* multiply them, as ordered */
195
196 switch( matorder ){
197 case MATORDER_SDU: aa = MAT_MUL(ss,dd) ; bb = uu ; break ;
198 case MATORDER_SUD: aa = MAT_MUL(ss,uu) ; bb = dd ; break ;
199 case MATORDER_DSU: aa = MAT_MUL(dd,ss) ; bb = uu ; break ;
200 case MATORDER_DUS: aa = MAT_MUL(dd,uu) ; bb = ss ; break ;
201 case MATORDER_USD: aa = MAT_MUL(uu,ss) ; bb = dd ; break ;
202 case MATORDER_UDS: aa = MAT_MUL(uu,dd) ; bb = ss ; break ;
203 }
204 mv_for.mm = MAT_MUL(aa,bb) ;
205
206 LOAD_FVEC3( vv , parvec[0] , parvec[1] , parvec[2] ) ;
207
208 switch( dcode ){
209 case DELTA_AFTER: mv_for.vv = vv ; break ;
210 case DELTA_BEFORE: mv_for.vv = MATVEC( mv_for.mm , vv ) ; break ;
211 }
212
213 mv_inv = INV_VECMAT( mv_for ) ;
214
215 #if 0
216 DUMP_VECMAT("mv_for",mv_for) ; DUMP_VECMAT("mv_inv",mv_inv) ;
217 #endif
218 }
219
220 /*--------------------------------------------------------------------------*/
221 /* For bilinear warps */
222
223 static float dd_fac = 1.0f ;
224 static float dd_for[3][3][3] , dd_inv[3][3][3] ;
225
warper_bilinear_for(float aa,float bb,float cc,float * p,float * q,float * r)226 void warper_bilinear_for( float aa , float bb , float cc ,
227 float *p , float *q , float *r )
228 {
229 THD_fvec3 v,w ; THD_mat33 dd,ee ;
230
231 LOAD_FVEC3( v , aa,bb,cc ) ;
232 w = VECMAT_VEC( mv_for , v ) ;
233
234 dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
235 dd.mat[0][1] = dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
236 dd.mat[0][2] = dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
237 dd.mat[1][0] = dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
238 dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
239 dd.mat[1][2] = dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
240 dd.mat[2][0] = dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
241 dd.mat[2][1] = dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
242 dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
243
244 ee = MAT_INV(dd) ;
245 v = MATVEC(ee,w) ;
246
247 *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
248 }
249
warper_bilinear_inv(float aa,float bb,float cc,float * p,float * q,float * r)250 void warper_bilinear_inv( float aa , float bb , float cc ,
251 float *p , float *q , float *r )
252 {
253 THD_fvec3 v,w ; THD_mat33 dd,ee ;
254
255 LOAD_FVEC3( v , aa,bb,cc ) ;
256 w = VECMAT_VEC( mv_inv , v ) ;
257
258 dd.mat[0][0] = 1.0f + dd_inv[0][0][0]*aa + dd_inv[0][0][1]*bb + dd_inv[0][0][2]*cc ;
259 dd.mat[0][1] = dd_inv[0][1][0]*aa + dd_inv[0][1][1]*bb + dd_inv[0][1][2]*cc ;
260 dd.mat[0][2] = dd_inv[0][2][0]*aa + dd_inv[0][2][1]*bb + dd_inv[0][2][2]*cc ;
261 dd.mat[1][0] = dd_inv[1][0][0]*aa + dd_inv[1][0][1]*bb + dd_inv[1][0][2]*cc ;
262 dd.mat[1][1] = 1.0f + dd_inv[1][1][0]*aa + dd_inv[1][1][1]*bb + dd_inv[1][1][2]*cc ;
263 dd.mat[1][2] = dd_inv[1][2][0]*aa + dd_inv[1][2][1]*bb + dd_inv[1][2][2]*cc ;
264 dd.mat[2][0] = dd_inv[2][0][0]*aa + dd_inv[2][0][1]*bb + dd_inv[2][0][2]*cc ;
265 dd.mat[2][1] = dd_inv[2][1][0]*aa + dd_inv[2][1][1]*bb + dd_inv[2][1][2]*cc ;
266 dd.mat[2][2] = 1.0f + dd_inv[2][2][0]*aa + dd_inv[2][2][1]*bb + dd_inv[2][2][2]*cc ;
267
268 ee = MAT_INV(dd) ;
269 v = MATVEC(ee,w) ;
270
271 *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
272 }
273
warper_bilinear_det(float ii,float jj,float kk)274 float warper_bilinear_det( float ii , float jj , float kk )
275 {
276 THD_fvec3 x,v,w ; THD_mat33 dd,ee ; int i,j ; float edet,adet,ddet , aa,bb,cc ;
277 static int first=1 ;
278
279 LOAD_FVEC3( x , ii,jj,kk ) ; v = VECMAT_VEC( ijk_to_xyz , x ) ;
280 UNLOAD_FVEC3( v , aa,bb,cc ) ; w = VECMAT_VEC( mv_for , v ) ;
281
282 dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
283 dd.mat[0][1] = dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
284 dd.mat[0][2] = dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
285 dd.mat[1][0] = dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
286 dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
287 dd.mat[1][2] = dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
288 dd.mat[2][0] = dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
289 dd.mat[2][1] = dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
290 dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
291
292 ddet = MAT_DET(dd) ;
293 if( first && fabs(ddet) < 0.001f ){
294 fprintf(stderr,"******* ddet=%g ii,jj,kk=%g %g %g aa,bb,cc=%g %g %g\n",
295 ddet,ii,jj,kk, aa,bb,cc ) ;
296 DUMP_MAT33("dd",dd) ;
297 DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
298 DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
299 first = 0 ;
300 }
301
302 ee = MAT_INV(dd) ; edet = MAT_DET(ee) ; v = MATVEC(ee,w) ;
303
304 for( i=0 ; i < 3 ; i++ )
305 for( j=0 ; j < 3 ; j++ )
306 dd.mat[i][j] = mv_for.mm.mat[i][j] - dd_for[i][0][j]*v.xyz[0]
307 - dd_for[i][1][j]*v.xyz[1]
308 - dd_for[i][2][j]*v.xyz[2] ;
309
310 adet = MAT_DET(dd) ;
311 return (adet*edet) ;
312 }
313
parset_bilinear(void)314 void parset_bilinear(void)
315 {
316 THD_mat33 ai ; THD_fvec3 df,di ; int i,j,k ;
317
318 parset_affine() ; /* sets up numerator matrices: mv_for and mv_inv */
319
320 /* load forward denominator 3-tensor */
321
322 dd_for[0][0][0] = parvec[12]; dd_for[0][0][1] = parvec[13]; dd_for[0][0][2] = parvec[14];
323 dd_for[0][1][0] = parvec[15]; dd_for[0][1][1] = parvec[16]; dd_for[0][1][2] = parvec[17];
324 dd_for[0][2][0] = parvec[18]; dd_for[0][2][1] = parvec[19]; dd_for[0][2][2] = parvec[20];
325 dd_for[1][0][0] = parvec[21]; dd_for[1][0][1] = parvec[22]; dd_for[1][0][2] = parvec[23];
326 dd_for[1][1][0] = parvec[24]; dd_for[1][1][1] = parvec[25]; dd_for[1][1][2] = parvec[26];
327 dd_for[1][2][0] = parvec[27]; dd_for[1][2][1] = parvec[28]; dd_for[1][2][2] = parvec[29];
328 dd_for[2][0][0] = parvec[30]; dd_for[2][0][1] = parvec[31]; dd_for[2][0][2] = parvec[32];
329 dd_for[2][1][0] = parvec[33]; dd_for[2][1][1] = parvec[34]; dd_for[2][1][2] = parvec[35];
330 dd_for[2][2][0] = parvec[36]; dd_for[2][2][1] = parvec[37]; dd_for[2][2][2] = parvec[38];
331
332 for( i=0 ; i < 3 ; i++ )
333 for( j=0 ; j < 3 ; j++ )
334 for( k=0 ; k < 3 ; k++ ) dd_for[i][j][k] *= dd_fac ; /* 18 Jul 2005 */
335
336 /* compute inverse denominator 3-tensor */
337
338 ai = mv_inv.mm ;
339 for( k=0 ; k < 3 ; k++ ){
340 for( j=0 ; j < 3 ; j++ ){
341 LOAD_FVEC3( df , -dd_for[0][k][j] , -dd_for[1][k][j] , -dd_for[2][k][j] ) ;
342 di = MATVEC( ai , df ) ;
343 UNLOAD_FVEC3( di , dd_inv[0][j][k] , dd_inv[1][j][k] , dd_inv[2][j][k] ) ;
344 }
345 }
346
347 #if 0
348 fprintf(stderr,"++++++++++++ parset_bilinear: dd_fac = %g\n",dd_fac) ;
349 for( k=0 ; k < 3 ; k++ ){
350 fprintf(stderr," +-------+ dd_for[][][%d] | dd_inv[][][%d]\n"
351 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n"
352 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n"
353 " %11.4g %11.4g %11.4g | %11.4g %11.4g %11.4g\n" ,
354 k , k ,
355 dd_for[0][0][k] , dd_for[0][1][k] , dd_for[0][2][k] ,
356 dd_inv[0][0][k] , dd_inv[0][1][k] , dd_inv[0][2][k] ,
357 dd_for[1][0][k] , dd_for[1][1][k] , dd_for[1][2][k] ,
358 dd_inv[1][0][k] , dd_inv[1][1][k] , dd_inv[1][2][k] ,
359 dd_for[2][0][k] , dd_for[2][1][k] , dd_for[2][2][k] ,
360 dd_inv[2][0][k] , dd_inv[2][1][k] , dd_inv[2][2][k] ) ;
361 }
362 #endif
363 }
364
365 /*--------------------------------------------------------------------------*/
366
main(int argc,char * argv[])367 int main( int argc , char * argv[] )
368 {
369 THD_3dim_dataset *inset=NULL, *outset=NULL, *baset=NULL, *wtset=NULL ;
370 MRI_warp3D_align_basis abas ;
371 char *prefix="warpdriven" ;
372 int warpdrive_code=-1 , nerr=0 , nx,ny,nz , nopt=1 ;
373 float dx,dy,dz , vp ;
374 int kim , nvals , kpar , np , nfree ;
375 MRI_IMAGE *qim , *tim , *fim ;
376 float clip_baset=0.0f , clip_inset=0.0f ;
377 char *W_1Dfile=NULL ; /* 04 Jan 2005 */
378 float **parsave=NULL ;
379 int output_float=0 ; /* 06 Jul 2005 */
380 char *base_idc=NULL , *wt_idc=NULL ;
381 int ctstart = NI_clock_time() ;
382 float i_xcm,i_ycm,i_zcm , b_xcm,b_ycm,b_zcm ; /* 26 Sep 2005 */
383 float sdif_before , sdif_after ; /* 28 Sep 2005 */
384 char *W_summfile=NULL ; FILE *summfp=NULL ;
385 int init_set=0 ; /* 06 Dec 2005 */
386 float init_val[12] ;
387 char *matrix_save_1D=NULL ; /* 25 Jul 2007 */
388 FILE *msfp=NULL ;
389 int null_output=0 ;
390
391 /*-- help? --*/
392
393 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
394 printf("\n"
395 "Usage: 3dWarpDrive [options] dataset\n"
396 "Warp a dataset to match another one (the base).\n"
397 "\n"
398 "* This program is a generalization of 3dvolreg. It tries to find\n"
399 " a spatial transformation that warps a given dataset to match an\n"
400 " input dataset (given by the -base option). It will be slow.\n"
401 "* Here, the spatical transformation is defined by a matrix; thus,\n"
402 " it is an affine warp.\n"
403 "* Program 3dAllineate can also compute such an affine transformation,\n"
404 " and it has more options for how the base and input (source) datasets\n"
405 " are to be matched. Thus, the usefulness of the older 3dWarpDrive\n"
406 " program is now limited. For future work, consider using 3dAllineate.\n"
407 "*****\n"
408 "***** For nonlinear spatial warping, see program 3dQwarp. *****\n"
409 "*****\n"
410 "\n"
411 " *** Also see the script align_epi_anat.py for a more general ***\n"
412 " ** alignment procedure, which does not require that the two **\n"
413 " ** datasets be defined on the same 3D grid. **\n"
414 " ** align_epi_anat.py uses program 3dAllineate. **\n"
415 "\n"
416 "--------------------------\n"
417 "Transform Defining Options: [exactly one of these must be used]\n"
418 "--------------------------\n"
419 " -shift_only = 3 parameters (shifts)\n"
420 " -shift_rotate = 6 parameters (shifts + angles)\n"
421 " -shift_rotate_scale = 9 parameters (shifts + angles + scale factors)\n"
422 " -affine_general = 12 parameters (3 shifts + 3x3 matrix)\n"
423 " -bilinear_general = 39 parameters (3 + 3x3 + 3x3x3)\n"
424 "\n"
425 " N.B.: At this time, the image intensity is NOT \n"
426 " adjusted for the Jacobian of the transformation.\n"
427 " N.B.: -bilinear_general is not yet implemented.\n"
428 " AND WILL NEVER BE.\n"
429 "\n"
430 "-------------\n"
431 "Other Options:\n"
432 "-------------\n"
433 " -linear }\n"
434 " -cubic } = Chooses spatial interpolation method.\n"
435 " -NN } = [default = linear; inaccurate but fast]\n"
436 " -quintic } [for accuracy, try '-cubic -final quintic']\n"
437 "\n"
438 " -base bbb = Load dataset 'bbb' as the base to which the\n"
439 " input dataset will be matched.\n"
440 " [This is a mandatory option]\n"
441 "\n"
442 " -verb = Print out lots of information along the way.\n"
443 " -prefix ppp = Sets the prefix of the output dataset.\n"
444 " If 'ppp' is 'NULL', no output dataset is written.\n"
445 " -input ddd = You can put the input dataset anywhere in the\n"
446 " command line option list by using the '-input'\n"
447 " option, instead of always putting it last.\n"
448 " -summ sss = Save summary of calculations into text file 'sss'.\n"
449 " (N.B.: If 'sss' is '-', summary goes to stdout.)\n"
450 "\n"
451 "-----------------\n"
452 "Technical Options:\n"
453 "-----------------\n"
454 " -maxite m = Allow up to 'm' iterations for convergence.\n"
455 " -delta d = Distance, in voxel size, used to compute\n"
456 " image derivatives using finite differences.\n"
457 " [Default=1.0]\n"
458 " -weight wset = Set the weighting applied to each voxel\n"
459 " proportional to the brick specified here.\n"
460 " [Default=computed by program from base]\n"
461 " -thresh t = Set the convergence parameter to be RMS 't' voxels\n"
462 " movement between iterations. [Default=0.03]\n"
463 " -twopass = Do the parameter estimation in two passes,\n"
464 " coarse-but-fast first, then fine-but-slow second\n"
465 " (much like the same option in program 3dvolreg).\n"
466 " This is useful if large-ish warping is needed to\n"
467 " align the volumes.\n"
468 " -final 'mode' = Set the final warp to be interpolated using 'mode'\n"
469 " instead of the spatial interpolation method used\n"
470 " to find the warp parameters.\n"
471 " -parfix n v = Fix the n'th parameter of the warp model to\n"
472 " the value 'v'. More than one -parfix option\n"
473 " can be used, to fix multiple parameters.\n"
474 " -1Dfile ename = Write out the warping parameters to the file\n"
475 " named 'ename'. Each sub-brick of the input\n"
476 " dataset gets one line in this file. Each\n"
477 " parameter in the model gets one column.\n"
478 " -float = Write output dataset in float format, even if\n"
479 " input dataset is short or byte.\n"
480 " -coarserot = Initialize shift+rotation parameters by a\n"
481 " brute force coarse search, as in the similar\n"
482 " 3dvolreg option.\n"
483 "\n"
484 " -1Dmatrix_save ff = Save base-to-input transformation matrices\n"
485 " in file 'ff' (1 row per sub-brick in the input\n"
486 " dataset). If 'ff' does NOT end in '.1D', then\n"
487 " the program will append '.aff12.1D' to 'ff' to\n"
488 " make the output filename.\n"
489 " *N.B.: This matrix is the coordinate transformation from base\n"
490 " to input DICOM coordinates. To get the inverse matrix\n"
491 " (input-to-base), use the cat_matvec program, as in\n"
492 " cat_matvec fred.aff12.1D -I\n"
493 "\n"
494 "----------------------\n"
495 "AFFINE TRANSFORMATIONS:\n"
496 "----------------------\n"
497 "The options below control how the affine tranformations\n"
498 "(-shift_rotate, -shift_rotate_scale, -affine_general)\n"
499 "are structured in terms of 3x3 matrices:\n"
500 "\n"
501 " -SDU or -SUD }= Set the order of the matrix multiplication\n"
502 " -DSU or -DUS }= for the affine transformations:\n"
503 " -USD or -UDS }= S = triangular shear (params #10-12)\n"
504 " D = diagonal scaling matrix (params #7-9)\n"
505 " U = rotation matrix (params #4-6)\n"
506 " Default order is '-SDU', which means that\n"
507 " the U matrix is applied first, then the\n"
508 " D matrix, then the S matrix.\n"
509 "\n"
510 " -Supper }= Set the S matrix to be upper or lower\n"
511 " -Slower }= triangular [Default=lower triangular]\n"
512 "\n"
513 " -ashift OR }= Apply the shift parameters (#1-3) after OR\n"
514 " -bshift }= before the matrix transformation. [Default=after]\n"
515 "\n"
516 "The matrices are specified in DICOM-ordered (x=-R+L,y=-A+P,z=-I+S)\n"
517 "coordinates as:\n"
518 "\n"
519 " [U] = [Rotate_y(param#6)] [Rotate_x(param#5)] [Rotate_z(param #4)]\n"
520 " (angles are in degrees)\n"
521 "\n"
522 " [D] = diag( param#7 , param#8 , param#9 )\n"
523 "\n"
524 " [ 1 0 0 ] [ 1 param#10 param#11 ]\n"
525 " [S] = [ param#10 1 0 ] OR [ 0 1 param#12 ]\n"
526 " [ param#11 param#12 1 ] [ 0 0 1 ]\n"
527 "\n"
528 " For example, the default (-SDU/-ashift/-Slower) has the warp\n"
529 " specified as [x]_warped = [S] [D] [U] [x]_in + [shift].\n"
530 " The shift vector comprises parameters #1, #2, and #3.\n"
531 "\n"
532 " The goal of the program is to find the warp parameters such that\n"
533 " I([x]_warped) = s * J([x]_in)\n"
534 " as closely as possible in a weighted least squares sense, where\n"
535 " 's' is a scaling factor (an extra, invisible, parameter), J(x)\n"
536 " is the base image, I(x) is the input image, and the weight image\n"
537 " is a blurred copy of J(x).\n"
538 "\n"
539 " Using '-parfix', you can specify that some of these parameters\n"
540 " are fixed. For example, '-shift_rotate_scale' is equivalent\n"
541 " '-affine_general -parfix 10 0 -parfix 11 0 -parfix 12 0'.\n"
542 " Don't attempt to use the '-parfix' option unless you understand\n"
543 " this example!\n"
544 "\n"
545 "-------------------------\n"
546 " RWCox - November 2004\n"
547 "-------------------------\n"
548 ) ;
549 PRINT_COMPILE_DATE ; exit(0) ;
550 }
551
552 /*-- startup mechanics --*/
553
554 #ifdef USING_MCW_MALLOC
555 enable_mcw_malloc() ;
556 #endif
557
558 mainENTRY("3dWarpDrive main"); machdep(); AFNI_logger("3dWarpDrive",argc,argv);
559 PRINT_VERSION("3dWarpDrive") ; AUTHOR("RW Cox") ;
560 THD_check_AFNI_version("3dWarpDrive") ;
561
562 /* initialize parameters of the alignment basis struct */
563
564 abas.nparam = 0 ;
565 abas.param = NULL ;
566 abas.scale_init = 1.0f ;
567 abas.delfac = 1.0f ;
568 abas.tolfac = 0.03f ;
569 abas.twoblur = 0.0f ;
570 abas.regmode = MRI_LINEAR ;
571 abas.regfinal = -1 ;
572 abas.verb = 0 ;
573 abas.max_iter = 0 ;
574 abas.wtproc = 1 ;
575 abas.imbase = NULL ;
576 abas.imwt = NULL ;
577 abas.vwfor = ijk_warp_for ;
578 abas.vwinv = ijk_warp_inv ;
579 abas.vwset = load_parvec ;
580 abas.vwdet = NULL ;
581
582 abas.xedge = abas.yedge = abas.zedge = -1 ;
583 abas.imww = abas.imap = abas.imps = abas.imsk = NULL ;
584 abas.imps_blur = NULL ;
585
586 nparfix = 0 ;
587
588 /*-- command line options --*/
589
590 while( nopt < argc && argv[nopt][0] == '-' ){
591
592 /*-----*/
593
594 if( strcmp(argv[nopt],"-summ") == 0 ){ /* 28 Sep 2005 */
595 if( ++nopt >= argc )
596 ERROR_exit("Need 1 parameter afer -summ!\n");
597 W_summfile = strdup( argv[nopt] ) ;
598 if( !THD_filename_ok(W_summfile) )
599 ERROR_exit("Name after -summ has bad characters!\n") ;
600 nopt++ ; continue ;
601 }
602
603 /*-----*/
604
605 if( strcmp(argv[nopt],"-float") == 0 ){ /* 06 Jul 2005 */
606 output_float = 1 ; nopt++ ; continue ;
607 }
608
609 /*-----*/
610
611 if( strcmp(argv[nopt],"-coarserot") == 0 ){ /* 05 Dec 2005 */
612 if( VL_coarse_rot ) INFO_message("-coarserot is already on") ;
613 VL_coarse_rot = 1 ; nopt++ ; continue ;
614 }
615
616 /*-----*/
617
618 if( strcmp(argv[nopt],"-1Dmatrix_save") == 0 ){
619 if( ++nopt >= argc )
620 ERROR_exit("Need 1 parameter afer -1Dmatrix_save!\n");
621 if( !THD_filename_ok(argv[nopt]) )
622 ERROR_exit("badly formed filename: %s '%s'",argv[nopt-1],argv[nopt]) ;
623 if( STRING_HAS_SUFFIX(argv[nopt],".1D") ){
624 matrix_save_1D = strdup(argv[nopt]) ;
625 } else {
626 matrix_save_1D = calloc(sizeof(char),strlen(argv[nopt])+16) ;
627 strcpy(matrix_save_1D,argv[nopt]); strcat(matrix_save_1D,".aff12.1D");
628 }
629 nopt++ ; continue ;
630 }
631
632 /*-----*/
633
634 if( strcmp(argv[nopt],"-1Dfile") == 0 ){
635 if( ++nopt >= argc )
636 ERROR_exit("Need 1 parameter afer -1Dfile!\n");
637 W_1Dfile = strdup( argv[nopt] ) ;
638 if( !THD_filename_ok(W_1Dfile) )
639 ERROR_exit("Name after -1Dfile has bad characters!\n") ;
640 nopt++ ; continue ;
641 }
642
643 /*-----*/
644
645 if( strcmp(argv[nopt],"-twopass") == 0 ){
646 float bbb = AFNI_numenv("AFNI_WARPDRIVE_TWOBLUR") ;
647 abas.twoblur = (bbb==0.0f) ? 3.0f : bbb ;
648 nopt++ ; continue ;
649 }
650
651 /*-----*/
652
653 if( strcmp(argv[nopt],"-SDU") == 0 ){
654 matorder = MATORDER_SDU ; nopt++ ; continue ;
655 }
656 if( strcmp(argv[nopt],"-SUD") == 0 ){
657 matorder = MATORDER_SUD ; nopt++ ; continue ;
658 }
659 if( strcmp(argv[nopt],"-DSU") == 0 ){
660 matorder = MATORDER_DSU ; nopt++ ; continue ;
661 }
662 if( strcmp(argv[nopt],"-DUS") == 0 ){
663 matorder = MATORDER_DUS ; nopt++ ; continue ;
664 }
665 if( strcmp(argv[nopt],"-USD") == 0 ){
666 matorder = MATORDER_USD ; nopt++ ; continue ;
667 }
668 if( strcmp(argv[nopt],"-UDS") == 0 ){
669 matorder = MATORDER_UDS ; nopt++ ; continue ;
670 }
671 if( strcmp(argv[nopt],"-ashift") == 0 ){
672 dcode = DELTA_AFTER ; nopt++ ; continue ;
673 }
674 if( strcmp(argv[nopt],"-bshift") == 0 ){
675 dcode = DELTA_BEFORE ; nopt++ ; continue ;
676 }
677 if( strcmp(argv[nopt],"-Slower") == 0 ){
678 smat = SMAT_LOWER ; nopt++ ; continue ;
679 }
680 if( strcmp(argv[nopt],"-Supper") == 0 ){
681 smat = SMAT_UPPER ; nopt++ ; continue ;
682 }
683
684 /*-----*/
685
686 if( strcmp(argv[nopt],"-final") == 0 ){
687 char *str ;
688
689 if( ++nopt >= argc )
690 ERROR_exit("Need 1 parameter afer -final!\n") ;
691 str = argv[nopt] ; if( *str == '-' ) str++ ;
692
693 if( strcmp(str,"cubic") == 0 ) abas.regfinal = MRI_CUBIC ;
694 else if( strcmp(str,"quintic") == 0 ) abas.regfinal = MRI_QUINTIC ;
695 else if( strcmp(str,"linear") == 0 ) abas.regfinal = MRI_LINEAR ;
696 else if( strcmp(str,"NN") == 0 ) abas.regfinal = MRI_NN ;
697 else
698 ERROR_exit("Illegal mode after -final\n");
699 nopt++ ; continue ;
700 }
701
702 /*-----*/
703
704 if( strcmp(argv[nopt],"-parfix") == 0 ){
705 int np , ip ; float vp ;
706 if( ++nopt >= argc-1 )
707 ERROR_exit("Need 2 parameters afer -parfix!\n");
708 np = strtol( argv[nopt] , NULL , 10 ) ; nopt++ ;
709 vp = strtod( argv[nopt] , NULL ) ;
710 if( np <= 0 || np > MAXPAR )
711 ERROR_exit("Param #%d after -parfix is illegal!\n",np) ;
712 for( ip=0 ; ip < nparfix ; ip++ ){
713 if( parfix[ip].np == np ){
714 WARNING_message("Ignoring later -parfix option for param #%d\n",ip);
715 break ;
716 }
717 }
718 if( ip == nparfix ) nparfix++ ;
719 parfix[ip].np = np ; parfix[ip].vp = vp ;
720 nopt++ ; continue ;
721 }
722
723 /*-----*/
724
725 if( strcmp(argv[nopt],"-shift_only") == 0 ){
726 warpdrive_code = WARPDRIVE_SHIFT ; nopt++ ; continue ;
727 }
728
729 if( strcmp(argv[nopt],"-shift_rotate") == 0 ){
730 warpdrive_code = WARPDRIVE_ROTATE ; nopt++ ; continue ;
731 }
732
733 if( strcmp(argv[nopt],"-shift_rotate_scale") == 0 ){
734 warpdrive_code = WARPDRIVE_SCALE ; nopt++ ; continue ;
735 }
736
737 if( strcmp(argv[nopt],"-affine_general") == 0 ){
738 warpdrive_code = WARPDRIVE_AFFINE ; nopt++ ; continue ;
739 }
740
741 if( strcmp(argv[nopt],"-bilinear_general") == 0 ){ /* not implemented */
742 #if 0
743 ERROR_exit("3dWarpDrive -bilinear_general NOT IMPLEMENTED!\n");
744 #else
745 WARNING_message("-bilinear_general transformation is experimental!") ;
746 #endif
747 warpdrive_code = WARPDRIVE_BILINEAR ; nopt++ ; continue ;
748 }
749
750 /*-----*/
751
752 if( strcmp(argv[nopt],"-NN") == 0 ){
753 abas.regmode = MRI_NN ; nopt++ ; continue ;
754 }
755 if( strcmp(argv[nopt],"-linear") == 0 ){
756 abas.regmode = MRI_LINEAR ; nopt++ ; continue ;
757 }
758 if( strcmp(argv[nopt],"-cubic") == 0 ){
759 abas.regmode = MRI_CUBIC ; nopt++ ; continue ;
760 }
761 if( strcmp(argv[nopt],"-quintic") == 0 ){
762 abas.regmode = MRI_QUINTIC ; nopt++ ; continue ;
763 }
764
765 /*-----*/
766
767 if( strcmp(argv[nopt],"-prefix") == 0 ){
768 if( ++nopt >= argc )
769 ERROR_exit("Need an argument after -prefix!\n");
770 if( !THD_filename_ok(argv[nopt]) )
771 ERROR_exit("-prefix argument is invalid!\n");
772 prefix = argv[nopt] ; if( strcmp(prefix,"NULL") == 0 ) null_output = 1 ;
773 nopt++ ; continue ;
774 }
775
776 /*-----*/
777
778 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
779 abas.verb++ ; nopt++ ; continue ;
780 }
781
782 /*-----*/
783
784 if( strcmp(argv[nopt],"-base") == 0 ){
785 if( ++nopt >= argc )
786 ERROR_exit("Need an argument after -base!\n");
787 baset = THD_open_dataset( argv[nopt] ) ;
788 if( baset == NULL )
789 ERROR_exit("Can't open -base dataset %s\n",argv[nopt]);
790 if( DSET_NVALS(baset) > 1 )
791 WARNING_message(
792 "-base dataset %s has %d sub-bricks; will only use #0\n",
793 argv[nopt],DSET_NVALS(baset) ) ;
794 nopt++ ; continue ;
795 }
796
797 /*-----*/
798
799 if( strcmp(argv[nopt],"-weight") == 0 ){
800 if( ++nopt >= argc )
801 ERROR_exit("Need an argument after -weight!\n");
802 if( wtset != NULL ){
803 WARNING_message("2nd -weight option replaces 1st") ;
804 DSET_delete(wtset) ;
805 }
806 wtset = THD_open_dataset( argv[nopt] ) ;
807 if( wtset == NULL )
808 ERROR_exit("Can't open -weight dataset %s\n",argv[nopt]);
809 if( DSET_NVALS(wtset) > 1 )
810 WARNING_message(
811 "-weight %s has %d sub-bricks; will only use #0\n",
812 argv[nopt],DSET_NVALS(wtset) ) ;
813 nopt++ ; continue ;
814 }
815
816 /*-----*/
817
818 if( strcmp(argv[nopt],"-input") == 0 || strcmp(argv[nopt],"-source") == 0 ){
819 if( ++nopt >= argc )
820 ERROR_exit("Need an argument after %s!",argv[nopt-1]);
821 inset = THD_open_dataset( argv[nopt] ) ;
822 if( inset == NULL )
823 ERROR_exit("Can't open %s dataset %s\n",argv[nopt-1],argv[nopt]);
824 nopt++ ; continue ;
825 }
826
827 /*-----*/
828
829 if( strcmp(argv[nopt],"-maxite") == 0 ){
830 int ival ;
831 if( ++nopt >= argc )
832 ERROR_exit("Need an argument after -maxite!\n");
833 ival = strtol( argv[nopt] , NULL , 10 ) ;
834 if( ival > 1 ) abas.max_iter = ival ;
835 nopt++ ; continue ;
836 }
837
838 /*-----*/
839
840 if( strcmp(argv[nopt],"-delta") == 0 ){
841 float val ;
842 if( ++nopt >= argc )
843 ERROR_exit("Need an argument after -delta!\n");
844 val = strtod( argv[nopt] , NULL ) ;
845 if( val > 0.0499 && val < 49.99 ) abas.delfac = val ;
846 nopt++ ; continue ;
847 }
848
849 /*-----*/
850
851 if( strcmp(argv[nopt],"-thresh") == 0 ){
852 float val ;
853 if( ++nopt >= argc )
854 ERROR_exit("Need an argument after -thresh!\n");
855 val = strtod( argv[nopt] , NULL ) ;
856 if( val > 0.001 && val < 3.01 ) abas.tolfac = val ;
857 nopt++ ; continue ;
858 }
859
860 /*-----*/
861
862 ERROR_exit("Unknown option %s\n",argv[nopt]) ;
863
864 } /*--- end of loop over command line options ---*/
865
866 /*-- 1 remaining argument should be a dataset --*/
867
868 if( inset == NULL && nopt != argc-1 )
869 ERROR_exit("Command line should have exactly 1 dataset!\n"
870 "** Whereas there seems to be %d of them!\n",
871 argc-nopt ) ;
872
873 /*-- input dataset header --*/
874
875 if( inset == NULL ){
876 inset = THD_open_dataset( argv[nopt] ) ;
877 if( !ISVALID_DSET(inset) )
878 ERROR_exit("Can't open dataset %s\n",argv[nopt]);
879 }
880
881 if( baset == NULL ){
882 ERROR_message("Need to specify a base dataset!\n") ;
883 nerr++ ;
884 }
885
886 nx = DSET_NX(inset) ; ny = DSET_NY(inset) ; nz = DSET_NZ(inset) ;
887 dx = DSET_DX(inset) ; dy = DSET_DY(inset) ; dz = DSET_DZ(inset) ;
888
889 if( DSET_NX(baset) != nx || DSET_NY(baset) != ny || DSET_NZ(baset) != nz ){
890 ERROR_message("base and input grid dimensions don't match!\n") ;
891 ERROR_message("base is %d X %d X %d voxels\n",
892 DSET_NX(baset),DSET_NY(baset),DSET_NZ(baset) ) ;
893 ERROR_exit ("input is %d X %d X %d voxels\n",nx,ny,nz) ;
894 }
895
896 /*- load datasets from disk; if can't do so, fatal errors all around -*/
897
898 if( abas.verb ) INFO_message("Loading datasets\n") ;
899
900 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
901 nvals = DSET_NVALS(inset) ;
902 if( nvals == 1 ){
903 clip_inset = THD_cliplevel( DSET_BRICK(inset,0) , 0.0 ) ;
904 if( DSET_BRICK_FACTOR(inset,0) > 0.0f )
905 clip_inset *= DSET_BRICK_FACTOR(inset,0) ;
906 mri_get_cmass_3D( DSET_BRICK(inset,0) , &i_xcm,&i_ycm,&i_zcm ) ;
907 } else {
908 qim = THD_median_brick( inset ) ;
909 clip_inset = THD_cliplevel( qim , 0.0 ) ;
910 mri_get_cmass_3D( qim , &i_xcm,&i_ycm,&i_zcm ) ;
911 mri_free(qim) ;
912 }
913
914 DSET_load(baset) ; CHECK_LOAD_ERROR(baset) ;
915 mri_get_cmass_3D( DSET_BRICK(baset,0) , &b_xcm,&b_ycm,&b_zcm ) ;
916 abas.imbase = mri_scale_to_float( DSET_BRICK_FACTOR(baset,0) ,
917 DSET_BRICK(baset,0) ) ;
918 clip_baset = THD_cliplevel( abas.imbase , 0.0 ) ;
919 base_idc = strdup(baset->idcode.str) ;
920 DSET_unload(baset) ;
921
922 if( wtset != NULL ){
923 DSET_load(wtset) ; CHECK_LOAD_ERROR(wtset) ;
924 abas.imwt = mri_scale_to_float( DSET_BRICK_FACTOR(wtset,0) ,
925 DSET_BRICK(wtset,0) ) ;
926 wt_idc = strdup(wtset->idcode.str) ;
927 DSET_unload(wtset) ;
928 }
929
930 if( abas.verb ) INFO_message("Checking inputs\n") ;
931
932 /*-- parameterize the warp model --*/
933
934 /*! Macro to add a parameter to the warp3D model.
935 - nm = name of parameter
936 - bb = min value allowed
937 - tt = max value allowed
938 - id = value for identity warp
939 - dd = delta to use for stepsize
940 - ll = tolerance for convergence test */
941
942 #define ADDPAR(nm,bb,tt,id,dd,ll) \
943 do{ int p=abas.nparam ; \
944 abas.param = (MRI_warp3D_param_def *) realloc( \
945 (void *)abas.param , \
946 sizeof(MRI_warp3D_param_def)*(p+1) ) ; \
947 abas.param[p].min = (bb) ; abas.param[p].max = (tt) ; \
948 abas.param[p].delta = (dd) ; abas.param[p].toler = (ll) ; \
949 abas.param[p].ident = abas.param[p].val_init = (id) ; \
950 strcpy( abas.param[p].name , (nm) ) ; \
951 abas.param[p].fixed = 0 ; \
952 abas.nparam = p+1 ; \
953 } while(0)
954
955 nerr = 0 ;
956 if( warpdrive_code <= 0 ){
957 ERROR_message("Need a transform-specifying option!\n");
958 nerr++ ;
959 } else if( WARPDRIVE_IS_AFFINE(warpdrive_code) ) {
960
961 char *lab09, *lab10, *lab11 ;
962 float xd,yd,zd , xl,yl,zl , xxx ;
963
964 xd = 3.33*fabs( dx * (i_xcm-b_xcm) ) ; xl = 0.333 * fabs( dx*nx ) ;
965 yd = 3.33*fabs( dy * (i_ycm-b_ycm) ) ; yl = 0.333 * fabs( dy*ny ) ;
966 zd = 3.33*fabs( dz * (i_zcm-b_zcm) ) ; zl = 0.333 * fabs( dz*nz ) ;
967 xd = MAX(xd,xl) ; yd = MAX(yd,yl) ; zd = MAX(zd,zl) ;
968 xxx = sqrt( xd*xd + yd*yd + zd*zd ) ;
969
970 /* add all 12 parameters (may ignore some, later) */
971
972 ADDPAR( "x-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
973 ADDPAR( "y-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
974 ADDPAR( "z-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
975
976 ADDPAR( "z-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ; /* degrees */
977 ADDPAR( "x-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ;
978 ADDPAR( "y-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ;
979
980 ADDPAR( "x-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ; /* identity */
981 ADDPAR( "y-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ; /* == 1.0 */
982 ADDPAR( "z-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
983
984 switch( smat ){
985 default:
986 case SMAT_LOWER:
987 lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
988 break ;
989
990 case SMAT_UPPER:
991 lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
992 break ;
993 }
994 ADDPAR( lab09 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
995 ADDPAR( lab10 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
996 ADDPAR( lab11 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
997
998 /* initialize transform parameter vector */
999
1000 for( kpar=0 ; kpar < 12 ; kpar++ )
1001 parvec[kpar] = abas.param[kpar].ident ;
1002
1003 /* initialize transformation function pointers */
1004
1005 if( warpdrive_code == WARPDRIVE_SHIFT ){
1006 warp_parset = parset_shift ;
1007 warp_for = warper_shift_for ;
1008 warp_inv = warper_shift_inv ;
1009 } else {
1010 warp_parset = parset_affine ;
1011 warp_for = warper_affine_for ;
1012 warp_inv = warper_affine_inv ;
1013 }
1014
1015 /* how many parameters to actually pay attention to */
1016
1017 switch( warpdrive_code ){
1018 case WARPDRIVE_SHIFT: abas.nparam = 3 ; break ;
1019 case WARPDRIVE_ROTATE: abas.nparam = 6 ; break ;
1020 case WARPDRIVE_SCALE: abas.nparam = 9 ; break ;
1021 case WARPDRIVE_AFFINE: abas.nparam = 12 ; break ;
1022 }
1023
1024 } else if( warpdrive_code == WARPDRIVE_BILINEAR ){
1025
1026 char *lab09, *lab10, *lab11 , labxx[16] ;
1027 float xr,yr,zr,rr ;
1028 float xd,yd,zd , xl,yl,zl , xxx ;
1029
1030 xd = 3.33*fabs( dx * (i_xcm-b_xcm) ) ; xl = 0.333 * fabs( dx*nx ) ;
1031 yd = 3.33*fabs( dy * (i_ycm-b_ycm) ) ; yl = 0.333 * fabs( dy*ny ) ;
1032 zd = 3.33*fabs( dz * (i_zcm-b_zcm) ) ; zl = 0.333 * fabs( dz*nz ) ;
1033 xd = MAX(xd,xl) ; yd = MAX(yd,yl) ; zd = MAX(zd,zl) ;
1034 xxx = sqrt( xd*xd + yd*yd + zd*zd ) ;
1035
1036 /* add all 39 parameters (may ignore some, later) */
1037
1038 ADDPAR( "x-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
1039 ADDPAR( "y-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
1040 ADDPAR( "z-shift" , -xxx , xxx , 0.0 , 0.0 , 0.0 ) ;
1041
1042 ADDPAR( "z-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ;
1043 ADDPAR( "x-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ;
1044 ADDPAR( "y-angle" , -120.0 , 120.0 , 0.0 , 0.0 , 0.0 ) ;
1045
1046 ADDPAR( "x-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
1047 ADDPAR( "y-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
1048 ADDPAR( "z-scale" , 0.618 , 1.618 , 1.0 , 0.0 , 0.0 ) ;
1049
1050 switch( smat ){
1051 default:
1052 case SMAT_LOWER:
1053 lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
1054 break ;
1055
1056 case SMAT_UPPER:
1057 lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
1058 break ;
1059 }
1060 ADDPAR( lab09 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
1061 ADDPAR( lab10 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
1062 ADDPAR( lab11 , -0.1666 , 0.1666 , 0.0 , 0.0 , 0.0 ) ;
1063
1064 xr = 0.5f*fabs(dx)*nx ; yr = 0.5f*fabs(dy)*ny ; zr = 0.5f*fabs(dz)*nz ;
1065 rr = MAX(xr,yr) ; rr = MAX(rr,zr) ; dd_fac = 1.0f / rr ;
1066 for( kpar=12 ; kpar < 39 ; kpar++ ){
1067 sprintf(labxx,"blin_%02d",kpar+1) ;
1068 ADDPAR( labxx , -0.1234 , 0.1234 , 0.0 , 0.0 , 0.0 ) ;
1069 }
1070
1071 /* initialize transform parameter vector */
1072
1073 for( kpar=0 ; kpar < 39 ; kpar++ )
1074 parvec[kpar] = abas.param[kpar].ident ;
1075
1076 abas.vwdet = warper_bilinear_det ;
1077
1078 /* initialize transformation function pointers */
1079
1080 warp_parset = parset_bilinear ;
1081 warp_for = warper_bilinear_for ;
1082 warp_inv = warper_bilinear_inv ;
1083
1084 /* how many parameters to actually pay attention to */
1085
1086 abas.nparam = 39 ;
1087
1088 } else {
1089 ERROR_message("Unimplemented transform model!\n") ;
1090 nerr++ ;
1091 }
1092
1093 /* Deal with -parfix options; nfree will be number of free parameters */
1094
1095 nfree = abas.nparam ;
1096 for( kpar=0 ; kpar < nparfix ; kpar++ ){
1097 np = parfix[kpar].np - 1 ; vp = parfix[kpar].vp ;
1098 if( np >= 0 && np < abas.nparam ){
1099 if( vp >= abas.param[np].min && vp <= abas.param[np].max ){
1100 abas.param[np].fixed = 1 ;
1101 abas.param[np].val_fixed = vp ;
1102 nfree -- ;
1103 } else { /* bad value */
1104 ERROR_message("-parfix for param #%d has illegal value!\n",np+1) ;
1105 nerr++ ;
1106 }
1107 } else { /* bad index */
1108 WARNING_message("-parfix for param #%d is out of range 1..%d\n",
1109 np+1 , abas.nparam+1 ) ;
1110 }
1111 }
1112 if( nfree <= 0 ){
1113 ERROR_message("No free parameters in transform model!\n") ;
1114 nerr++ ;
1115 }
1116
1117 /*-- default number of iterations allowed --*/
1118
1119 if( abas.max_iter <= 0 ) abas.max_iter = 11*nfree+5 ;
1120
1121 /*-- other checks for good set of inputs --*/
1122
1123 if( nerr ) exit(1) ; /** bad user!! **/
1124
1125 if( abas.verb ) INFO_message("Creating empty output dataset\n") ;
1126
1127 outset = EDIT_empty_copy( inset ) ;
1128
1129 EDIT_dset_items( outset , ADN_prefix,prefix , ADN_none ) ;
1130
1131 if( output_float ){
1132 EDIT_dset_items( outset , ADN_datum_all , MRI_float , ADN_none ) ;
1133 for( kim=0 ; kim < nvals ; kim++ )
1134 EDIT_BRICK_FACTOR( outset , kim , 0.0 ) ;
1135 }
1136
1137 if( THD_deathcon() && THD_is_file( DSET_HEADNAME(outset) ) )
1138 ERROR_exit("Output file %s already exists!\n",DSET_HEADNAME(outset) ) ;
1139
1140 tross_Copy_History( inset , outset ) ;
1141 tross_Make_History( "3dWarpDrive" , argc,argv , outset ) ;
1142
1143 outset->daxes->xxorg = inset->daxes->xxorg ;
1144 outset->daxes->yyorg = inset->daxes->yyorg ;
1145 outset->daxes->zzorg = inset->daxes->zzorg ;
1146
1147 /*-- more checks --*/
1148
1149 /* the following aren't fatal errors, but merit a warning slap in the face */
1150 /* [modified to be less stringent [11 Feb 2018] -- RWC */
1151
1152 if( FLDIF(DSET_DX(baset),dx) > 0.002f ||
1153 FLDIF(DSET_DY(baset),dy) > 0.002f || FLDIF(DSET_DZ(baset),dz) > 0.002f ){
1154 WARNING_message("base and input grid spacings don't match!\n"
1155 /*"If this is not by design or if registration\n"
1156 "fails, use 3dresample with -master option:\n"
1157 " 3dresample -master %s \\\n"
1158 " -input %s \\\n"
1159 " -prefix RS\n"
1160 "to sesample the input data and rerun the\n"
1161 "registration with RS as the input.\n", DSET_BRIKNAME(baset), DSET_BRIKNAME(inset) */) ;
1162 WARNING_message("base grid = %.6f X %.6f X %.6f mm\n",
1163 DSET_DX(baset),DSET_DY(baset),DSET_DZ(baset) ) ;
1164 WARNING_message("input grid = %.6f X %.6f X %.6f mm\n",
1165 DSET_DX(inset),DSET_DY(inset),DSET_DZ(inset) ) ;
1166 }
1167
1168 if( FLDIF(DSET_XORG(baset),DSET_XORG(inset)) > 0.02f ||
1169 FLDIF(DSET_YORG(baset),DSET_YORG(inset)) > 0.02f ||
1170 FLDIF(DSET_ZORG(baset),DSET_ZORG(inset)) > 0.02f ){
1171 WARNING_message("base and input grid offsets don't match!\n"
1172 /*"If this is not by design or if registration\n"
1173 "fails, use 3dresample with -master option:\n"
1174 " 3dresample -master %s\\\n"
1175 " -input %s \\\n"
1176 " -prefix RS\n"
1177 "to sesample the input data and rerun the\n"
1178 "registration with RS as the input.\n", DSET_BRIKNAME(baset), DSET_BRIKNAME(inset) */) ;
1179 WARNING_message("base offsets = %.6f X %.6f X %.6f mm\n",
1180 DSET_XORG(baset),DSET_YORG(baset),DSET_ZORG(baset) ) ;
1181 WARNING_message("input offsets = %.6f X %.6f X %.6f mm\n",
1182 DSET_XORG(inset),DSET_YORG(inset),DSET_ZORG(inset) ) ;
1183 }
1184
1185 if( baset->daxes->xxorient != inset->daxes->xxorient ||
1186 baset->daxes->yyorient != inset->daxes->yyorient ||
1187 baset->daxes->zzorient != inset->daxes->zzorient ){
1188 WARNING_message("base and input orientations don't match!\n"
1189 /*"If this is not by design or if registration\n"
1190 "fails, use 3dresample with -master option:\n"
1191 " 3dresample -master %s \\\n"
1192 " -input %s \\\n"
1193 " -prefix RS\n"
1194 "to sesample the input data and rerun the\n"
1195 "registration with RS as the input.\n", DSET_BRIKNAME(baset), DSET_BRIKNAME(inset) */) ;
1196 WARNING_message("base = %s X %s X %s\n",
1197 ORIENT_shortstr[baset->daxes->xxorient] ,
1198 ORIENT_shortstr[baset->daxes->yyorient] ,
1199 ORIENT_shortstr[baset->daxes->zzorient] ) ;
1200 WARNING_message("input = %s X %s X %s\n",
1201 ORIENT_shortstr[inset->daxes->xxorient] ,
1202 ORIENT_shortstr[inset->daxes->yyorient] ,
1203 ORIENT_shortstr[inset->daxes->zzorient] ) ;
1204 }
1205
1206 /*- however, this is a fatal error -*/
1207
1208 if( wtset != NULL &&
1209 (DSET_NX(wtset) != nx || DSET_NY(wtset) != ny || DSET_NZ(wtset) != nz) ){
1210 ERROR_message("weight and input grid dimensions don't match!\n") ;
1211 ERROR_message("weight is %d X %d X %d voxels\n",
1212 DSET_NX(wtset),DSET_NY(wtset),DSET_NZ(wtset) ) ;
1213 ERROR_exit ("input is %d X %d X %d voxels\n",nx,ny,nz) ;
1214 }
1215
1216 /*-- set up (x,y,z) <-> (i,j,k) transformations ---*/
1217
1218 { THD_vecmat ijk_to_inset_xyz , xyz_to_dicom ;
1219
1220 LOAD_DIAG_MAT( ijk_to_inset_xyz.mm ,
1221 inset->daxes->xxdel ,
1222 inset->daxes->yydel , inset->daxes->zzdel );
1223
1224 if( warpdrive_code == WARPDRIVE_BILINEAR ){
1225 /* define (x,y,z)=(0,0,0) at mid-point of dataset 3D array */
1226
1227 LOAD_FVEC3 ( ijk_to_inset_xyz.vv ,
1228 -0.5*(nx-1)*inset->daxes->xxdel ,
1229 -0.5*(ny-1)*inset->daxes->yydel ,
1230 -0.5*(nz-1)*inset->daxes->zzdel ) ;
1231
1232 } else {
1233 /* define (x,y,z) based strictly on dataset coords */
1234
1235 LOAD_FVEC3 ( ijk_to_inset_xyz.vv ,
1236 DSET_XORG(inset) , DSET_YORG(inset), DSET_ZORG(inset) ) ;
1237 }
1238
1239 xyz_to_dicom.mm = inset->daxes->to_dicomm ;
1240 LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
1241
1242 ijk_to_xyz = MUL_VECMAT( xyz_to_dicom , ijk_to_inset_xyz ) ;
1243 xyz_to_ijk = INV_VECMAT( ijk_to_xyz ) ;
1244
1245 #if 1
1246 if( abas.verb ){
1247 DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
1248 DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
1249 }
1250 #endif
1251 }
1252
1253 /*-- make the shell of the new dataset --*/
1254
1255 if( clip_baset > 0.0f && clip_inset > 0.0f ){
1256 float fac = clip_inset / clip_baset ;
1257 abas.scale_init = fac ;
1258 if( abas.verb || fac >= 100.0 || fac <= 0.01 )
1259 INFO_message("Initial scale factor set to %.2f/%.2f=%.2g\n",
1260 clip_baset , clip_inset , fac ) ;
1261 }
1262
1263 /*===== do the hard work =====*/
1264
1265 /*-- 05 Dec 2005: initialize shift + rotation parameters? --*/
1266
1267 if( VL_coarse_rot && abas.nparam >= 6 ){
1268 float roll, pitch , yaw , sx , sy , sz ;
1269 int dxp , dyp , dzp ;
1270 MRI_IMAGE *vim ;
1271 THD_fvec3 dxyz , sxyz ;
1272
1273 if( abas.verb ) fprintf(stderr,"++ Coarse initialize ") ;
1274
1275 vim = THD_median_brick( inset ) ;
1276 (void)get_best_shiftrot( &abas , abas.imbase , vim ,
1277 &roll,&pitch,&yaw , &dxp,&dyp,&dzp ) ;
1278 mri_free(vim) ;
1279
1280 LOAD_FVEC3( dxyz , -dxp,-dyp,-dzp ) ;
1281 sxyz = MATVEC( ijk_to_xyz.mm , dxyz ) ;
1282 UNLOAD_FVEC3( sxyz , sx,sy,sz ) ;
1283
1284 init_set = 6 ;
1285 init_val[0] = sx ; init_val[1] = sy ; init_val[2] = sz ;
1286 init_val[3] = roll ; init_val[4] = pitch ; init_val[5] = yaw ;
1287
1288 if( abas.verb )
1289 fprintf(stderr,
1290 "\n++ dx=%.1f dy=%.1f dz=%.1f roll=%.1f pitch=%.1f yaw=%.1f\n",
1291 sx,sy,sz , roll,pitch,yaw ) ;
1292 }
1293
1294 if( abas.verb ) INFO_message("Beginning alignment setup\n") ;
1295
1296 /* 04 Jan 2005: set up to save the computed parameters */
1297
1298 parsave = (float **)malloc( sizeof(float *) * abas.nparam ) ;
1299 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
1300 parsave[kpar] = (float *)calloc( sizeof(float) , nvals ) ;
1301
1302 mri_warp3D_align_setup( &abas ) ;
1303
1304 /*-- open summ file, if desired --*/
1305
1306 if( W_summfile != NULL ){
1307 if( THD_is_file(W_summfile) )
1308 WARNING_message("Over-writing -summ file '%s'",W_summfile) ;
1309 if( strcmp(W_summfile,"-") == 0 ) summfp = stdout ;
1310 else summfp = fopen( W_summfile, "w" ) ;
1311 if( summfp == NULL )
1312 ERROR_message("Can't open -summ file '%s'",W_summfile) ;
1313 }
1314
1315 if( abas.verb ) INFO_message("Beginning alignment loop\n") ;
1316
1317 /**--------------- loop over input sub-bricks ---------------**/
1318
1319 for( kim=0 ; kim < nvals ; kim++ ){
1320
1321 for( kpar=0 ; kpar < abas.nparam ; kpar++ ){ /** init params **/
1322 if( abas.param[kpar].fixed )
1323 abas.param[kpar].val_init = abas.param[kpar].val_fixed ;
1324 else if( kpar < init_set )
1325 abas.param[kpar].val_init = init_val[kpar] ; /* 06 Dec 2005 */
1326 else
1327 abas.param[kpar].val_init = abas.param[kpar].ident ;
1328 }
1329
1330 /** create copy of input brick into qim
1331 then warp-align it, with result into tim **/
1332
1333 qim = mri_scale_to_float( DSET_BRICK_FACTOR(inset,kim) ,
1334 DSET_BRICK(inset,kim) ) ;
1335
1336 sdif_before = mri_scaled_diff( abas.imbase , qim , abas.imsk ) ;
1337
1338 tim = mri_warp3D_align_one( &abas , qim ) ;
1339 mri_free( qim ) ; DSET_unload_one( inset , kim ) ;
1340
1341 sdif_after = mri_scaled_diff( abas.imbase , tim , abas.imsk ) ;
1342
1343 #if 0
1344 if( abas.verb ){ DUMP_VECMAT( "end mv_for" , mv_for ) ; }
1345 #endif
1346
1347 if( abas.verb )
1348 INFO_message("#%d RMS_diff: before=%g after=%g",kim,sdif_before,sdif_after) ;
1349
1350 if( summfp != NULL )
1351 fprintf(summfp,"RMS[%d] = %g %g ITER = %d/%d\n",
1352 kim , sdif_before , sdif_after , abas.num_iter , abas.max_iter ) ;
1353
1354 /** save output parameters for later **/
1355
1356 for( kpar=0 ; kpar < abas.nparam ; kpar++ ) {
1357 parsave[kpar][kim] = abas.param[kpar].val_out ; /* 04 Jan 2005 */
1358 }
1359 /** convert output image from float to whatever **/
1360
1361 if( !null_output ){
1362 switch( DSET_BRICK_TYPE(outset,kim) ){
1363
1364 default:
1365 ERROR_message("Can't store bricks of type %s\n",
1366 MRI_TYPE_name[DSET_BRICK_TYPE(inset,kim)] ) ;
1367 /* fall thru on purpose */
1368
1369 case MRI_float:
1370 EDIT_substitute_brick( outset, kim, MRI_float, MRI_FLOAT_PTR(tim) );
1371 mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
1372 break ;
1373
1374 case MRI_short:{
1375 double fac = DSET_BRICK_FACTOR(inset,kim) ;
1376 fac = (fac > 0.0) ? 1.0/fac : 1.0 ;
1377 fim = mri_to_short(fac,tim) ; mri_free( tim ) ;
1378 EDIT_substitute_brick( outset, kim, MRI_short, MRI_SHORT_PTR(fim) );
1379 fac = (fac != 1.0) ? 1.0/fac : 0.0 ;
1380 EDIT_BRICK_FACTOR( outset , kim , fac ) ;
1381 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
1382 }
1383 break ;
1384
1385 case MRI_byte:
1386 vp = mri_min(tim) ;
1387 if( vp < 0.0f ){
1388 WARNING_message(
1389 "output sub-brick #%d is byte, but has negative values\n",
1390 kim ) ;
1391 }
1392 fim = mri_to_byte(tim) ; mri_free( tim ) ;
1393 EDIT_substitute_brick( outset, kim, MRI_byte, MRI_BYTE_PTR(fim) ) ;
1394 mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
1395 break ;
1396 }
1397 } else {
1398 mri_free(tim) ; tim = NULL ; /* 25 Jul 2007 */
1399 }
1400
1401 } /* end of loop over input sub-bricks */
1402 DSET_unload( inset ) ;
1403
1404 if( summfp != NULL && summfp != stdout ){ fclose(summfp); summfp = NULL; }
1405
1406 /*===== hard work is done =====*/
1407
1408 /* mri_warp3D_align_cleanup( &abas ) ; ZSS: This one's moved below, just in case it messes with abas which is still in use. */
1409
1410 /*-- 06 Jul 2005:
1411 write the affine transform matrices into the output header --*/
1412
1413 THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_IDCODE" ,
1414 inset->idcode.str ) ;
1415 THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_NAME" ,
1416 DSET_HEADNAME(inset) ) ;
1417 if( base_idc != NULL )
1418 THD_set_string_atr( outset->dblk , "WARPDRIVE_BASE_IDCODE" , base_idc ) ;
1419 if( wt_idc != NULL )
1420 THD_set_string_atr( outset->dblk , "WARPDRIVE_WEIGHT_IDCODE" , wt_idc ) ;
1421
1422 if( WARPDRIVE_IS_AFFINE(warpdrive_code) ){ /* can't do bilinear here */
1423 float matar[12] ; char anam[64] ;
1424
1425 for( kpar=0 ; kpar < 12 ; kpar++ ) parvec[kpar] = 0.0 ;
1426
1427 if( matrix_save_1D != NULL ){
1428 msfp = fopen(matrix_save_1D,"w") ;
1429 fprintf(msfp,
1430 "# 3dWarpDrive matrices (DICOM-to-DICOM, row-by-row):\n") ;
1431 }
1432
1433 for( kim=0 ; kim < nvals ; kim++ ){
1434 for( kpar=0 ; kpar < abas.nparam ; kpar++ ) {/* load params */
1435 parvec[kpar] = parsave[kpar][kim] ;
1436 }
1437 parset_affine() ; /* compute matrices */
1438
1439 UNLOAD_MAT(mv_for.mm,matar[0],matar[1],matar[2],
1440 matar[4],matar[5],matar[6],
1441 matar[8],matar[9],matar[10] ) ;
1442 UNLOAD_FVEC3(mv_for.vv,matar[3],matar[7],matar[11]) ;
1443 sprintf(anam,"WARPDRIVE_MATVEC_FOR_%06d",kim) ;
1444 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
1445
1446 if( msfp != NULL ) /* 25 Jul 2007 */
1447 fprintf(msfp,"%13.6g %13.6g %13.6g %13.6g "
1448 "%13.6g %13.6g %13.6g %13.6g "
1449 "%13.6g %13.6g %13.6g %13.6g\n" ,
1450 matar[0],matar[1],matar[ 2],matar[ 3],
1451 matar[4],matar[5],matar[ 6],matar[ 7],
1452 matar[8],matar[9],matar[10],matar[11] ) ;
1453
1454 UNLOAD_MAT(mv_inv.mm,matar[0],matar[1],matar[2],
1455 matar[4],matar[5],matar[6],
1456 matar[8],matar[9],matar[10] ) ;
1457 UNLOAD_FVEC3(mv_inv.vv,matar[3],matar[7],matar[11]) ;
1458 sprintf(anam,"WARPDRIVE_MATVEC_INV_%06d",kim) ;
1459 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
1460
1461 /* redo with just the shift+rotation parameters [22 Aug 2006] */
1462
1463 for( kpar=0 ; kpar < 3 ; kpar++ ) parvec[kpar] = 0.0 ;
1464 for( kpar=6 ; kpar < abas.nparam ; kpar++ ) parvec[kpar] = 0.0 ;
1465 parset_affine() ;
1466
1467 UNLOAD_MAT(mv_for.mm,matar[0],matar[1],matar[2],
1468 matar[4],matar[5],matar[6],
1469 matar[8],matar[9],matar[10] ) ;
1470 UNLOAD_FVEC3(mv_for.vv,matar[3],matar[7],matar[11]) ;
1471 sprintf(anam,"WARPDRIVE_ROTMAT_FOR_%06d",kim) ;
1472 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
1473
1474 UNLOAD_MAT(mv_inv.mm,matar[0],matar[1],matar[2],
1475 matar[4],matar[5],matar[6],
1476 matar[8],matar[9],matar[10] ) ;
1477 UNLOAD_FVEC3(mv_inv.vv,matar[3],matar[7],matar[11]) ;
1478 sprintf(anam,"WARPDRIVE_ROTMAT_INV_%06d",kim) ;
1479 THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
1480 }
1481 if( msfp != NULL ) fclose(msfp) ;
1482 }
1483
1484 /*-- write the results to disk for all of history to see --*/
1485
1486 if( !null_output ){
1487 DSET_write( outset ) ; DSET_unload( outset ) ;
1488 if( abas.verb ) WROTE_DSET(outset) ;
1489 } else {
1490 INFO_message("-prefix NULL was given ==> no output dataset written") ;
1491 }
1492
1493 if( W_1Dfile != NULL ){
1494 FILE *fp ;
1495 if( abas.verb ) INFO_message("Writing 1Dfile: %s\n",W_1Dfile) ;
1496 if( THD_is_file(W_1Dfile) )
1497 WARNING_message("Overwriting file %s\n",W_1Dfile) ;
1498
1499 fp = fopen( W_1Dfile , "w" ) ;
1500 if( fp != NULL ){
1501
1502 fprintf(fp,"#") ;
1503 for( kim=0 ; kim < argc ; kim++ ) fprintf(fp," %s",argv[kim]) ;
1504 fprintf(fp,"\n") ;
1505
1506 fprintf(fp,"#") ;
1507 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
1508 fprintf(fp," %-13.13s",abas.param[kpar].name) ;
1509 fprintf(fp,"\n") ;
1510
1511 fprintf(fp,"#") ;
1512 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
1513 fprintf(fp," -------------") ;
1514 fprintf(fp,"\n") ;
1515
1516 for( kim=0 ; kim < nvals ; kim++ ){
1517 for( kpar=0 ; kpar < abas.nparam ; kpar++ )
1518 fprintf(fp," %13.6g",parsave[kpar][kim]) ;
1519 fprintf(fp,"\n") ;
1520 }
1521 fclose(fp) ;
1522 }
1523
1524 }
1525
1526 if( abas.verb ){
1527 double tt = (NI_clock_time()-ctstart)*0.001 ;
1528 INFO_message("Total elapsed time = %.2f s\n",tt) ;
1529 }
1530
1531 mri_warp3D_align_cleanup( &abas ) ;
1532
1533 exit(0) ;
1534 }
1535
1536 /*==========================================================================*/
1537 /*======= Code to find best initial shift_rotate parameters, crudely =======*/
1538 /*==========================================================================*/
1539
1540 /*--------------------------------------------------------------------
1541 Calculate
1542 ( [ 2 ] )
1543 min ( sum [ {a v(i-dx,j-dy,k-dz) - b(i,j,k)} w(i,j,k) ] )
1544 a ( ijk [ ] )
1545
1546 where the sum is taken over voxels at least 'edge' in from the edge.
1547 'edge' must be bigger than max(|dx|,|dy|,|dz|). The weight w may
1548 be NULL, in which case it is taken to be identically 1.
1549 ----------------------------------------------------------------------*/
1550
1551 #define B(i,j,k) b[(i)+(j)*nx+(k)*nxy]
1552 #define V(i,j,k) v[(i)+(j)*nx+(k)*nxy]
1553 #define W(i,j,k) w[(i)+(j)*nx+(k)*nxy]
1554
voldif(int nx,int ny,int nz,float * b,int dx,int dy,int dz,float * v,int edge,float * w)1555 static float voldif( int nx, int ny, int nz, float *b,
1556 int dx, int dy, int dz, float *v, int edge , float *w )
1557 {
1558 int nxy=nx*ny, nxtop=nx-edge, nytop=ny-edge, nztop=nz-edge , ii,jj,kk ;
1559 float bbsum=0.0f , vvsum=0.0f , bvsum=0.0f , bb,vv,ww ;
1560
1561 if( w == NULL ){ /** no weight given **/
1562 for( kk=edge ; kk < nztop ; kk++ ){
1563 for( jj=edge ; jj < nytop ; jj++ ){
1564 for( ii=edge ; ii < nxtop ; ii++ ){
1565 bb = B(ii,jj,kk) ; vv = V(ii-dx,jj-dy,kk-dz) ;
1566 bbsum += bb*bb ; vvsum += vv*vv ; bvsum += bb*vv ;
1567 }}}
1568 } else { /** use given weight **/
1569 for( kk=edge ; kk < nztop ; kk++ ){
1570 for( jj=edge ; jj < nytop ; jj++ ){
1571 for( ii=edge ; ii < nxtop ; ii++ ){
1572 ww = W(ii,jj,kk) ;
1573 if( ww > 0.0f ){
1574 bb = B(ii,jj,kk) ; vv = V(ii-dx,jj-dy,kk-dz) ;
1575 bbsum += ww*bb*bb ; vvsum += ww*vv*vv ; bvsum += ww*bb*vv ;
1576 }
1577 }}}
1578 }
1579
1580 if( vvsum > 0.0f ) bbsum -= bvsum*bvsum/vvsum ;
1581 return bbsum ;
1582 }
1583
1584 /*---------------------------------------------------------------------
1585 Do some shifts to find the best starting point for registration
1586 (globals VL_coarse_del and VL_coarse_num control operations).
1587 -----------------------------------------------------------------------*/
1588
get_best_shift(int nx,int ny,int nz,float * b,float * v,float * w,int * dxp,int * dyp,int * dzp)1589 static float get_best_shift( int nx, int ny, int nz,
1590 float *b, float *v , float *w ,
1591 int *dxp , int *dyp , int *dzp )
1592 {
1593 int bdx=0 , bdy=0 , bdz=0 , dx,dy,dz , nxyz=nx*ny*nz ;
1594 float bsum , sum ;
1595
1596 int shift = VL_coarse_del, numsh = VL_coarse_num,
1597 shtop = shift*numsh , edge = shtop+shift , sqtop = shtop*shtop ;
1598
1599 bsum = 0.0 ;
1600 for( dx=0 ; dx < nxyz ; dx++ ) bsum += b[dx]*b[dx] ;
1601
1602 for( dz=-shtop ; dz <= shtop ; dz+=shift ){
1603 for( dy=-shtop ; dy <= shtop ; dy+=shift ){
1604 for( dx=-shtop ; dx <= shtop ; dx+=shift ){
1605 if( dx*dx+dy*dy+dz*dz > sqtop ) continue ;
1606 sum = voldif( nx,ny,nz , b , dx,dy,dz , v , edge , w ) ;
1607 if( sum < bsum ){ bsum = sum; bdx = dx; bdy = dy; bdz = dz; }
1608 }}}
1609
1610 *dxp = bdx ; *dyp = bdy ; *dzp = bdz ; return bsum ;
1611 }
1612
1613 /*----------------------------------------------------------------------*/
1614 /* Find best angles AND best shifts all at once't. Will only be
1615 called if VL_coarse_rot is nonzero.
1616 ------------------------------------------------------------------------*/
1617
1618 #define DANGLE 9.0f
1619 #define NROLL 1
1620 #define NPITCH 2
1621 #define NYAW 1
1622
get_best_shiftrot(MRI_warp3D_align_basis * bas,MRI_IMAGE * base,MRI_IMAGE * vol,float * roll,float * pitch,float * yaw,int * dxp,int * dyp,int * dzp)1623 static float get_best_shiftrot( MRI_warp3D_align_basis *bas ,
1624 MRI_IMAGE *base , MRI_IMAGE *vol ,
1625 float *roll , float *pitch , float *yaw ,
1626 int *dxp , int *dyp , int *dzp )
1627 {
1628 int ii,jj,kk ;
1629 float r,p,y , br=0.0f , bp=0.0f , by=0.0f ;
1630 float bsum=1.e+38 , sum ;
1631 MRI_IMAGE *tim , *wim=NULL , *bim, *vim , *msk ;
1632 float *bar , *tar , *var , dif , *www=NULL , wtop , param[12] ;
1633 byte *mmm ;
1634 int nx,ny,nz , sx,sy,sz , bsx=0,bsy=0,bsz=0 , nxy,nxyz , subd=0 ;
1635 THD_vecmat ijk_to_xyz_save , xyz_to_ijk_save ;
1636 THD_mat33 mat ;
1637
1638 *roll = *pitch = *yaw = 0.0f ; /* in case of sudden death */
1639 *dxp = *dyp = *dzp = 0 ;
1640 memset(&xyz_to_ijk_save,0,sizeof(THD_vecmat));
1641 memset(&ijk_to_xyz_save,0,sizeof(THD_vecmat));
1642 if( bas->nparam < 6 ) return 0.0f ; /* no rotations allowed? */
1643
1644 for( ii=0 ; ii < 12 ; ii++ ) param[ii] = 0.0f ;
1645 param[6] = param[7] = param[8] = 1.0f ;
1646
1647 nx = base->nx ; ny = base->ny ; nz = base->nz ; nxy = nx*ny ;
1648
1649 /** if image volume is big enough, sub-sample by 2 for speedup **/
1650
1651 bim = base ; vim = vol ;
1652
1653 #if 1
1654 if( nx >= 120 && ny >= 120 && nz >= 120 ){
1655 int hnx=(nx+1)/2 , hny=(ny+1)/2 , hnz=(nz+1)/2 , hnxy=hnx*hny ;
1656
1657 if( bas->verb ) fprintf(stderr,"x") ;
1658
1659 /* copy and blur base, then subsample it into new image bim */
1660
1661 tim = mri_copy(base) ; tar = MRI_FLOAT_PTR(tim) ;
1662 FIR_blur_volume( nx,ny,nz , 1.0f,1.0f,1.0f , tar , 1.0f ) ;
1663 bim = mri_new_vol( hnx,hny,hnz , MRI_float ) ; bar = MRI_FLOAT_PTR(bim) ;
1664 for( kk=0 ; kk < hnz ; kk++ ) /* subsampling */
1665 for( jj=0 ; jj < hny ; jj++ )
1666 for( ii=0 ; ii < hnx ; ii++ )
1667 bar[ii+jj*hnx+kk*hnxy] = tar[2*(ii+jj*nx+kk*nxy)] ;
1668 mri_free(tim) ;
1669
1670 /* copy and blur vol, then subsample it into a new image vim */
1671
1672 tim = mri_copy(vol) ; tar = MRI_FLOAT_PTR(tim) ;
1673 FIR_blur_volume( nx,ny,nz , 1.0f,1.0f,1.0f , tar , 1.0f ) ;
1674 vim = mri_new_vol( hnx,hny,hnz , MRI_float ) ; var = MRI_FLOAT_PTR(vim) ;
1675 for( kk=0 ; kk < hnz ; kk++ ) /* subsampling */
1676 for( jj=0 ; jj < hny ; jj++ )
1677 for( ii=0 ; ii < hnx ; ii++ )
1678 var[ii+jj*hnx+kk*hnxy] = tar[2*(ii+jj*nx+kk*nxy)] ;
1679 mri_free(tim) ;
1680
1681 /* adjust grid spacing in new images */
1682
1683 bim->dx = vim->dx = 2.0f * base->dx ;
1684 bim->dy = vim->dy = 2.0f * base->dy ;
1685 bim->dz = vim->dz = 2.0f * base->dz ;
1686
1687 nx = hnx; ny = hny; nz = hnz; nxy = hnxy; subd = 2;
1688 VL_coarse_del /= 2 ;
1689
1690 /* scale up the index-to-coord transformation */
1691
1692 ijk_to_xyz_save = ijk_to_xyz ;
1693 xyz_to_ijk_save = xyz_to_ijk ;
1694 mat = ijk_to_xyz.mm ;
1695 ijk_to_xyz.mm = MAT_SCALAR( mat , 2.0f ) ;
1696 xyz_to_ijk = INV_VECMAT( ijk_to_xyz ) ;
1697 }
1698 #endif
1699
1700 /* make a weighting image (blurred & masked copy of base) */
1701
1702 if( bas->verb ) fprintf(stderr,"w") ;
1703
1704 wim = mri_copy(bim) ; www = MRI_FLOAT_PTR(wim) ; nxyz = nx*ny*nz ;
1705 for( ii=0 ; ii < nxyz ; ii++ ) www[ii] = fabsf(www[ii]) ;
1706 FIR_blur_volume( nx,ny,nz , 1.0f,1.0f,1.0f , www , 1.0f ) ;
1707 wtop = 0.0f ;
1708 for( ii=0 ; ii < nxyz ; ii++ ) wtop = MAX(wtop,www[ii]) ;
1709 wtop = 1.0f / wtop ;
1710 for( ii=0 ; ii < nxyz ; ii++ ){
1711 www[ii] *= wtop ; if( www[ii] < 0.05 ) www[ii] = 0.0f ;
1712 }
1713 mmm = mri_automask_image( wim ) ;
1714 for( ii=0 ; ii < nxyz ; ii++ ) if( mmm[ii] == 0 ) www[ii] = 0.0f ;
1715 msk = mri_empty_conforming( bim , MRI_byte ) ;
1716 mri_fix_data_pointer( mmm , msk ) ;
1717
1718 if( bas->verb )
1719 fprintf(stderr,"[%.1f%%]" , (100.0*THD_countmask(nxyz,mmm))/nxyz ) ;
1720
1721 /* prepare to rotate and shift the night away */
1722
1723 mri_warp3D_set_womask( msk ) ;
1724 mri_warp3D_method( MRI_NN ) ;
1725 bar = MRI_FLOAT_PTR(bim) ;
1726
1727 for( kk=-NROLL ; kk <= NROLL ; kk++ ){
1728 for( jj=-NPITCH ; jj <= NPITCH ; jj++ ){
1729 for( ii=-NYAW ; ii <= NYAW ; ii++ ){
1730 r = kk*DANGLE ; p = jj*DANGLE ; y = ii*DANGLE ;
1731
1732 if( r == 0.0f && p == 0.0f && y == 0.0f ){ /* no rotate */
1733 tim = vim ;
1734 } else { /* rotate vim */
1735 param[3] = r ; param[4] = p ; param[5] = y ;
1736 bas->vwset( 12 , param ) ;
1737 tim = mri_warp3D( vim , 0,0,0 , bas->vwfor ) ;
1738 }
1739 tar = MRI_FLOAT_PTR(tim) ;
1740 sum = get_best_shift( nx,ny,nz, bar, tar, www, &sx,&sy,&sz ) ;
1741 if( bas->verb ) fprintf(stderr,"%s",(sum<bsum)?"*":".") ;
1742 if( sum < bsum ){
1743 br=r ; bp=p ; by=y ; bsx=sx ; bsy=sy; bsz=sz ; bsum=sum ;
1744 }
1745 if( tim != vim ) mri_free(tim) ;
1746 }}}
1747
1748 /* cleanup and exeunt */
1749
1750 mri_free(wim) ; mri_free(msk) ; mri_warp3D_set_womask(NULL) ;
1751 if( subd ){
1752 mri_free(bim); mri_free(vim);
1753 bsx *= 2; bsy *= 2; bsz *= 2; VL_coarse_del *=2;
1754 ijk_to_xyz = ijk_to_xyz_save ;
1755 xyz_to_ijk = xyz_to_ijk_save ;
1756 }
1757
1758 *roll = br ; *pitch = bp ; *yaw = by ;
1759 *dxp = bsx; *dyp = bsy; *dzp = bsz ;
1760
1761 return bsum ;
1762 }
1763