00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*** NOT 7D SAFE ***/
00010 
00011 #define REF_FLOAT_SINGLE
00012 #undef  VOX_SHORT
00013 #include "pcor.h"
00014 
00015 #define DFILT_SIGMA  (4.0*0.42466090)  /* 4.0 = FWHM */
00016 #define MAX_ITER     5
00017 #define DXY_THRESH   0.15         /* pixels */
00018 #define PHI_THRESH   0.45         /* degrees */
00019 #define DFAC         (PI/180.0)
00020 
00021 /** FINE_FIT:
00022     experimental code for a less smoothed fit at the end
00023     RW Cox, 17 July 1995
00024  **/
00025 
00026 #define FINE_FIT
00027 
00028 #ifdef FINE_FIT
00029 #  define FINE_SIGMA   (1.0*0.42466090)
00030 #  define FINE_DXY_THRESH  0.07
00031 #  define FINE_PHI_THRESH  0.21
00032 #endif
00033 
00034 #define USE_DELAYED_FIT
00035 
00036 /*********************************************************************
00037   05 Nov 1997: make the parameters that control the iterations
00038                be alterable.
00039 **********************************************************************/
00040 
00041 static float dfilt_sigma     = DFILT_SIGMA ,
00042              dxy_thresh      = DXY_THRESH ,
00043              phi_thresh      = PHI_THRESH ,
00044              fine_sigma      = FINE_SIGMA ,
00045              fine_dxy_thresh = FINE_DXY_THRESH ,
00046              fine_phi_thresh = FINE_PHI_THRESH  ;
00047 
00048 static int max_iter = MAX_ITER ;
00049 
00050 void mri_align_params( int maxite,
00051                        float sig , float dxy , float dph ,
00052                        float fsig, float fdxy, float fdph )
00053 {
00054    if( maxite > 0   ) max_iter    = maxite ; else max_iter    = MAX_ITER    ;
00055    if( sig    > 0.0 ) dfilt_sigma = sig    ; else dfilt_sigma = DFILT_SIGMA ;
00056    if( dxy    > 0.0 ) dxy_thresh  = dxy    ; else dxy_thresh  = DXY_THRESH  ;
00057    if( dph    > 0.0 ) phi_thresh  = dph    ; else phi_thresh  = PHI_THRESH  ;
00058 
00059    fine_sigma = fsig ;
00060    if( fdxy > 0.0 ) fine_dxy_thresh = fdxy ; else fine_dxy_thresh = FINE_DXY_THRESH ;
00061    if( fdph > 0.0 ) fine_phi_thresh = fdph ; else fine_phi_thresh = FINE_PHI_THRESH ;
00062 
00063    return ;
00064 }
00065 
00066 /*-------------------------------------------------------------------------------*/
00067 static int almode_coarse = MRI_BICUBIC ;  /* 1 Oct 1998 */
00068 static int almode_fine   = MRI_BICUBIC ;
00069 static int almode_reg    = MRI_BICUBIC ;
00070 
00071 #define MRI_ROTA_COARSE(a,b,c,d) mri_rota_variable(almode_coarse,(a),(b),(c),(d))
00072 #define MRI_ROTA_FINE(a,b,c,d)   mri_rota_variable(almode_fine  ,(a),(b),(c),(d))
00073 #define MRI_ROTA_REG(a,b,c,d)    mri_rota_variable(almode_reg   ,(a),(b),(c),(d))
00074 
00075 void mri_align_method( int coarse , int fine , int reg )  /* 1 Oct 1998 */
00076 {
00077    if( coarse > 0 ) almode_coarse = coarse ;
00078    if( fine   > 0 ) almode_fine   = fine   ;
00079    if( reg    > 0 ) almode_reg    = reg    ;
00080    return ;
00081 }
00082 /*-------------------------------------------------------------------------------*/
00083 
00084 /*---------------------------------------------------------------------
00085    Inputs:  imbase = base image that others will align to
00086             imwt   = image of weight factors to align to
00087                        (if NULL, will generate one internally)
00088             ims    = array of images to align to imbase
00089             code   = inclusive OR of ALIGN_* codes, which are:
00090                        ALIGN_NOITER_CODE   --> no iteration
00091                        ALIGN_VERBOSE_CODE  --> print progress reports
00092                        ALIGN_REGISTER_CODE --> return images registered
00093                        ALIGN_DEBUG_CODE    --> print debugging info
00094 
00095    Outputs: dx[i] , dy[i] , phi[i] , i=0..(ims->num - 1) =
00096             inputs to mri_rota that will bring ims->imarr[i] into
00097             alignment with imbase (dx,dy in pixels, phi in radians).
00098 
00099             if( (code & ALIGN_REGISTER_CODE) != 0 ) THEN
00100             the return value from this function will be an array of
00101             registered images;  otherwise, the return value is NULL.
00102 
00103    Method:  Iterated differential alignment.  Only works for small
00104             displacements (up to 4 pixels, say).
00105 -----------------------------------------------------------------------*/
00106 
00107 MRI_IMARR * mri_align_dfspace( MRI_IMAGE * imbase , MRI_IMAGE * imwt ,
00108                                MRI_IMARR * ims ,
00109                                int code , float dx[] , float dy[] , float phi[] )
00110 {
00111    MRI_IMAGE * im1 , *bim,*xim,*yim,*tim , *bim2 , * im2 , *imww ;
00112    MRI_IMARR * fitim ;
00113    float * fit , *tar,*xar,*yar , *dfit ;
00114    int nx,ny, ii,jj, joff, iter, good, kim, debug, verbose ;
00115    float hnx,hny ;
00116    double * chol_fitim=NULL ;
00117 
00118 #ifdef FINE_FIT
00119    MRI_IMARR * fine_fitim  =NULL ;
00120    MRI_IMAGE * fine_imww   =NULL ;
00121    double * chol_fine_fitim=NULL ;
00122    int use_fine_fit = (fine_sigma > 0.0) ;
00123 #endif
00124 
00125    if( ! MRI_IS_2D(imbase) ){
00126       fprintf(stderr,"\n*** mri_align_dfspace: cannot use nD images!\a\n") ;
00127       EXIT(1) ;
00128    }
00129 
00130    debug    = (code & ALIGN_DEBUG_CODE)    != 0 ;
00131    verbose  = (code & ALIGN_VERBOSE_CODE)  != 0 && !debug ;
00132 
00133    if( verbose ){printf("-- mri_align_dfspace");fflush(stdout);}
00134 
00135    /** create the fitting images **/
00136 
00137    if( debug ){
00138       printf("-- mri_align_dfspace: code=%d\n",code);
00139       if( imbase->name != NULL )
00140         printf("-- imbase name = %s\n",imbase->name);
00141    }
00142 
00143    im1 = mri_to_float( imbase ) ;
00144    nx  = im1->nx ;  hnx = 0.5 * nx ;
00145    ny  = im1->ny ;  hny = 0.5 * ny ;
00146 
00147    bim = mri_filt_fft( im1 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;  /* smooth */
00148    xim = mri_filt_fft( im1 , dfilt_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ;  /* d/dx */
00149    yim = mri_filt_fft( im1 , dfilt_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ;  /* d/dy */
00150 
00151    tim = mri_new( nx , ny , MRI_float ) ;    /* x * d/dy - y * d/dx */
00152    tar = mri_data_pointer( tim ) ;           /* which is d/d(theta) */
00153    xar = mri_data_pointer( xim ) ;
00154    yar = mri_data_pointer( yim ) ;
00155    for( jj=0 ; jj < ny ; jj++ ){
00156       joff = jj * nx ;
00157       for( ii=0 ; ii < nx ; ii++ ){
00158          tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00159                                 - (jj-hny) * xar[ii+joff] ) ;
00160       }
00161    }
00162    INIT_IMARR ( fitim ) ;
00163    ADDTO_IMARR( fitim , bim ) ;
00164    ADDTO_IMARR( fitim , xim ) ;
00165    ADDTO_IMARR( fitim , yim ) ;
00166    ADDTO_IMARR( fitim , tim ) ;
00167 
00168    if( imwt == NULL ) imww = mri_to_float( bim ) ;  /* 28 Oct 1996 */
00169    else               imww = mri_to_float( imwt ) ;
00170 
00171    tar = MRI_FLOAT_PTR(imww) ;
00172    for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;  /* 16 Nov 1998 */
00173 
00174 #ifdef USE_DELAYED_FIT
00175    chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
00176 #endif
00177 
00178    /*** create the FINE_FIT analog of the above, if required ***/
00179 #ifdef FINE_FIT
00180    if( use_fine_fit ){
00181      bim = mri_filt_fft( im1 , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;  /* smooth */
00182      xim = mri_filt_fft( im1 , fine_sigma , 1 , 0 , FILT_FFT_WRAPAROUND ) ;  /* d/dx */
00183      yim = mri_filt_fft( im1 , fine_sigma , 0 , 1 , FILT_FFT_WRAPAROUND ) ;  /* d/dy */
00184 
00185      tim = mri_new( nx , ny , MRI_float ) ;    /* x * d/dy - y * d/dx */
00186      tar = mri_data_pointer( tim ) ;           /* which is d/d(theta) */
00187      xar = mri_data_pointer( xim ) ;
00188      yar = mri_data_pointer( yim ) ;
00189      for( jj=0 ; jj < ny ; jj++ ){
00190         joff = jj * nx ;
00191         for( ii=0 ; ii < nx ; ii++ ){
00192            tar[ii+joff] = DFAC * (  (ii-hnx) * yar[ii+joff]
00193                                   - (jj-hny) * xar[ii+joff] ) ;
00194         }
00195      }
00196      INIT_IMARR ( fine_fitim ) ;
00197      ADDTO_IMARR( fine_fitim , bim ) ;
00198      ADDTO_IMARR( fine_fitim , xim ) ;
00199      ADDTO_IMARR( fine_fitim , yim ) ;
00200      ADDTO_IMARR( fine_fitim , tim ) ;
00201 
00202      if( imwt == NULL ) fine_imww = mri_to_float( bim ) ;  /* 03 Oct 1997 */
00203      else               fine_imww = mri_to_float( imwt ) ;
00204 
00205      tar = MRI_FLOAT_PTR(fine_imww) ;
00206      for( ii=0 ; ii < nx*ny ; ii++ ) tar[ii] = fabs(tar[ii]) ;
00207 
00208 #ifdef USE_DELAYED_FIT
00209      chol_fine_fitim = mri_startup_lsqfit( fine_fitim , fine_imww ) ;
00210 #endif
00211    }
00212 #endif  /* FINE_FIT */
00213 
00214    mri_free( im1 ) ;
00215 
00216    /** fit each image to the fitting images **/
00217 
00218    for( kim=0 ; kim < ims->num ; kim++ ){
00219 
00220       if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
00221 
00222       if( debug ) printf("-- Start image %d: %s\n",kim,
00223                          (ims->imarr[kim]->name==NULL)?" ":ims->imarr[kim]->name);
00224 
00225       im2  = mri_to_float( ims->imarr[kim] ) ;
00226       bim2 = mri_filt_fft( im2 , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00227 #ifdef USE_DELAYED_FIT
00228       fit  = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00229 #else
00230       fit  = mri_lsqfit( bim2 , fitim , imww ) ;
00231 #endif
00232       mri_free( bim2 ) ;
00233 
00234       if( debug ) printf("   fit = %13.6g %13.6g %13.6g\n",
00235                          fit[1],fit[2],fit[3] ) ;
00236 
00237       iter = 0 ;
00238       good = ( fabs(fit[1]) > dxy_thresh ||
00239                fabs(fit[2]) > dxy_thresh || fabs(fit[3]) > phi_thresh ) &&
00240              ( (code & ALIGN_NOITER_CODE) == 0 ) ;
00241 
00242       while( good ){
00243          tim  = MRI_ROTA_COARSE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00244          bim2 = mri_filt_fft( tim , dfilt_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00245 #ifdef USE_DELAYED_FIT
00246          dfit = mri_delayed_lsqfit( bim2 , fitim , chol_fitim ) ;
00247 #else
00248          dfit = mri_lsqfit( bim2 , fitim , imww ) ;
00249 #endif
00250          mri_free( bim2 ) ; mri_free( tim ) ;
00251 
00252          if( debug )
00253             printf(" Cdfit = %13.6g %13.6g %13.6g\n",
00254                    dfit[1],dfit[2],dfit[3] ) ;
00255 
00256          fit[1] += dfit[1] ;
00257          fit[2] += dfit[2] ;
00258          fit[3] += dfit[3] ;
00259 
00260          good = (++iter < max_iter) &&
00261                   ( fabs(dfit[1]) > dxy_thresh ||
00262                     fabs(dfit[2]) > dxy_thresh || fabs(dfit[3]) > phi_thresh ) ;
00263 
00264          free(dfit) ; dfit = NULL ;
00265       } /* end while */
00266 
00267       /*** perform fine adjustments (always use bicubic interpolation) ***/
00268 #ifdef FINE_FIT
00269       if( use_fine_fit && iter < max_iter && (code & ALIGN_NOITER_CODE) == 0 ){
00270          good = 1 ;
00271          while( good ){
00272             tim  = MRI_ROTA_FINE( im2 , fit[1] , fit[2] , fit[3]*DFAC ) ;
00273             bim2 = mri_filt_fft( tim , fine_sigma , 0 , 0 , FILT_FFT_WRAPAROUND ) ;
00274 #ifdef USE_DELAYED_FIT
00275             dfit = mri_delayed_lsqfit( bim2 , fine_fitim , chol_fine_fitim ) ;
00276 #else
00277             dfit = mri_lsqfit( bim2 , fine_fitim , fine_imww ) ;
00278 #endif
00279             mri_free( bim2 ) ; mri_free( tim ) ;
00280 
00281             if( debug )
00282                printf(" Fdfit = %13.6g %13.6g %13.6g\n",
00283                       dfit[1],dfit[2],dfit[3] ) ;
00284 
00285             fit[1] += dfit[1] ;
00286             fit[2] += dfit[2] ;
00287             fit[3] += dfit[3] ;
00288 
00289             good = (++iter < max_iter) &&
00290                      ( fabs(dfit[1]) > fine_dxy_thresh ||
00291                        fabs(dfit[2]) > fine_dxy_thresh || fabs(dfit[3]) > fine_phi_thresh ) ;
00292 
00293             free(dfit) ; dfit = NULL ;
00294          } /* end while */
00295       }
00296 #endif
00297 
00298       dx[kim]  = fit[1] ;
00299       dy[kim]  = fit[2] ;
00300       phi[kim] = fit[3]*DFAC ;
00301 
00302       if( debug )
00303          printf("   FIT = %13.6g %13.6g %13.6g\n",
00304                 fit[1],fit[2],fit[3] ) ;
00305 
00306       free(fit) ; fit = NULL ; mri_free( im2 ) ;
00307    }
00308 
00309    if( verbose ){printf("\n");fflush(stdout);}
00310 
00311    DESTROY_IMARR( fitim ) ;
00312 #ifdef USE_DELAYED_FIT
00313    free(chol_fitim) ; chol_fitim = NULL ;
00314 #endif
00315 
00316 #ifdef FINE_FIT
00317    if( use_fine_fit ){
00318      DESTROY_IMARR( fine_fitim ) ;
00319      mri_free( fine_imww ) ;
00320 #ifdef USE_DELAYED_FIT
00321      free(chol_fine_fitim) ; chol_fine_fitim = NULL ;
00322 #endif
00323    }
00324 #endif
00325 
00326    mri_free( imww ) ;
00327 
00328    if( (code & ALIGN_REGISTER_CODE) == 0 ) return NULL ;
00329 
00330    /** do the actual registration, if ordered (is always bicubic) **/
00331 
00332    INIT_IMARR( fitim ) ;
00333 
00334    if( verbose ){printf("-- registering");fflush(stdout);}
00335 
00336    for( kim=0 ; kim < ims->num ; kim++ ){
00337       tim = MRI_ROTA_REG( ims->imarr[kim] , dx[kim],dy[kim],phi[kim] ) ;
00338       ADDTO_IMARR( fitim , tim ) ;
00339       if( verbose && kim%5 == 0 ){printf(".");fflush(stdout);}
00340    }
00341 
00342    if( verbose ){printf("\n");fflush(stdout);}
00343 
00344    return fitim ;
00345 }
00346 
00347 #ifdef ALLOW_DFTIME
00348 /*----------------------------------------------------------------------------
00349   Similar routine to above, but uses temporal filtering to determine how
00350   much to alter each pixel.  Note that if ALIGN_REGISTER_CODE is not turned
00351   on, this routine will produce identical outputs to the dfspace routine.
00352 ------------------------------------------------------------------------------*/
00353 
00354 #define NREF 5  /* 1st 3 refs = dx,dy,phi; others are polynomials */
00355 
00356 MRI_IMARR * mri_align_dftime( MRI_IMAGE * imbase , MRI_IMAGE * imwt ,
00357                               MRI_IMARR * ims ,
00358                               int code , float dx[] , float dy[] , float phi[] )
00359 {
00360    int kim , ii , nim,npix , my_code , nx,ny , jj , detrend , doboth , verbose,debug ;
00361    references * xypref ;
00362    voxel_corr * xypcor ;
00363    float ref[NREF] , val , sum ;
00364    MRI_IMAGE * tim ;
00365    MRI_IMARR * fitim , * outim , * lims ;
00366    float * fitar[NREF] , * tar ;
00367 
00368    /** run the dfspace code to get the displacements **/
00369 
00370    debug   = (code & ALIGN_DEBUG_CODE) != 0 ;
00371    verbose = (code & ALIGN_VERBOSE_CODE) != 0 && !debug ;
00372 
00373    doboth  = (code & ALIGN_DOBOTH_CODE) != 0  && (code & ALIGN_REGISTER_CODE) != 0 ;
00374 
00375    if( !doboth ){
00376       my_code = code & ~ALIGN_REGISTER_CODE ;
00377    } else {
00378       my_code = code ;
00379    }
00380    lims = mri_align_dfspace( imbase , imwt , ims , my_code , dx,dy,phi ) ;
00381    if( (code & ALIGN_REGISTER_CODE) == 0 ) return NULL ;
00382 
00383    /** do some real work now **/
00384 
00385    detrend = (code & ALIGN_DETREND_CODE) != 0 ;
00386    if( !doboth ) lims = ims ;
00387 
00388    /** run the time series code to get the pixel-by-pixel
00389        correlation between the displacements and the pixel intensities **/
00390 
00391    nim  = lims->num ;
00392    nx   = imbase->nx ; ny = imbase->ny ; npix = nx * ny ;
00393 
00394    xypref = new_references( NREF ) ;
00395    xypcor = new_voxel_corr( npix , NREF ) ;
00396 
00397    if( verbose ){printf("-- temporal filtering");fflush(stdout);}
00398 
00399    for( kim=0 ; kim < nim ; kim++ ){
00400       ref[0] = dx[kim] ; ref[1] = dy[kim] ; ref[2] = phi[kim] / DFAC ;
00401 
00402       val = 1.0 ;
00403       for( jj=3 ; jj < NREF ; jj++ ){
00404          ref[jj] = val ;
00405          val    *= (kim-0.5*nim) ;
00406       }
00407 
00408       tim = lims->imarr[kim] ;
00409       if( tim->kind != MRI_float ) tim = mri_to_float( lims->imarr[kim] ) ;
00410 
00411       update_references( ref , xypref ) ;
00412       update_voxel_corr( MRI_FLOAT_PTR(tim) , xypref , xypcor ) ;
00413 
00414       if( tim != lims->imarr[kim] ) mri_free( tim ) ;
00415 
00416       if( verbose && kim%10==0 ){printf(".");fflush(stdout);}
00417    }
00418 
00419    INIT_IMARR( fitim ) ;
00420    for( jj=0 ; jj < 3 ; jj++ ){
00421       tim = mri_new( nx , ny , MRI_float ) ;
00422       ADDTO_IMARR( fitim , tim ) ;
00423       fitar[jj] = MRI_FLOAT_PTR(tim) ;
00424    }
00425 
00426    if( detrend ){
00427       for( jj=3 ; jj < NREF ; jj++ ){
00428          tim = mri_new( nx , ny , MRI_float ) ;
00429          ADDTO_IMARR( fitim , tim ) ;
00430          fitar[jj] = MRI_FLOAT_PTR(tim) ;
00431       }
00432    } else {
00433       for( jj=3 ; jj < NREF ; jj++ ) fitar[jj] = NULL ;
00434    }
00435 
00436    if( verbose ){printf("!");fflush(stdout);}
00437 
00438    get_lsqfit( xypref , xypcor , fitar ) ;
00439 
00440    free_references( xypref ) ;
00441    free_voxel_corr( xypcor ) ;
00442 
00443    /** apply the fits to each image **/
00444 
00445    INIT_IMARR( outim ) ;
00446 
00447    for( kim=0 ; kim < nim ; kim++ ){
00448 
00449       tim = mri_to_float( lims->imarr[kim] ) ;
00450       tar = MRI_FLOAT_PTR(tim) ;
00451 
00452       if( detrend ){
00453          for( ii=0 ; ii < npix ; ii++ ){
00454             sum = tar[ii] - (  fitar[0][ii] * dx[kim]
00455                              + fitar[1][ii] * dy[kim] + fitar[2][ii] * phi[kim] / DFAC ) ;
00456             for( jj=3 ; jj < NREF ; jj++ ) sum -= fitar[jj][ii] * ref[jj] ;
00457             tar[ii] = sum ;
00458          }
00459       } else {
00460          for( ii=0 ; ii < npix ; ii++ )
00461             tar[ii] -= (  fitar[0][ii] * dx[kim]
00462                         + fitar[1][ii] * dy[kim] + fitar[2][ii] * phi[kim] / DFAC ) ;
00463       }
00464 
00465       ADDTO_IMARR( outim , tim ) ;
00466 
00467       if( verbose && kim%10==0 ){printf(".");fflush(stdout);}
00468 
00469       if( lims != ims ) mri_free( lims->imarr[kim] ) ;
00470    }
00471 
00472    if( verbose ){printf("\n");fflush(stdout);}
00473 
00474    DESTROY_IMARR( fitim ) ;
00475    if( lims != ims ) FREE_IMARR( lims ) ;
00476 
00477    return outim ;
00478 }
00479 #endif /* ALLOW_DFTIME */