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 /*** NOT 7D SAFE ***/
10
11 #define REF_FLOAT_SINGLE
12 #undef VOX_SHORT
13 #include "pcor.h"
14
15 #define DFILT_SIGMA (4.0*0.42466090) /* 4.0 = FWHM */
16 #define MAX_ITER 5
17 #define DXY_THRESH 0.15 /* pixels */
18 #define PHI_THRESH 0.45 /* degrees */
19 #define DFAC (PI/180.0)
20
21 /** FINE_FIT:
22 experimental code for a less smoothed fit at the end
23 RW Cox, 17 July 1995
24 **/
25
26 #define FINE_FIT
27
28 #ifdef FINE_FIT
29 # define FINE_SIGMA (1.0*0.42466090)
30 # define FINE_DXY_THRESH 0.07
31 # define FINE_PHI_THRESH 0.21
32 #endif
33
34 #define USE_DELAYED_FIT
35
36 /*********************************************************************
37 05 Nov 1997: make the parameters that control the iterations
38 be alterable.
39 **********************************************************************/
40
41 static float dfilt_sigma = DFILT_SIGMA ,
42 dxy_thresh = DXY_THRESH ,
43 phi_thresh = PHI_THRESH ,
44 fine_sigma = FINE_SIGMA ,
45 fine_dxy_thresh = FINE_DXY_THRESH ,
46 fine_phi_thresh = FINE_PHI_THRESH ;
47
48 static int max_iter = MAX_ITER ;
49
mri_align_params(int maxite,float sig,float dxy,float dph,float fsig,float fdxy,float fdph)50 void mri_align_params( int maxite,
51 float sig , float dxy , float dph ,
52 float fsig, float fdxy, float fdph )
53 {
54 if( maxite > 0 ) max_iter = maxite ; else max_iter = MAX_ITER ;
55 if( sig > 0.0 ) dfilt_sigma = sig ; else dfilt_sigma = DFILT_SIGMA ;
56 if( dxy > 0.0 ) dxy_thresh = dxy ; else dxy_thresh = DXY_THRESH ;
57 if( dph > 0.0 ) phi_thresh = dph ; else phi_thresh = PHI_THRESH ;
58
59 fine_sigma = fsig ;
60 if( fdxy > 0.0 ) fine_dxy_thresh = fdxy ; else fine_dxy_thresh = FINE_DXY_THRESH ;
61 if( fdph > 0.0 ) fine_phi_thresh = fdph ; else fine_phi_thresh = FINE_PHI_THRESH ;
62
63 return ;
64 }
65
66 /*-------------------------------------------------------------------------------*/
67 static int almode_coarse = MRI_BICUBIC ; /* 1 Oct 1998 */
68 static int almode_fine = MRI_BICUBIC ;
69 static int almode_reg = MRI_BICUBIC ;
70
71 #define MRI_ROTA_COARSE(a,b,c,d) mri_rota_variable(almode_coarse,(a),(b),(c),(d))
72 #define MRI_ROTA_FINE(a,b,c,d) mri_rota_variable(almode_fine ,(a),(b),(c),(d))
73 #define MRI_ROTA_REG(a,b,c,d) mri_rota_variable(almode_reg ,(a),(b),(c),(d))
74
mri_align_method(int coarse,int fine,int reg)75 void mri_align_method( int coarse , int fine , int reg ) /* 1 Oct 1998 */
76 {
77 if( coarse > 0 ) almode_coarse = coarse ;
78 if( fine > 0 ) almode_fine = fine ;
79 if( reg > 0 ) almode_reg = reg ;
80 return ;
81 }
82 /*-------------------------------------------------------------------------------*/
83
84 /*---------------------------------------------------------------------
85 Inputs: imbase = base image that others will align to
86 imwt = image of weight factors to align to
87 (if NULL, will generate one internally)
88 ims = array of images to align to imbase
89 code = inclusive OR of ALIGN_* codes, which are:
90 ALIGN_NOITER_CODE --> no iteration
91 ALIGN_VERBOSE_CODE --> print progress reports
92 ALIGN_REGISTER_CODE --> return images registered
93 ALIGN_DEBUG_CODE --> print debugging info
94
95 Outputs: dx[i] , dy[i] , phi[i] , i=0..(ims->num - 1) =
96 inputs to mri_rota that will bring ims->imarr[i] into
97 alignment with imbase (dx,dy in pixels, phi in radians).
98
99 if( (code & ALIGN_REGISTER_CODE) != 0 ) THEN
100 the return value from this function will be an array of
101 registered images; otherwise, the return value is NULL.
102
103 Method: Iterated differential alignment. Only works for small
104 displacements (up to 4 pixels, say).
105 -----------------------------------------------------------------------*/
106
mri_align_dfspace(MRI_IMAGE * imbase,MRI_IMAGE * imwt,MRI_IMARR * ims,int code,float dx[],float dy[],float phi[])107 MRI_IMARR * mri_align_dfspace( MRI_IMAGE * imbase , MRI_IMAGE * imwt ,
108 MRI_IMARR * ims ,
109 int code , float dx[] , float dy[] , float phi[] )
110 {
111 MRI_IMAGE * im1 , *bim,*xim,*yim,*tim , *bim2 , * im2 , *imww ;
112 MRI_IMARR * fitim ;
113 float * fit , *tar,*xar,*yar , *dfit ;
114 int nx,ny, ii,jj, joff, iter, good, kim, debug, verbose ;
115 float hnx,hny ;
116 double * chol_fitim=NULL ;
117
118 #ifdef FINE_FIT
119 MRI_IMARR * fine_fitim =NULL ;
120 MRI_IMAGE * fine_imww =NULL ;
121 double * chol_fine_fitim=NULL ;
122 int use_fine_fit = (fine_sigma > 0.0) ;
123 #endif
124
125 if( ! MRI_IS_2D(imbase) ){
126 fprintf(stderr,"\n*** mri_align_dfspace: cannot use nD images!\a\n") ;
127 EXIT(1) ;
128 }
129
130 debug = (code & ALIGN_DEBUG_CODE) != 0 ;
131 verbose = (code & ALIGN_VERBOSE_CODE) != 0 && !debug ;
132
133 if( verbose ){printf("-- mri_align_dfspace");fflush(stdout);}
134
135 /** create the fitting images **/
136
137 if( debug ){
138 printf("-- mri_align_dfspace: code=%d\n",code);
139 if( imbase->name != NULL )
140 printf("-- imbase name = %s\n",imbase->name);
141 }
142
143 im1 = mri_to_float( imbase ) ;
144 nx = im1->nx ; hnx = 0.5 * nx ;
145 ny = im1->ny ; hny = 0.5 * ny ;
146
147 bim = mri_filt_fft( im1 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ; /* smooth */
148 xim = mri_filt_fft( im1 , dfilt_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ; /* d/dx */
149 yim = mri_filt_fft( im1 , dfilt_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ; /* d/dy */
150
151 tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */
152 tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */
153 xar = mri_data_pointer( xim ) ;
154 yar = mri_data_pointer( yim ) ;
155 for( jj=0 ; jj < ny ; jj++ ){
156 joff = jj * nx ;
157 for( ii=0 ; ii < nx ; ii++ ){
158 tar[ii+joff] = DFAC * ( (ii-hnx) * yar[ii+joff]
159 - (jj-hny) * xar[ii+joff] ) ;
160 }
161 }
162 INIT_IMARR ( fitim ) ;
163 ADDTO_IMARR( fitim , bim ) ;
164 ADDTO_IMARR( fitim , xim ) ;
165 ADDTO_IMARR( fitim , yim ) ;
166 ADDTO_IMARR( fitim , tim ) ;
167
168 if( imwt == NULL ) imww = mri_to_float( bim ) ; /* 28 Oct 1996 */
169 else imww = mri_to_float( imwt ) ;
170
171 tar = MRI_FLOAT_PTR(imww) ;
172 for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ; /* 16 Nov 1998 */
173
174 #ifdef USE_DELAYED_FIT
175 chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
176 #endif
177
178 /*** create the FINE_FIT analog of the above, if required ***/
179 #ifdef FINE_FIT
180 if( use_fine_fit ){
181 bim = mri_filt_fft( im1 , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ; /* smooth */
182 xim = mri_filt_fft( im1 , fine_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ; /* d/dx */
183 yim = mri_filt_fft( im1 , fine_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ; /* d/dy */
184
185 tim = mri_new( nx , ny , MRI_float ) ; /* x * d/dy - y * d/dx */
186 tar = mri_data_pointer( tim ) ; /* which is d/d(theta) */
187 xar = mri_data_pointer( xim ) ;
188 yar = mri_data_pointer( yim ) ;
189 for( jj=0 ; jj < ny ; jj++ ){
190 joff = jj * nx ;
191 for( ii=0 ; ii < nx ; ii++ ){
192 tar[ii+joff] = DFAC * ( (ii-hnx) * yar[ii+joff]
193 - (jj-hny) * xar[ii+joff] ) ;
194 }
195 }
196 INIT_IMARR ( fine_fitim ) ;
197 ADDTO_IMARR( fine_fitim , bim ) ;
198 ADDTO_IMARR( fine_fitim , xim ) ;
199 ADDTO_IMARR( fine_fitim , yim ) ;
200 ADDTO_IMARR( fine_fitim , tim ) ;
201
202 if( imwt == NULL ) fine_imww = mri_to_float( bim ) ; /* 03 Oct 1997 */
203 else fine_imww = mri_to_float( imwt ) ;
204
205 tar = MRI_FLOAT_PTR(fine_imww) ;
206 for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;
207
208 #ifdef USE_DELAYED_FIT
209 chol_fine_fitim = mri_startup_lsqfit( fine_fitim , fine_imww ) ;
210 #endif
211 }
212 #endif /* FINE_FIT */
213
214 mri_free( im1 ) ;
215
216 /** fit each image to the fitting images **/
217
218 for( kim=0 ; kim < ims->num ; kim++ ){
219
220 if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
221
222 if( debug ) printf("-- Start image %d: %s\n",kim,
223 (ims->imarr[kim]->name==NULL)?" ":ims->imarr[kim]->name);
224
225 im2 = mri_to_float( ims->imarr[kim] ) ;
226 bim2 = mri_filt_fft( im2 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
227 #ifdef USE_DELAYED_FIT
228 fit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
229 #else
230 fit = mri_lsqfit( bim2 , fitim , imww ) ;
231 #endif
232 mri_free( bim2 ) ;
233
234
235 if( debug ) printf(" fit = %13.6g %13.6g %13.6g\n",
236 fit[1],fit[2],fit[3] ) ;
237
238 iter = 0 ;
239 good = ( fabs(fit[1]) > dxy_thresh ||
240 fabs(fit[2]) > dxy_thresh || fabs(fit[3]) > phi_thresh ) &&
241 ( (code & ALIGN_NOITER_CODE) == 0 ) ;
242
243 while( good ){
244 tim = MRI_ROTA_COARSE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
245 bim2 = mri_filt_fft( tim , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
246 /* fprintf(stderr,"nx=%d, ny=%d\n", bim2->nx, bim2->ny); image is returned square */
247
248 #ifdef USE_DELAYED_FIT
249 dfit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
250 #else
251 dfit = mri_lsqfit( bim2 , fitim , imww ) ;
252 #endif
253 mri_free( bim2 ) ; mri_free( tim ) ;
254
255 if( debug )
256 printf(" Cdfit = %13.6g %13.6g %13.6g\n",
257 dfit[1],dfit[2],dfit[3] ) ;
258
259 fit[1] += dfit[1] ;
260 fit[2] += dfit[2] ;
261 fit[3] += dfit[3] ;
262
263 good = (++iter < max_iter) &&
264 ( fabs(dfit[1]) > dxy_thresh ||
265 fabs(dfit[2]) > dxy_thresh || fabs(dfit[3]) > phi_thresh ) ;
266
267 free(dfit) ; dfit = NULL ;
268 } /* end while */
269
270 /*** perform fine adjustments (always use bicubic interpolation) ***/
271 #ifdef FINE_FIT
272 if( use_fine_fit && iter < max_iter && (code & ALIGN_NOITER_CODE) == 0 ){
273 good = 1 ;
274 while( good ){
275 tim = MRI_ROTA_FINE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
276 bim2 = mri_filt_fft( tim , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
277 #ifdef USE_DELAYED_FIT
278 dfit = mri_delayed_lsqfit( bim2 , fine_fitim , chol_fine_fitim ) ;
279 #else
280 dfit = mri_lsqfit( bim2 , fine_fitim , fine_imww ) ;
281 #endif
282 mri_free( bim2 ) ; mri_free( tim ) ;
283
284 if( debug )
285 printf(" Fdfit = %13.6g %13.6g %13.6g\n",
286 dfit[1],dfit[2],dfit[3] ) ;
287
288 fit[1] += dfit[1] ;
289 fit[2] += dfit[2] ;
290 fit[3] += dfit[3] ;
291
292 good = (++iter < max_iter) &&
293 ( fabs(dfit[1]) > fine_dxy_thresh ||
294 fabs(dfit[2]) > fine_dxy_thresh || fabs(dfit[3]) > fine_phi_thresh ) ;
295
296 free(dfit) ; dfit = NULL ;
297 } /* end while */
298 }
299 #endif
300
301 dx[kim] = fit[1] ;
302 dy[kim] = fit[2] ;
303 phi[kim] = fit[3]*DFAC ;
304
305 if( debug )
306 printf(" FIT = %13.6g %13.6g %13.6g\n",
307 fit[1],fit[2],fit[3] ) ;
308
309 free(fit) ; fit = NULL ; mri_free( im2 ) ;
310 }
311
312 if( verbose ){printf("\n");fflush(stdout);}
313
314 DESTROY_IMARR( fitim ) ;
315 #ifdef USE_DELAYED_FIT
316 free(chol_fitim) ; chol_fitim = NULL ;
317 #endif
318
319 #ifdef FINE_FIT
320 if( use_fine_fit ){
321 DESTROY_IMARR( fine_fitim ) ;
322 mri_free( fine_imww ) ;
323 #ifdef USE_DELAYED_FIT
324 free(chol_fine_fitim) ; chol_fine_fitim = NULL ;
325 #endif
326 }
327 #endif
328
329 mri_free( imww ) ;
330
331 if( (code & ALIGN_REGISTER_CODE) == 0 ) return NULL ;
332
333 /** do the actual registration, if ordered (is always bicubic) **/
334
335 INIT_IMARR( fitim ) ;
336
337 if( verbose ){printf("-- registering");fflush(stdout);}
338
339 for( kim=0 ; kim < ims->num ; kim++ ){
340 tim = MRI_ROTA_REG( ims->imarr[kim] , dx[kim],dy[kim],phi[kim] ) ;
341 ADDTO_IMARR( fitim , tim ) ;
342 if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
343 }
344
345 if( verbose ){printf("\n");fflush(stdout);}
346
347 return fitim ;
348 }
349
350 #ifdef ALLOW_DFTIME
351 /*----------------------------------------------------------------------------
352 Similar routine to above, but uses temporal filtering to determine how
353 much to alter each pixel. Note that if ALIGN_REGISTER_CODE is not turned
354 on, this routine will produce identical outputs to the dfspace routine.
355 ------------------------------------------------------------------------------*/
356
357 #define NREF 5 /* 1st 3 refs = dx,dy,phi; others are polynomials */
358
mri_align_dftime(MRI_IMAGE * imbase,MRI_IMAGE * imwt,MRI_IMARR * ims,int code,float dx[],float dy[],float phi[])359 MRI_IMARR * mri_align_dftime( MRI_IMAGE * imbase , MRI_IMAGE * imwt ,
360 MRI_IMARR * ims ,
361 int code , float dx[] , float dy[] , float phi[] )
362 {
363 int kim , ii , nim,npix , my_code , nx,ny , jj , detrend , doboth , verbose,debug ;
364 references * xypref ;
365 voxel_corr * xypcor ;
366 float ref[NREF] , val , sum ;
367 MRI_IMAGE * tim ;
368 MRI_IMARR * fitim , * outim , * lims ;
369 float * fitar[NREF] , * tar ;
370
371 /** run the dfspace code to get the displacements **/
372
373 debug = (code & ALIGN_DEBUG_CODE) != 0 ;
374 verbose = (code & ALIGN_VERBOSE_CODE) != 0 && !debug ;
375
376 doboth = (code & ALIGN_DOBOTH_CODE) != 0 && (code & ALIGN_REGISTER_CODE) != 0 ;
377
378 if( !doboth ){
379 my_code = code & ~ALIGN_REGISTER_CODE ;
380 } else {
381 my_code = code ;
382 }
383 lims = mri_align_dfspace( imbase , imwt , ims , my_code , dx,dy,phi ) ;
384 if( (code & ALIGN_REGISTER_CODE) == 0 ) return NULL ;
385
386 /** do some real work now **/
387
388 detrend = (code & ALIGN_DETREND_CODE) != 0 ;
389 if( !doboth ) lims = ims ;
390
391 /** run the time series code to get the pixel-by-pixel
392 correlation between the displacements and the pixel intensities **/
393
394 nim = lims->num ;
395 nx = imbase->nx ; ny = imbase->ny ; npix = nx * ny ;
396
397 xypref = new_references( NREF ) ;
398 xypcor = new_voxel_corr( npix , NREF ) ;
399
400 if( verbose ){printf("-- temporal filtering");fflush(stdout);}
401
402 for( kim=0 ; kim < nim ; kim++ ){
403 ref[0] = dx[kim] ; ref[1] = dy[kim] ; ref[2] = phi[kim] / DFAC ;
404
405 val = 1.0 ;
406 for( jj=3 ; jj < NREF ; jj++ ){
407 ref[jj] = val ;
408 val *= (kim-0.5*nim) ;
409 }
410
411 tim = lims->imarr[kim] ;
412 if( tim->kind != MRI_float ) tim = mri_to_float( lims->imarr[kim] ) ;
413
414 update_references( ref , xypref ) ;
415 update_voxel_corr( MRI_FLOAT_PTR(tim) , xypref , xypcor ) ;
416
417 if( tim != lims->imarr[kim] ) mri_free( tim ) ;
418
419 if( verbose && kim%10==0 ){printf(".");fflush(stdout);}
420 }
421
422 INIT_IMARR( fitim ) ;
423 for( jj=0 ; jj < 3 ; jj++ ){
424 tim = mri_new( nx , ny , MRI_float ) ;
425 ADDTO_IMARR( fitim , tim ) ;
426 fitar[jj] = MRI_FLOAT_PTR(tim) ;
427 }
428
429 if( detrend ){
430 for( jj=3 ; jj < NREF ; jj++ ){
431 tim = mri_new( nx , ny , MRI_float ) ;
432 ADDTO_IMARR( fitim , tim ) ;
433 fitar[jj] = MRI_FLOAT_PTR(tim) ;
434 }
435 } else {
436 for( jj=3 ; jj < NREF ; jj++ ) fitar[jj] = NULL ;
437 }
438
439 if( verbose ){printf("!");fflush(stdout);}
440
441 get_lsqfit( xypref , xypcor , fitar ) ;
442
443 free_references( xypref ) ;
444 free_voxel_corr( xypcor ) ;
445
446 /** apply the fits to each image **/
447
448 INIT_IMARR( outim ) ;
449
450 for( kim=0 ; kim < nim ; kim++ ){
451
452 tim = mri_to_float( lims->imarr[kim] ) ;
453 tar = MRI_FLOAT_PTR(tim) ;
454
455 if( detrend ){
456 for( ii=0 ; ii < npix ; ii++ ){
457 sum = tar[ii] - ( fitar[0][ii] * dx[kim]
458 + fitar[1][ii] * dy[kim] + fitar[2][ii] * phi[kim] / DFAC ) ;
459 for( jj=3 ; jj < NREF ; jj++ ) sum -= fitar[jj][ii] * ref[jj] ;
460 tar[ii] = sum ;
461 }
462 } else {
463 for( ii=0 ; ii < npix ; ii++ )
464 tar[ii] -= ( fitar[0][ii] * dx[kim]
465 + fitar[1][ii] * dy[kim] + fitar[2][ii] * phi[kim] / DFAC ) ;
466 }
467
468 ADDTO_IMARR( outim , tim ) ;
469
470 if( verbose && kim%10==0 ){printf(".");fflush(stdout);}
471
472 if( lims != ims ) mri_free( lims->imarr[kim] ) ;
473 }
474
475 if( verbose ){printf("\n");fflush(stdout);}
476
477 DESTROY_IMARR( fitim ) ;
478 if( lims != ims ) FREE_IMARR( lims ) ;
479
480 return outim ;
481 }
482 #endif /* ALLOW_DFTIME */
483