1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 #undef ASSIF
10 #define ASSIF(b,a) if( (b) != NULL ) *(b) = (a)
11
12 /*** NOT 7D SAFE ***/
13
14 /*************************************************************************
15 ** 3D version of mri_2dalign.c -- RWCox -- October 1998 **
16 *************************************************************************/
17
18 #define MAX_ITER 5
19 #define DXY_THRESH 0.07 /* pixels */
20 #define PHI_THRESH 0.21 /* degrees */
21 #define DFAC (PI/180.0)
22
23 static float dxy_thresh = DXY_THRESH ,
24 phi_thresh = PHI_THRESH ,
25 delfac = 1.5 ;
26
27 static int max_iter = MAX_ITER ;
28 static int ax1 = 0 , ax2 = 1 , ax3 = 2 ;
29 static int dcode = -1 ;
30
31 static int wproc = 0 ; /* 06 Jun 2002: process imwt? */
32 static int wtrim = 0 ; /* 06 Jun 2002: trimming stuff */
33
34 static float sinit = 1.0 ; /* 22 Mar 2004: init scale? */
35
36 #define DOTRIM (basis->xa >= 0)
37
38 #define IMTRIM(qqq) mri_cut_3D( qqq , basis->xa,basis->xb , \
39 basis->ya,basis->yb , \
40 basis->za,basis->zb )
41
42 /*! Macro to replace an image with a trimmed copy, if needed. */
43
44 #define TRIM(imq) \
45 do{ if( DOTRIM ){ \
46 MRI_IMAGE *qim = IMTRIM(imq) ; \
47 mri_free(imq) ; imq = qim ; \
48 } } while(0)
49
50 /*! Macro to print out a volume range in verbose mode. */
51
52 #define VRANG(str,imq) \
53 do{ if(verbose){ \
54 double tt=mri_max(imq), bb=mri_min(imq) ; \
55 fprintf(stderr," %s range: min=%g max=%g\n",str,bb,tt); } } while(0)
56
57 /*--------------------------------------------------------------------*/
58
mri_3dalign_wtrimming(int ttt)59 void mri_3dalign_wtrimming( int ttt ){ wtrim = ttt; } /* 06 Jun 2002 */
60
mri_3dalign_wproccing(int ttt)61 void mri_3dalign_wproccing( int ttt ){ wproc = ttt; } /* 06 Jun 2002 */
62
mri_3dalign_scaleinit(float ttt)63 void mri_3dalign_scaleinit( float ttt ) /* 22 Mar 2004 */
64 {
65 if( ttt > 0.0 ) sinit = ttt ;
66 }
67
68 /*--------------------------------------------------------------------*/
69
mri_3dalign_params(int maxite,float dxy,float dph,float dfac,int bx1,int bx2,int bx3,int dc)70 void mri_3dalign_params( int maxite ,
71 float dxy , float dph , float dfac ,
72 int bx1 , int bx2 , int bx3 , int dc )
73 {
74 if( maxite > 0 ) max_iter = maxite ; else max_iter = MAX_ITER ;
75 if( dxy > 0.0 ) dxy_thresh = dxy ; else dxy_thresh = DXY_THRESH ;
76 if( dph > 0.0 ) phi_thresh = dph ; else phi_thresh = PHI_THRESH ;
77 if( dfac > 0.0 ) delfac = dfac ;
78
79 if( bx1 >= 0 && bx1 <= 2 ) ax1 = bx1 ;
80 if( bx2 >= 0 && bx2 <= 2 ) ax2 = bx2 ;
81 if( bx3 >= 0 && bx3 <= 2 ) ax3 = bx3 ;
82
83 dcode = dc ;
84 return ;
85 }
86
87 /*--------------------------------------------------------------------*/
88
89 static float init_dth1=0.0 , init_dth2=0.0 , init_dth3=0.0 ;
90 static float init_dx =0.0 , init_dy =0.0 , init_dz =0.0 ;
91
92 #define CLEAR_INITVALS mri_3dalign_initvals(0.0,0.0,0.0,0.0,0.0,0.0)
93
94 #define NONZERO_INITVALS \
95 ( init_dth1 != 0.0 || init_dth2 != 0.0 || init_dth3 != 0.0 || \
96 init_dx != 0.0 || init_dy != 0.0 || init_dz != 0.0 )
97
mri_3dalign_initvals(float th1,float th2,float th3,float dx,float dy,float dz)98 void mri_3dalign_initvals( float th1,float th2,float th3 ,
99 float dx ,float dy ,float dz )
100 {
101 init_dth1 = th1 ; init_dth2 = th2 ; init_dth3 = th3 ; /* degrees */
102 init_dx = dx ; init_dy = dy ; init_dz = dz ; /* mm */
103 }
104
105 /*--------------------------------------------------------------------*/
106
107 static int regmode = MRI_QUINTIC ;
108 static int verbose = 0 ;
109 static int noreg = 0 ;
110 static int clipit = 0 ;
111
mri_3dalign_method(int rmode,int verb,int norgg,int clip)112 void mri_3dalign_method( int rmode , int verb , int norgg , int clip )
113 {
114 regmode = rmode ;
115 verbose = verb ;
116 noreg = norgg ;
117 clipit = clip ;
118 return ;
119 }
120
121 /*-------------------------------------------------------------------*/
122
123 static float blurit = 0.0 ;
mri_3dalign_blurring(float bl)124 void mri_3dalign_blurring( float bl ){ blurit = bl ; return ; }
125
126 static int final_regmode = -1 ; /* 20 Nov 1998 */
mri_3dalign_final_regmode(int frm)127 void mri_3dalign_final_regmode( int frm )
128 {
129 final_regmode = frm ;
130 return ;
131 }
132
133 /*-------------------------------------------------------------------*/
134
135 static int xedge=-1 , yedge=-1 , zedge=-1 ;
136 static int xfade , yfade , zfade ;
137
138 static int force_edging=0 ;
139
mri_3dalign_edging(int x,int y,int z)140 void mri_3dalign_edging( int x , int y , int z ) /* 10 Dec 2000 */
141 {
142 xedge = x ; yedge = y ; zedge = z ;
143 }
144
mri_3dalign_force_edging(int n)145 void mri_3dalign_force_edging( int n )
146 {
147 force_edging = n ;
148 }
149
mri_3dalign_edging_default(int nx,int ny,int nz)150 void mri_3dalign_edging_default( int nx , int ny , int nz )
151 {
152 char *ef=my_getenv("AFNI_VOLREG_EDGING") , *eq ;
153
154 if( ef == NULL ){ /* the 5% solution */
155 xfade = (int)(0.05*nx+0.5) ;
156 yfade = (int)(0.05*ny+0.5) ;
157 zfade = (int)(0.05*nz+0.5) ;
158 } else {
159 float ff = strtod(ef,&eq) ;
160 if( ff < 0 ){ /* again */
161 xfade = (int)(0.05*nx+0.5) ;
162 yfade = (int)(0.05*ny+0.5) ;
163 zfade = (int)(0.05*nz+0.5) ;
164 } else {
165 if( *eq == '%' ){ /* the whatever % solution */
166 xfade = (int)(0.01*ff*nx+0.5) ;
167 yfade = (int)(0.01*ff*ny+0.5) ;
168 zfade = (int)(0.01*ff*nz+0.5) ;
169 } else { /* the fixed value solution */
170 xfade = (int)( MIN(0.25*nx,ff) ) ;
171 yfade = (int)( MIN(0.25*ny,ff) ) ;
172 zfade = (int)( MIN(0.25*nz,ff) ) ;
173 }
174 }
175 }
176 }
177
178 /*--------------------------------------------------------------------
179 Inputs: imbase = base image for alignment
180 imwt = image of weight factors to align to
181 (if NULL, will generate one internally)
182
183 Output: pointer to a MRI_3dalign_basis struct, for later use.
184 The malloc-ed data in there can be freed using
185 routine MRI_3dalign_cleanup.
186 ----------------------------------------------------------------------*/
187
mri_3dalign_setup(MRI_IMAGE * imbase,MRI_IMAGE * imwt)188 MRI_3dalign_basis * mri_3dalign_setup( MRI_IMAGE *imbase , MRI_IMAGE *imwt )
189 {
190 MRI_IMAGE *bim , *pim , *mim , *dim , *imww , *cim ;
191 float *dar , *par , *mar ;
192 float delta , dx,dy,dz ;
193 int ii ;
194 MRI_IMARR * fitim =NULL;
195 double * chol_fitim=NULL ;
196 MRI_3dalign_basis *basis = NULL ;
197
198 ENTRY("mri_3dalign_setup") ;
199
200 if( !MRI_IS_3D(imbase) ){
201 fprintf(stderr,"\n*** mri_3dalign_setup: cannot use nD images!\a\n") ;
202 RETURN( NULL );
203 }
204
205 /*--- create output struct ---*/
206
207 /* from malloc 12 Feb 2009 [lesstif patrol] */
208 basis = (MRI_3dalign_basis *) calloc( 1, sizeof(MRI_3dalign_basis) ) ;
209
210 /*-- local copy of input image --*/
211
212 cim = mri_to_float( imbase ) ;
213
214 dx = fabs(cim->dx) ; if( dx == 0.0 ) dx = 1.0 ;
215 dy = fabs(cim->dy) ; if( dy == 0.0 ) dy = 1.0 ;
216 dz = fabs(cim->dz) ; if( dz == 0.0 ) dz = 1.0 ;
217
218 /*--- get the weighting image ---*/
219
220 if( imwt != NULL &&
221 (imwt->nx != cim->nx || imwt->ny != cim->ny || imwt->nz != cim->nz) ){
222
223 fprintf(stderr,"*** WARNING: in mri_3dalign_setup, weight image mismatch!\n") ;
224 imwt = NULL ;
225 }
226
227 /* make weight up from the base if it isn't supplied */
228
229 if( imwt == NULL ){
230 int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny , nxyz=nxy*nz ;
231 int ii ;
232 float * f , clip ;
233
234 /* copy base image */
235
236 imww = mri_to_float( cim ) ; f = MRI_FLOAT_PTR(imww) ;
237
238 if( verbose ) fprintf(stderr," initializing weight") ;
239
240 for( ii=0 ; ii < nxyz ; ii++ ) f[ii] = fabs(f[ii]) ; /* 16 Nov 1998 */
241
242 #if 1
243 EDIT_blur_volume_3d( nx,ny,nz , dx,dy,dz ,
244 MRI_float , f , 3.0*dx , 3.0*dy , 3.0*dz ) ;
245 #else
246 MRI_5blur_inplace_3D( imww ) ; /* 07 Jun 2002 */
247 #endif
248
249 if( verbose ) fprintf(stderr,":") ;
250
251 clip = 0.025 * mri_max(imww) ;
252 for( ii=0 ; ii < nxyz ; ii++ ) if( f[ii] < clip ) f[ii] = 0.0 ;
253
254 } else {
255 imww = mri_to_float( imwt ) ; /* just copy input weight image */
256
257 if( wproc ){ /* 06 Jun 2002: process input weight */
258 int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny , nxyz=nxy*nz ;
259 int ii ;
260 float * f , clip ;
261
262 if( verbose ) fprintf(stderr," processing weight") ;
263 f = MRI_FLOAT_PTR(imww) ;
264 for( ii=0 ; ii < nxyz ; ii++ ) f[ii] = fabs(f[ii]) ; /* 16 Nov 1998 */
265 #if 1
266 EDIT_blur_volume_3d( nx,ny,nz , dx,dy,dz ,
267 MRI_float , f , 3.0*dx , 3.0*dy , 3.0*dz ) ;
268 #else
269 MRI_5blur_inplace_3D( imww ) ; /* 07 Jun 2002 */
270 #endif
271 if( verbose ) fprintf(stderr,":") ;
272 clip = 0.025 * mri_max(imww) ;
273 for( ii=0 ; ii < nxyz ; ii++ ) if( f[ii] < clip ) f[ii] = 0.0 ;
274 }
275 }
276
277 /*-- 10 Dec 2000: user-controlled fade out around the edges --*/
278
279 if( imwt == NULL || force_edging ){
280 int ff , ii,jj,kk ;
281 int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny ;
282 float *f = MRI_FLOAT_PTR(imww) ;
283
284 xfade = xedge ; yfade = yedge ; zfade = zedge ; /* static variables */
285
286 if( xfade < 0 || yfade < 0 || zfade < 0 )
287 mri_3dalign_edging_default(nx,ny,nz) ; /* reassign fades */
288
289 #define FF(i,j,k) f[(i)+(j)*nx+(k)*nxy]
290
291 for( jj=0 ; jj < ny ; jj++ )
292 for( ii=0 ; ii < nx ; ii++ )
293 for( ff=0 ; ff < zfade ; ff++ )
294 FF(ii,jj,ff) = FF(ii,jj,nz-1-ff) = 0.0 ;
295
296 for( kk=0 ; kk < nz ; kk++ )
297 for( jj=0 ; jj < ny ; jj++ )
298 for( ff=0 ; ff < xfade ; ff++ )
299 FF(ff,jj,kk) = FF(nx-1-ff,jj,kk) = 0.0 ;
300
301 for( kk=0 ; kk < nz ; kk++ )
302 for( ii=0 ; ii < nx ; ii++ )
303 for( ff=0 ; ff < yfade ; ff++ )
304 FF(ii,ff,kk) = FF(ii,ny-1-ff,kk) = 0.0 ;
305 }
306
307 /*-- 06 Jun 2002: compute wtrimmed volume size --*/
308
309 basis->xa = -1 ; /* flag for no wtrim */
310
311 if( wtrim ){
312 int xa=-1,xb , ya,yb , za,zb ;
313 MRI_autobbox( imww , &xa,&xb , &ya,&yb , &za,&zb ) ;
314 if( xa >= 0 ){
315 float nxyz = imww->nx * imww->ny * imww->nz ;
316 float nttt = (xb-xa+1)*(yb-ya+1)*(zb-za+1) ;
317 float trat = 100.0 * nttt / nxyz ;
318
319 if( verbose )
320 fprintf(stderr," wtrim: [%d..%d]x[%d..%d]x[%d..%d]"
321 " = %d voxels kept, out of %d (%.1f%%)\n" ,
322 xa,xb , ya,yb , za,zb , (int)nttt , (int)nxyz , trat ) ;
323
324 /* keep trimming if saves at least 10% per volume */
325
326 if( trat < 90.0 ){
327 basis->xa = xa ; basis->xb = xb ;
328 basis->ya = ya ; basis->yb = yb ;
329 basis->za = za ; basis->zb = zb ;
330
331 TRIM(imww) ;
332 } else if( verbose ){
333 fprintf(stderr," skipping use of trim - too little savings\n");
334 }
335 }
336 }
337
338 VRANG("weight",imww) ;
339 if( verbose ){
340 float *f=MRI_FLOAT_PTR(imww) , perc ;
341 int ii , nxyz = imww->nvox , nzer=0 ;
342 for( ii=0 ; ii < nxyz ; ii++ ) nzer += (f[ii] == 0.0) ;
343 perc = (100.0*nzer)/nxyz ;
344 fprintf(stderr," # zero weights=%d out of %d (%.1f%%)\n",nzer,nxyz,perc);
345 }
346
347 /*-- base image --*/
348
349 INIT_IMARR( fitim ) ; /* array of fitting images */
350
351 if( DOTRIM ) bim = IMTRIM(cim) ; /* a trimmed duplicate */
352 else bim = cim ; /* base image */
353 ADDTO_IMARR( fitim , bim ) ;
354
355 THD_rota_method( regmode ) ;
356
357 #ifndef MEGA
358 #define MEGA (1024*1024)
359 #endif
360 if( verbose ) fprintf(stderr ,
361 " mri_3dalign: using %d Mbytes of workspace\n" ,
362 (int)(10 * bim->nvox * bim->pixel_size / MEGA) ) ;
363
364 /*-- d/d(th1) image [angles in degrees here] --*/
365
366 if( verbose ) fprintf(stderr," initializing d/d(th1)") ;
367
368 delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
369 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
370
371 pim = THD_rota3D( cim , ax1,delta , ax2,0.0 , ax3,0.0 ,
372 dcode , 0.0 , 0.0 , 0.0 ) ;
373 if( verbose ) fprintf(stderr,":") ;
374
375 mim = THD_rota3D( cim , ax1,-delta , ax2,0.0 , ax3,0.0 ,
376 dcode , 0.0 , 0.0 , 0.0 ) ;
377 if( verbose ) fprintf(stderr,":") ;
378
379 dim = mri_new_conforming( cim , MRI_float ) ;
380 delta = sinit * 0.5 * DFAC / delta ;
381 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
382 for( ii=0 ; ii < dim->nvox ; ii++ )
383 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
384 mri_free(pim) ; mri_free(mim) ;
385 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
386 VRANG("d/d(th1)",dim) ;
387
388 /*-- d/d(th2) image --*/
389
390 if( verbose ) fprintf(stderr," initializing d/d(th2)") ;
391
392 delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
393 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
394
395 pim = THD_rota3D( cim , ax1,0.0 , ax2,delta , ax3,0.0 ,
396 dcode , 0.0 , 0.0 , 0.0 ) ;
397 if( verbose ) fprintf(stderr,":") ;
398
399 mim = THD_rota3D( cim , ax1,0.0 , ax2,-delta , ax3,0.0 ,
400 dcode , 0.0 , 0.0 , 0.0 ) ;
401 if( verbose ) fprintf(stderr,":") ;
402
403 dim = mri_new_conforming( cim , MRI_float ) ;
404 delta = sinit * 0.5 * DFAC / delta ;
405 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
406 for( ii=0 ; ii < dim->nvox ; ii++ )
407 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
408 mri_free(pim) ; mri_free(mim) ;
409 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
410 VRANG("d/d(th2)",dim) ;
411
412 /*-- d/d(th3) image --*/
413
414 if( verbose ) fprintf(stderr," initializing d/d(th3)") ;
415
416 delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
417 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
418
419 pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,delta ,
420 dcode , 0.0 , 0.0 , 0.0 ) ;
421 if( verbose ) fprintf(stderr,":") ;
422
423 mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,-delta ,
424 dcode , 0.0 , 0.0 , 0.0 ) ;
425 if( verbose ) fprintf(stderr,":") ;
426
427 dim = mri_new_conforming( cim , MRI_float ) ;
428 delta = sinit * 0.5 * DFAC / delta ;
429 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
430 for( ii=0 ; ii < dim->nvox ; ii++ )
431 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
432 mri_free(pim) ; mri_free(mim) ;
433 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
434 VRANG("d/d(th3)",dim) ;
435
436 /*-- d/dx image --*/
437
438 if( verbose ) fprintf(stderr," initializing d/dx") ;
439
440 delta = delfac * dx ;
441 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
442
443 pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
444 dcode , delta , 0.0 , 0.0 ) ;
445 if( verbose ) fprintf(stderr,":") ;
446
447 mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
448 dcode , -delta , 0.0 , 0.0 ) ;
449 if( verbose ) fprintf(stderr,":") ;
450
451 dim = mri_new_conforming( cim , MRI_float ) ;
452 delta = sinit * 0.5 / delta ;
453 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
454 for( ii=0 ; ii < dim->nvox ; ii++ )
455 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
456 mri_free(pim) ; mri_free(mim) ;
457 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
458 VRANG("d/dx",dim) ;
459
460 /*-- d/dy image --*/
461
462 if( verbose ) fprintf(stderr," initializing d/dy") ;
463
464 delta = delfac * dy ;
465 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
466
467 pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
468 dcode , 0.0 , delta , 0.0 ) ;
469 if( verbose ) fprintf(stderr,":") ;
470
471 mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
472 dcode , 0.0 , -delta , 0.0 ) ;
473 if( verbose ) fprintf(stderr,":") ;
474
475 dim = mri_new_conforming( cim , MRI_float ) ;
476 delta = sinit * 0.5 / delta ;
477 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
478 for( ii=0 ; ii < dim->nvox ; ii++ )
479 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
480 mri_free(pim) ; mri_free(mim) ;
481 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
482 VRANG("d/dy",dim) ;
483
484 /*-- d/dz image --*/
485
486 if( verbose ) fprintf(stderr," initializing d/dz") ;
487
488 delta = delfac * dz ;
489 if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
490
491 pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
492 dcode , 0.0 , 0.0 , delta ) ;
493 if( verbose ) fprintf(stderr,":") ;
494
495 mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
496 dcode , 0.0 , 0.0 , -delta ) ;
497 if( verbose ) fprintf(stderr,":") ;
498
499 dim = mri_new_conforming( cim , MRI_float ) ;
500 delta = sinit * 0.5 / delta ;
501 dar = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
502 for( ii=0 ; ii < dim->nvox ; ii++ )
503 dar[ii] = delta * ( mar[ii] - par[ii] ) ;
504 mri_free(pim) ; mri_free(mim) ;
505 TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
506 VRANG("d/dz",dim) ;
507
508 /*-- done with input copy, unless it is same as base for lsqfit --*/
509
510 if( cim != bim ) mri_free(cim) ;
511
512 /*-- initialize linear least squares --*/
513
514 if( verbose ) fprintf(stderr," initializing least squares\n") ;
515
516 chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
517
518 mri_free(imww) ;
519
520 /*-- save stuff --*/
521
522 basis->fitim = fitim ;
523 basis->chol_fitim = chol_fitim ;
524
525 RETURN( basis );
526 }
527
528 /*-----------------------------------------------------------------------
529 Input: basis = MRI_3dalign_basis * return from setup routine above.
530 im = MRI_IMAGE * to align to base image
531
532 Output: Return value is aligned image;
533 *dx, *dy, *dz, *th1, *th2, *th3 are set to estimated
534 alignment parameters. Note that returned image is floats.
535
536 20 May 2009: modified to return a complex-valued image if the input
537 is complex-valued (estimation is done on abs image)
538 -------------------------------------------------------------------------*/
539
mri_3dalign_one(MRI_3dalign_basis * basis,MRI_IMAGE * im,float * th1,float * th2,float * th3,float * dx,float * dy,float * dz)540 MRI_IMAGE * mri_3dalign_one( MRI_3dalign_basis *basis , MRI_IMAGE *im ,
541 float *th1 , float *th2 , float *th3 ,
542 float *dx , float *dy , float *dz )
543 {
544 MRI_IMARR *fitim ;
545 double *chol_fitim=NULL ;
546 float *fit , *dfit ;
547 int iter , good , ii ;
548 float dxt , dyt , dzt , ftop,fbot ;
549 MRI_IMAGE *tim=NULL , *fim ;
550
551 ENTRY("mri_3dalign_one") ;
552
553 fitim = basis->fitim ;
554 chol_fitim = basis->chol_fitim ;
555
556 ASSIF(th1,0.0f) ; /* initialize output params */
557 ASSIF(th2,0.0f) ;
558 ASSIF(th3,0.0f) ;
559 ASSIF(dx ,0.0f) ;
560 ASSIF(dy ,0.0f) ;
561 ASSIF(dz ,0.0f) ;
562
563 /* use original image if possible */
564
565 if( im->kind == MRI_float ) fim = im ;
566 else fim = mri_to_float( im ) ;
567
568 iter = 0 ;
569
570 THD_rota_method( regmode ) ;
571
572 /* convert displacement threshold from voxels to mm in each direction */
573
574 dxt = (im->dx != 0.0) ? (fabs(im->dx) * dxy_thresh) : dxy_thresh ;
575 dyt = (im->dy != 0.0) ? (fabs(im->dy) * dxy_thresh) : dxy_thresh ;
576 dzt = (im->dz != 0.0) ? (fabs(im->dz) * dxy_thresh) : dxy_thresh ;
577
578 if( NONZERO_INITVALS ){ /* 04 Sep 2000 */
579 fit = (float *) malloc(sizeof(float)*7) ;
580 fit[0] = 1.0 ;
581 fit[1] = init_dth1; fit[2] = init_dth2; fit[3] = init_dth3; /* degrees */
582 fit[4] = init_dx ; fit[5] = init_dy ; fit[6] = init_dz ; /* mm */
583
584 good = 1 ;
585 } else {
586
587 /* 06 Jun 2002: do initial fit with trimmed image, if ordered */
588
589 if( DOTRIM ){
590 tim = IMTRIM(fim) ;
591 fit = mri_delayed_lsqfit( tim , fitim , chol_fitim ) ;
592 mri_free( tim ) ; tim = NULL ;
593 } else { /* L2 fit input image */
594 fit = mri_delayed_lsqfit( fim , fitim , chol_fitim ) ;
595 }
596
597 good = ( 10.0*fabs(fit[4]) > dxt || 10.0*fabs(fit[5]) > dyt ||
598 10.0*fabs(fit[6]) > dzt || 10.0*fabs(fit[1]) > phi_thresh ||
599 10.0*fabs(fit[2]) > phi_thresh || 10.0*fabs(fit[3]) > phi_thresh ) ;
600 }
601
602 if( verbose )
603 fprintf(stderr,
604 "\nFirst fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
605 fit[0] , fit[1] , fit[2] , fit[3] , fit[4] , fit[5] , fit[6] ) ;
606
607 /*-- iterate fit --*/
608
609 while( good ){
610 tim = THD_rota3D( fim ,
611 ax1,fit[1]*DFAC , ax2,fit[2]*DFAC , ax3,fit[3]*DFAC ,
612 dcode , fit[4],fit[5],fit[6] ) ;
613
614 TRIM(tim) ; /* 06 Jun 2002: trim it if ordered to */
615
616 dfit = mri_delayed_lsqfit( tim , fitim , chol_fitim ) ; /* delta angle/shift */
617 mri_free( tim ) ; tim = NULL ;
618
619 fit[1] += dfit[1] ; fit[2] += dfit[2] ; fit[3] += dfit[3] ; /* accumulate */
620 fit[4] += dfit[4] ; fit[5] += dfit[5] ; fit[6] += dfit[6] ; /* angle/shift */
621
622 if( verbose ){
623 fprintf(stderr,
624 "Delta fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
625 dfit[0], dfit[1], dfit[2], dfit[3], dfit[4], dfit[5], dfit[6] ) ;
626 fprintf(stderr,
627 "Total fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
628 dfit[0], fit[1], fit[2], fit[3], fit[4], fit[5], fit[6] ) ;
629 }
630
631 good = (++iter < max_iter) &&
632 ( fabs(dfit[4]) > dxt || fabs(dfit[5]) > dyt ||
633 fabs(dfit[6]) > dzt || fabs(dfit[1]) > phi_thresh ||
634 fabs(dfit[2]) > phi_thresh || fabs(dfit[3]) > phi_thresh ) ;
635
636 free(dfit) ; dfit = NULL ;
637 } /* end while */
638
639 if( verbose ) fprintf(stderr,"Iteration complete at %d steps\n",iter) ;
640
641 /*-- save final alignment parameters --*/
642
643 ASSIF(th1,fit[1]*DFAC) ; /* convert to radians */
644 ASSIF(th2,fit[2]*DFAC) ;
645 ASSIF(th3,fit[3]*DFAC) ;
646 ASSIF(dx ,fit[4] ) ;
647 ASSIF(dy ,fit[5] ) ;
648 ASSIF(dz ,fit[6] ) ;
649
650 /*-- do the actual realignment --*/
651
652 tim = NULL ;
653 if( ! noreg ){
654 if( final_regmode < 0 ) final_regmode = regmode ; /* 20 Nov 1998 */
655 THD_rota_method( final_regmode ) ;
656
657 if( im->kind == MRI_complex ){ /* 20 May 2009: special case! */
658 MRI_IMARR *impair ; MRI_IMAGE *rim,*iim , *xim,*yim ;
659 impair = mri_complex_to_pair(im) ;
660 if( impair == NULL ){
661 ERROR_message("mri_complex_to_pair fails in mri_3dalign_one!") ;
662 } else {
663 rim = IMAGE_IN_IMARR(impair,0) ;
664 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
665 xim = THD_rota3D( rim ,
666 ax1,fit[1]*DFAC, ax2,fit[2]*DFAC, ax3,fit[3]*DFAC,
667 dcode , fit[4],fit[5],fit[6] ) ; mri_free(rim) ;
668 yim = THD_rota3D( iim ,
669 ax1,fit[1]*DFAC, ax2,fit[2]*DFAC, ax3,fit[3]*DFAC,
670 dcode , fit[4],fit[5],fit[6] ) ; mri_free(iim) ;
671 tim = mri_pair_to_complex(xim,yim) ; mri_free(xim) ; mri_free(yim) ;
672 }
673
674 } else { /* real-valued input */
675 tim = THD_rota3D( fim ,
676 ax1,fit[1]*DFAC , ax2,fit[2]*DFAC , ax3,fit[3]*DFAC ,
677 dcode , fit[4],fit[5],fit[6] ) ;
678 }
679 }
680
681 if( tim != NULL && tim->kind == MRI_float && clipit &&
682 (final_regmode == MRI_QUINTIC || final_regmode==MRI_CUBIC ||
683 final_regmode == MRI_HEPTIC || final_regmode==MRI_FOURIER ) ){
684
685 register int ii ;
686 register float ftop , fbot , *tar ;
687
688 ftop = mri_max( fim ); fbot = mri_min( fim );
689 tar = MRI_FLOAT_PTR(tim) ;
690 for( ii=0 ; ii < tim->nvox ; ii++ ){
691 if( tar[ii] < fbot ) tar[ii] = fbot ;
692 else if( tar[ii] > ftop ) tar[ii] = ftop ;
693 }
694 }
695
696 if( fim != im ) mri_free(fim) ; /* if it was a copy, junk it */
697
698 RETURN( tim ); /* 10-4, good buddy */
699 }
700
701 /*---------------------------------------------------------------------------*/
702 /* Register the first volume in imar to the base, then rotate/shift all
703 the others in the same way.
704 *//*-------------------------------------------------------------------------*/
705
mri_3dalign_oneplus(MRI_3dalign_basis * basis,MRI_IMARR * imar,float * th1,float * th2,float * th3,float * dx,float * dy,float * dz)706 MRI_IMARR * mri_3dalign_oneplus( MRI_3dalign_basis *basis, MRI_IMARR *imar ,
707 float *th1 , float *th2 , float *th3 ,
708 float *dx , float *dy , float *dz )
709 {
710 int nim = IMARR_COUNT(imar) , kk ;
711 MRI_IMAGE *bim , *outim ;
712 MRI_IMARR *outar ;
713 float dth1,dth2,dth3 , ddx,ddy,ddz ;
714
715 ENTRY("mri_3dalign_oneplus") ;
716
717 bim = IMARR_SUBIM(imar,0) ;
718 outim = mri_3dalign_one( basis , bim ,
719 &dth1,&dth2,&dth3 , &ddx,&ddy,&ddz ) ;
720
721 ASSIF(th1,dth1) ; ASSIF(th2,dth2) ; ASSIF(th3,dth3) ;
722 ASSIF(dx ,ddx ) ; ASSIF(dy ,ddy ) ; ASSIF(dz ,ddz ) ;
723
724 if( outim == NULL ) RETURN(NULL) ;
725
726 INIT_IMARR (outar) ;
727 ADDTO_IMARR(outar,outim) ;
728
729 for( kk=1 ; kk < nim ; kk++ ){
730 bim = IMARR_SUBIM(imar,kk) ;
731
732 if( bim->kind == MRI_complex ){ /* special case! */
733 MRI_IMARR *impair ; MRI_IMAGE *rim,*iim , *xim,*yim ;
734 impair = mri_complex_to_pair(bim) ;
735 if( impair == NULL ){
736 ERROR_message("mri_complex_to_pair fails in mri_3dalign_oneplus! ") ;
737 } else {
738 rim = IMAGE_IN_IMARR(impair,0) ;
739 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
740 xim = THD_rota3D( rim ,
741 ax1,dth1, ax2,dth2, ax3,dth3,
742 dcode , ddx,ddy,ddz ) ; mri_free(rim) ;
743 yim = THD_rota3D( iim ,
744 ax1,dth1, ax2,dth2, ax3,dth3,
745 dcode , ddx,ddy,ddz ) ; mri_free(iim) ;
746 outim = mri_pair_to_complex(xim,yim) ; mri_free(xim) ; mri_free(yim) ;
747 }
748
749 } else { /* real-valued input */
750 outim = THD_rota3D( bim ,
751 ax1,dth1, ax2,dth2, ax3,dth3, dcode , ddx,ddy,ddz ) ;
752
753 if( outim != NULL && outim->kind == MRI_float && clipit &&
754 (final_regmode == MRI_QUINTIC || final_regmode==MRI_CUBIC ||
755 final_regmode == MRI_HEPTIC || final_regmode==MRI_FOURIER ) ){
756
757 register int ii ; register float ftop, fbot, *tar ;
758 ftop = mri_max(bim); fbot = mri_min(bim); tar = MRI_FLOAT_PTR(outim);
759 for( ii=0 ; ii < outim->nvox ; ii++ ){
760 if( tar[ii] < fbot ) tar[ii] = fbot ;
761 else if( tar[ii] > ftop ) tar[ii] = ftop ;
762 }
763 }
764 }
765
766 ADDTO_IMARR(outar,outim) ;
767 }
768
769 RETURN(outar) ;
770 }
771
772
773 /*---------------------------------------------------------------------------*/
774 /* apply rotate/shift parameters to list of images 27 Oct 2016 [rickr]
775
776 Like mri_3dalign_oneplus(), but expect that mri_3dalign_one() has already
777 happened, so parameters are known.
778
779 If keep_datum and im->kind is MRI_byte or MRI_short, convert back.
780 21 Feb 2017 [rickr]
781 *//*-------------------------------------------------------------------------*/
782
mri_3dalign_apply(MRI_3dalign_basis * basis,MRI_IMARR * imar,float dth1,float dth2,float dth3,float ddx,float ddy,float ddz,int keep_datum)783 MRI_IMARR * mri_3dalign_apply( MRI_3dalign_basis *basis, MRI_IMARR *imar ,
784 float dth1 , float dth2 , float dth3 ,
785 float ddx , float ddy , float ddz ,
786 int keep_datum )
787 {
788 int nim = IMARR_COUNT(imar) , kk ;
789 MRI_IMAGE *bim=NULL , *outim=NULL , *keepim=NULL ;
790 MRI_IMARR *outar=NULL ;
791
792 ENTRY("mri_3dalign_apply") ;
793
794 INIT_IMARR (outar) ;
795
796 for( kk=0 ; kk < nim ; kk++ ){
797 bim = IMARR_SUBIM(imar,kk) ;
798
799 if( bim->kind == MRI_complex ){ /* special case! */
800 MRI_IMARR *impair ; MRI_IMAGE *rim,*iim , *xim,*yim ;
801 impair = mri_complex_to_pair(bim) ;
802 if( impair == NULL ){
803 ERROR_message("mri_complex_to_pair fails in mri_3dalign_apply! ") ;
804 } else {
805 rim = IMAGE_IN_IMARR(impair,0) ;
806 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
807 xim = THD_rota3D( rim ,
808 ax1,dth1, ax2,dth2, ax3,dth3,
809 dcode , ddx,ddy,ddz ) ; mri_free(rim) ;
810 yim = THD_rota3D( iim ,
811 ax1,dth1, ax2,dth2, ax3,dth3,
812 dcode , ddx,ddy,ddz ) ; mri_free(iim) ;
813 outim = mri_pair_to_complex(xim,yim) ; mri_free(xim) ; mri_free(yim) ;
814 }
815
816 } else { /* real-valued input */
817 outim = THD_rota3D( bim ,
818 ax1,dth1, ax2,dth2, ax3,dth3, dcode , ddx,ddy,ddz ) ;
819
820 /* if keeping shorts or bytes, convert immediately 21 Feb 2017 */
821 if( outim && keep_datum ) {
822 if( bim->kind == MRI_byte ) {
823 keepim = mri_to_byte(outim); mri_free(outim); outim=keepim;
824 } else if( bim->kind == MRI_short ) {
825 keepim = mri_to_short(1.0, outim); mri_free(outim); outim=keepim;
826 }
827 }
828
829 if( outim != NULL && outim->kind == MRI_float && clipit &&
830 (final_regmode == MRI_QUINTIC || final_regmode==MRI_CUBIC ||
831 final_regmode == MRI_HEPTIC || final_regmode==MRI_FOURIER ) ){
832
833 register int ii ; register float ftop, fbot, *tar ;
834 ftop = mri_max(bim); fbot = mri_min(bim); tar = MRI_FLOAT_PTR(outim);
835 for( ii=0 ; ii < outim->nvox ; ii++ ){
836 if( tar[ii] < fbot ) tar[ii] = fbot ;
837 else if( tar[ii] > ftop ) tar[ii] = ftop ;
838 }
839 }
840 }
841
842 ADDTO_IMARR(outar,outim) ;
843 }
844
845 RETURN(outar) ;
846 }
847
848 /*---------------------------------------------------------------------------*/
849
mri_3dalign_many(MRI_IMAGE * im,MRI_IMAGE * imwt,MRI_IMARR * ims,float * th1,float * th2,float * th3,float * dx,float * dy,float * dz)850 MRI_IMARR * mri_3dalign_many( MRI_IMAGE *im , MRI_IMAGE *imwt , MRI_IMARR *ims ,
851 float *th1 , float *th2 , float *th3 ,
852 float *dx , float *dy , float *dz )
853 {
854 int kim ;
855 MRI_IMAGE *tim ;
856 MRI_IMARR *alim ;
857 MRI_3dalign_basis *basis ;
858
859 ENTRY("mri_3dalign_many") ;
860
861 basis = mri_3dalign_setup( im , imwt ) ;
862 if( basis == NULL ) RETURN( NULL );
863
864 INIT_IMARR( alim ) ;
865
866 #define PK(x) ( (x!=NULL) ? (x+kim) : NULL )
867
868 for( kim=0 ; kim < ims->num ; kim++ ){
869 tim = mri_3dalign_one( basis , ims->imarr[kim] ,
870 PK(th1), PK(th2), PK(th3),
871 PK(dx) , PK(dy) , PK(dz) ) ;
872 ADDTO_IMARR(alim,tim) ;
873 }
874
875 mri_3dalign_cleanup( basis ) ;
876 RETURN( alim );
877 }
878
879 /*--------------------------------------------------------------------*/
880
mri_3dalign_cleanup(MRI_3dalign_basis * basis)881 void mri_3dalign_cleanup( MRI_3dalign_basis * basis )
882 {
883 ENTRY("mri_3dalign_cleanup") ;
884 if( basis == NULL ) EXRETURN ;
885
886 if( basis->fitim != NULL ){ DESTROY_IMARR( basis->fitim ) ; }
887 if( basis->chol_fitim != NULL ){ free(basis->chol_fitim) ; }
888
889 free(basis) ; EXRETURN ;
890 }
891