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 /*************************************************************************
00012  **       3D version of mri_2dalign.c -- RWCox -- October 1998          **
00013  *************************************************************************/
00014 
00015 #define MAX_ITER     5
00016 #define DXY_THRESH   0.07         /* pixels */
00017 #define PHI_THRESH   0.21         /* degrees */
00018 #define DFAC         (PI/180.0)
00019 
00020 static float dxy_thresh = DXY_THRESH ,
00021              phi_thresh = PHI_THRESH ,
00022              delfac     = 1.5 ;
00023 
00024 static int max_iter = MAX_ITER ;
00025 static int ax1 = 0 , ax2 = 1 , ax3 = 2 ;
00026 static int dcode = -1 ;
00027 
00028 static int wproc = 0 ;   /* 06 Jun 2002: process imwt? */
00029 static int wtrim = 0 ;   /* 06 Jun 2002: trimming stuff */
00030 
00031 static float sinit = 1.0 ; /* 22 Mar 2004: init scale? */
00032 
00033 #define DOTRIM (basis->xa >= 0)
00034 
00035 #define IMTRIM(qqq) mri_cut_3D( qqq , basis->xa,basis->xb ,  \
00036                                       basis->ya,basis->yb ,  \
00037                                       basis->za,basis->zb  )
00038 
00039 /*! Macro to replace an image with a trimmed copy, if needed. */
00040 
00041 #define TRIM(imq)                      \
00042   do{ if( DOTRIM ){                    \
00043         MRI_IMAGE *qim = IMTRIM(imq) ; \
00044         mri_free(imq) ; imq = qim ;    \
00045   } } while(0)
00046 
00047 /*! Macro to print out a volume range in verbose mode. */
00048 
00049 #define VRANG(str,imq)                            \
00050  do{ if(verbose){                                 \
00051        double tt=mri_max(imq), bb=mri_min(imq) ;  \
00052        fprintf(stderr,"  %s range: min=%g  max=%g\n",str,bb,tt); } } while(0)
00053 
00054 /*--------------------------------------------------------------------*/
00055 
00056 void mri_3dalign_wtrimming( int ttt ){ wtrim = ttt; } /* 06 Jun 2002 */
00057 
00058 void mri_3dalign_wproccing( int ttt ){ wproc = ttt; } /* 06 Jun 2002 */
00059 
00060 void mri_3dalign_scaleinit( float ttt )               /* 22 Mar 2004 */
00061 {
00062   if( ttt > 0.0 ) sinit = ttt ;
00063 }
00064 
00065 /*--------------------------------------------------------------------*/
00066 
00067 void mri_3dalign_params( int maxite ,
00068                          float dxy , float dph , float dfac ,
00069                          int bx1 , int bx2 , int bx3 , int dc )
00070 {
00071    if( maxite > 0   ) max_iter    = maxite ; else max_iter    = MAX_ITER    ;
00072    if( dxy    > 0.0 ) dxy_thresh  = dxy    ; else dxy_thresh  = DXY_THRESH  ;
00073    if( dph    > 0.0 ) phi_thresh  = dph    ; else phi_thresh  = PHI_THRESH  ;
00074    if( dfac   > 0.0 ) delfac      = dfac   ;
00075 
00076    if( bx1 >= 0 && bx1 <= 2 ) ax1 = bx1 ;
00077    if( bx2 >= 0 && bx2 <= 2 ) ax2 = bx2 ;
00078    if( bx3 >= 0 && bx3 <= 2 ) ax3 = bx3 ;
00079 
00080    dcode = dc ;
00081    return ;
00082 }
00083 
00084 /*--------------------------------------------------------------------*/
00085 
00086 static float init_dth1=0.0 , init_dth2=0.0 , init_dth3=0.0 ;
00087 static float init_dx  =0.0 , init_dy  =0.0 , init_dz  =0.0 ;
00088 
00089 #define CLEAR_INITVALS mri_3dalign_initvals(0.0,0.0,0.0,0.0,0.0,0.0)
00090 
00091 #define NONZERO_INITVALS                                        \
00092  ( init_dth1 != 0.0 || init_dth2 != 0.0 || init_dth3 != 0.0 ||  \
00093    init_dx   != 0.0 || init_dy   != 0.0 || init_dz   != 0.0   )
00094 
00095 void mri_3dalign_initvals( float th1,float th2,float th3 ,
00096                            float dx ,float dy ,float dz   )
00097 {
00098    init_dth1 = th1 ; init_dth2 = th2 ; init_dth3 = th3 ;  /* degrees */
00099    init_dx   = dx  ; init_dy   = dy  ; init_dz   = dz  ;  /* mm      */
00100 }
00101 
00102 /*--------------------------------------------------------------------*/
00103 
00104 static int regmode = MRI_QUINTIC ;
00105 static int verbose = 0 ;
00106 static int noreg   = 0 ;
00107 static int clipit  = 0 ;
00108 
00109 void mri_3dalign_method( int rmode , int verb , int norgg , int clip )
00110 {
00111    regmode = rmode ;
00112    verbose = verb ;
00113    noreg   = norgg ;
00114    clipit  = clip ;
00115    return ;
00116 }
00117 
00118 /*-------------------------------------------------------------------*/
00119 
00120 static float blurit = 0.0 ;
00121 void mri_3dalign_blurring( float bl ){ blurit = bl ; return ; }
00122 
00123 static int final_regmode = -1 ;            /* 20 Nov 1998 */
00124 void mri_3dalign_final_regmode( int frm )
00125 {
00126    final_regmode = frm ;
00127    return ;
00128 }
00129 
00130 /*-------------------------------------------------------------------*/
00131 
00132 static int xedge=-1 , yedge=-1 , zedge=-1 ;
00133 static int xfade    , yfade    , zfade    ;
00134 
00135 static int force_edging=0 ;
00136 
00137 void mri_3dalign_edging( int x , int y , int z )  /* 10 Dec 2000 */
00138 {
00139    xedge = x ; yedge = y ; zedge = z ;
00140 }
00141 
00142 void mri_3dalign_force_edging( int n )
00143 {
00144    force_edging = n ;
00145 }
00146 
00147 void mri_3dalign_edging_default( int nx , int ny , int nz )
00148 {
00149    char *ef=my_getenv("AFNI_VOLREG_EDGING") , *eq ;
00150 
00151    if( ef == NULL ){                  /* the 5% solution */
00152       xfade = (int)(0.05*nx+0.5) ;
00153       yfade = (int)(0.05*ny+0.5) ;
00154       zfade = (int)(0.05*nz+0.5) ;
00155    } else {
00156       float ff = strtod(ef,&eq) ;
00157       if( ff < 0 ){                   /* again */
00158          xfade = (int)(0.05*nx+0.5) ;
00159          yfade = (int)(0.05*ny+0.5) ;
00160          zfade = (int)(0.05*nz+0.5) ;
00161       } else {
00162          if( *eq == '%' ){            /* the whatever % solution */
00163             xfade = (int)(0.01*ff*nx+0.5) ;
00164             yfade = (int)(0.01*ff*ny+0.5) ;
00165             zfade = (int)(0.01*ff*nz+0.5) ;
00166          } else {                     /* the fixed value solution */
00167             xfade = (int)( MIN(0.25*nx,ff) ) ;
00168             yfade = (int)( MIN(0.25*ny,ff) ) ;
00169             zfade = (int)( MIN(0.25*nz,ff) ) ;
00170          }
00171       }
00172    }
00173 }
00174 
00175 /*--------------------------------------------------------------------
00176    Inputs: imbase = base image for alignment
00177            imwt   = image of weight factors to align to
00178                       (if NULL, will generate one internally)
00179 
00180    Output: pointer to a MRI_3dalign_basis struct, for later use.
00181            The malloc-ed data in there can be freed using
00182            routine MRI_3dalign_cleanup.
00183 ----------------------------------------------------------------------*/
00184 
00185 MRI_3dalign_basis * mri_3dalign_setup( MRI_IMAGE *imbase , MRI_IMAGE *imwt )
00186 {
00187    MRI_IMAGE *bim , *pim , *mim , *dim , *imww , *cim ;
00188    float *dar , *par , *mar ;
00189    float delta , dx,dy,dz ;
00190    int ii ;
00191    MRI_IMARR * fitim  =NULL;
00192    double * chol_fitim=NULL ;
00193    MRI_3dalign_basis *basis = NULL ;
00194 
00195 ENTRY("mri_3dalign_setup") ;
00196 
00197    if( !MRI_IS_3D(imbase) ){
00198      fprintf(stderr,"\n*** mri_3dalign_setup: cannot use nD images!\a\n") ;
00199      RETURN( NULL );
00200    }
00201 
00202    /*--- create output struct ---*/
00203 
00204    basis = (MRI_3dalign_basis *) malloc( sizeof(MRI_3dalign_basis) ) ;
00205 
00206    /*-- local copy of input image --*/
00207 
00208    cim = mri_to_float( imbase ) ;
00209 
00210    dx = fabs(cim->dx) ; if( dx == 0.0 ) dx = 1.0 ;
00211    dy = fabs(cim->dy) ; if( dy == 0.0 ) dy = 1.0 ;
00212    dz = fabs(cim->dz) ; if( dz == 0.0 ) dz = 1.0 ;
00213 
00214    /*--- get the weighting image ---*/
00215 
00216    if( imwt != NULL &&
00217        (imwt->nx != cim->nx || imwt->ny != cim->ny || imwt->nz != cim->nz) ){
00218 
00219       fprintf(stderr,"*** WARNING: in mri_3dalign_setup, weight image mismatch!\n") ;
00220       imwt = NULL ;
00221    }
00222 
00223    /* make weight up from the base if it isn't supplied */
00224 
00225    if( imwt == NULL ){
00226       int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny , nxyz=nxy*nz ;
00227       int ii ;
00228       float * f , clip ;
00229 
00230       /* copy base image */
00231 
00232       imww = mri_to_float( cim ) ; f = MRI_FLOAT_PTR(imww) ;
00233 
00234       if( verbose ) fprintf(stderr,"  initializing weight") ;
00235 
00236       for( ii=0 ; ii < nxyz ; ii++ ) f[ii] = fabs(f[ii]) ;  /* 16 Nov 1998 */
00237 
00238 #if 1
00239       EDIT_blur_volume_3d( nx,ny,nz , dx,dy,dz ,
00240                            MRI_float , f , 3.0*dx , 3.0*dy , 3.0*dz ) ;
00241 #else
00242       MRI_5blur_inplace_3D( imww ) ;  /* 07 Jun 2002 */
00243 #endif
00244 
00245       if( verbose ) fprintf(stderr,":") ;
00246 
00247       clip  = 0.025 * mri_max(imww) ;
00248       for( ii=0 ; ii < nxyz ; ii++ ) if( f[ii] < clip ) f[ii] = 0.0 ;
00249 
00250    } else {
00251       imww = mri_to_float( imwt ) ;  /* just copy input weight image */
00252 
00253       if( wproc ){  /* 06 Jun 2002: process input weight */
00254         int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny , nxyz=nxy*nz ;
00255         int ii ;
00256         float * f , clip ;
00257 
00258         if( verbose ) fprintf(stderr,"  processing weight") ;
00259         f = MRI_FLOAT_PTR(imww) ;
00260         for( ii=0 ; ii < nxyz ; ii++ ) f[ii] = fabs(f[ii]) ;  /* 16 Nov 1998 */
00261 #if 1
00262         EDIT_blur_volume_3d( nx,ny,nz , dx,dy,dz ,
00263                              MRI_float , f , 3.0*dx , 3.0*dy , 3.0*dz ) ;
00264 #else
00265         MRI_5blur_inplace_3D( imww ) ;  /* 07 Jun 2002 */
00266 #endif
00267         if( verbose ) fprintf(stderr,":") ;
00268         clip  = 0.025 * mri_max(imww) ;
00269         for( ii=0 ; ii < nxyz ; ii++ ) if( f[ii] < clip ) f[ii] = 0.0 ;
00270       }
00271    }
00272 
00273    /*-- 10 Dec 2000: user-controlled fade out around the edges --*/
00274 
00275    if( imwt == NULL || force_edging ){
00276      int ff , ii,jj,kk ;
00277      int nx=cim->nx , ny=cim->ny , nz=cim->nz , nxy = nx*ny ;
00278      float *f = MRI_FLOAT_PTR(imww) ;
00279 
00280      xfade = xedge ; yfade = yedge ; zfade = zedge ;  /* static variables */
00281 
00282      if( xfade < 0 || yfade < 0 || zfade < 0 )
00283         mri_3dalign_edging_default(nx,ny,nz) ;        /* reassign fades */
00284 
00285 #define FF(i,j,k) f[(i)+(j)*nx+(k)*nxy]
00286 
00287       for( jj=0 ; jj < ny ; jj++ )
00288          for( ii=0 ; ii < nx ; ii++ )
00289             for( ff=0 ; ff < zfade ; ff++ )
00290                FF(ii,jj,ff) = FF(ii,jj,nz-1-ff) = 0.0 ;
00291 
00292       for( kk=0 ; kk < nz ; kk++ )
00293          for( jj=0 ; jj < ny ; jj++ )
00294             for( ff=0 ; ff < xfade ; ff++ )
00295                FF(ff,jj,kk) = FF(nx-1-ff,jj,kk) = 0.0 ;
00296 
00297       for( kk=0 ; kk < nz ; kk++ )
00298          for( ii=0 ; ii < nx ; ii++ )
00299             for( ff=0 ; ff < yfade ; ff++ )
00300                FF(ii,ff,kk) = FF(ii,ny-1-ff,kk) = 0.0 ;
00301    }
00302 
00303    /*-- 06 Jun 2002: compute wtrimmed volume size --*/
00304 
00305    basis->xa = -1 ;  /* flag for no wtrim */
00306 
00307    if( wtrim ){
00308      int xa=-1,xb , ya,yb , za,zb ;
00309      MRI_autobbox( imww , &xa,&xb , &ya,&yb , &za,&zb ) ;
00310      if( xa >= 0 ){
00311        float nxyz = imww->nx * imww->ny * imww->nz ;
00312        float nttt = (xb-xa+1)*(yb-ya+1)*(zb-za+1) ;
00313        float trat = 100.0 * nttt / nxyz ;
00314 
00315        if( verbose )
00316          fprintf(stderr,"  wtrim: [%d..%d]x[%d..%d]x[%d..%d]"
00317                         " = %d voxels kept, out of %d (%.1f%%)\n" ,
00318                  xa,xb , ya,yb , za,zb , (int)nttt , (int)nxyz , trat ) ;
00319 
00320        /* keep trimming if saves at least 10% per volume */
00321 
00322        if( trat < 90.0 ){
00323          basis->xa = xa ; basis->xb = xb ;
00324          basis->ya = ya ; basis->yb = yb ;
00325          basis->za = za ; basis->zb = zb ;
00326 
00327          TRIM(imww) ;
00328        } else if( verbose ){
00329          fprintf(stderr,"  skipping use of trim - too little savings\n");
00330        }
00331      }
00332    }
00333 
00334    VRANG("weight",imww) ;
00335    if( verbose ){
00336      float *f=MRI_FLOAT_PTR(imww) , perc ;
00337      int ii , nxyz = imww->nvox , nzer=0 ;
00338      for( ii=0 ; ii < nxyz ; ii++ ) nzer += (f[ii] == 0.0) ;
00339      perc = (100.0*nzer)/nxyz ;
00340      fprintf(stderr,"  # zero weights=%d out of %d (%.1f%%)\n",nzer,nxyz,perc);
00341    }
00342 
00343    /*-- base image --*/
00344 
00345    INIT_IMARR( fitim ) ;             /* array of fitting images */
00346 
00347    if( DOTRIM ) bim = IMTRIM(cim) ;  /* a trimmed duplicate */
00348    else         bim = cim ;          /* base image */
00349    ADDTO_IMARR( fitim , bim ) ;
00350 
00351    THD_rota_method( regmode ) ;
00352 
00353 #ifndef MEGA
00354 #define MEGA (1024*1024)
00355 #endif
00356    if( verbose ) fprintf(stderr ,
00357                          "  mri_3dalign: using %d Mbytes of workspace\n" ,
00358                          10 * bim->nvox * bim->pixel_size / MEGA ) ;
00359 
00360    /*-- d/d(th1) image [angles in degrees here] --*/
00361 
00362    if( verbose ) fprintf(stderr,"  initializing d/d(th1)") ;
00363 
00364    delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
00365    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00366 
00367    pim = THD_rota3D( cim , ax1,delta , ax2,0.0 , ax3,0.0 ,
00368                      dcode , 0.0 , 0.0 , 0.0 ) ;
00369    if( verbose ) fprintf(stderr,":") ;
00370 
00371    mim = THD_rota3D( cim , ax1,-delta , ax2,0.0 , ax3,0.0 ,
00372                      dcode , 0.0 , 0.0 , 0.0 ) ;
00373    if( verbose ) fprintf(stderr,":") ;
00374 
00375    dim   = mri_new_conforming( cim , MRI_float ) ;
00376    delta = sinit * 0.5 * DFAC / delta ;
00377    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00378    for( ii=0 ; ii < dim->nvox ; ii++ )
00379       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00380    mri_free(pim) ; mri_free(mim) ;
00381    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00382    VRANG("d/d(th1)",dim) ;
00383 
00384    /*-- d/d(th2) image --*/
00385 
00386    if( verbose ) fprintf(stderr,"  initializing d/d(th2)") ;
00387 
00388    delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
00389    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00390 
00391    pim = THD_rota3D( cim , ax1,0.0 , ax2,delta , ax3,0.0 ,
00392                      dcode , 0.0 , 0.0 , 0.0 ) ;
00393    if( verbose ) fprintf(stderr,":") ;
00394 
00395    mim = THD_rota3D( cim , ax1,0.0 , ax2,-delta , ax3,0.0 ,
00396                      dcode , 0.0 , 0.0 , 0.0 ) ;
00397    if( verbose ) fprintf(stderr,":") ;
00398 
00399    dim   = mri_new_conforming( cim , MRI_float ) ;
00400    delta = sinit * 0.5 * DFAC / delta ;
00401    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00402    for( ii=0 ; ii < dim->nvox ; ii++ )
00403       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00404    mri_free(pim) ; mri_free(mim) ;
00405    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00406    VRANG("d/d(th2)",dim) ;
00407 
00408    /*-- d/d(th3) image --*/
00409 
00410    if( verbose ) fprintf(stderr,"  initializing d/d(th3)") ;
00411 
00412    delta = 2.0*delfac/( cim->nx + cim->ny + cim->nz ) ;
00413    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00414 
00415    pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,delta ,
00416                      dcode , 0.0 , 0.0 , 0.0 ) ;
00417    if( verbose ) fprintf(stderr,":") ;
00418 
00419    mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,-delta ,
00420                      dcode , 0.0 , 0.0 , 0.0 ) ;
00421    if( verbose ) fprintf(stderr,":") ;
00422 
00423    dim   = mri_new_conforming( cim , MRI_float ) ;
00424    delta = sinit * 0.5 * DFAC / delta ;
00425    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00426    for( ii=0 ; ii < dim->nvox ; ii++ )
00427       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00428    mri_free(pim) ; mri_free(mim) ;
00429    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00430    VRANG("d/d(th3)",dim) ;
00431 
00432    /*-- d/dx image --*/
00433 
00434    if( verbose ) fprintf(stderr,"  initializing d/dx") ;
00435 
00436    delta = delfac * dx ;
00437    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00438 
00439    pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00440                      dcode , delta , 0.0 , 0.0 ) ;
00441    if( verbose ) fprintf(stderr,":") ;
00442 
00443    mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00444                      dcode , -delta , 0.0 , 0.0 ) ;
00445    if( verbose ) fprintf(stderr,":") ;
00446 
00447    dim   = mri_new_conforming( cim , MRI_float ) ;
00448    delta = sinit * 0.5 / delta ;
00449    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00450    for( ii=0 ; ii < dim->nvox ; ii++ )
00451       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00452    mri_free(pim) ; mri_free(mim) ;
00453    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00454    VRANG("d/dx",dim) ;
00455 
00456    /*-- d/dy image --*/
00457 
00458    if( verbose ) fprintf(stderr,"  initializing d/dy") ;
00459 
00460    delta = delfac * dy ;
00461    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00462 
00463    pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00464                      dcode , 0.0 , delta , 0.0 ) ;
00465    if( verbose ) fprintf(stderr,":") ;
00466 
00467    mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00468                      dcode , 0.0 , -delta , 0.0 ) ;
00469    if( verbose ) fprintf(stderr,":") ;
00470 
00471    dim   = mri_new_conforming( cim , MRI_float ) ;
00472    delta = sinit * 0.5 / delta ;
00473    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00474    for( ii=0 ; ii < dim->nvox ; ii++ )
00475       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00476    mri_free(pim) ; mri_free(mim) ;
00477    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00478    VRANG("d/dy",dim) ;
00479 
00480    /*-- d/dz image --*/
00481 
00482    if( verbose ) fprintf(stderr,"  initializing d/dz") ;
00483 
00484    delta = delfac * dz ;
00485    if( verbose ) fprintf(stderr,"[delta=%g]",delta) ;
00486 
00487    pim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00488                      dcode , 0.0 , 0.0 , delta ) ;
00489    if( verbose ) fprintf(stderr,":") ;
00490 
00491    mim = THD_rota3D( cim , ax1,0.0 , ax2,0.0 , ax3,0.0 ,
00492                      dcode , 0.0 , 0.0 , -delta ) ;
00493    if( verbose ) fprintf(stderr,":") ;
00494 
00495    dim   = mri_new_conforming( cim , MRI_float ) ;
00496    delta = sinit * 0.5 / delta ;
00497    dar   = MRI_FLOAT_PTR(dim) ; par = MRI_FLOAT_PTR(pim) ; mar = MRI_FLOAT_PTR(mim) ;
00498    for( ii=0 ; ii < dim->nvox ; ii++ )
00499       dar[ii] = delta * ( mar[ii] - par[ii] ) ;
00500    mri_free(pim) ; mri_free(mim) ;
00501    TRIM( dim ) ; ADDTO_IMARR( fitim , dim ) ;
00502    VRANG("d/dz",dim) ;
00503 
00504    /*-- done with input copy, unless it is same as base for lsqfit --*/
00505 
00506    if( cim != bim ) mri_free(cim) ;
00507 
00508    /*-- initialize linear least squares --*/
00509 
00510    if( verbose ) fprintf(stderr,"  initializing least squares\n") ;
00511 
00512    chol_fitim = mri_startup_lsqfit( fitim , imww ) ;
00513 
00514    mri_free(imww) ;
00515 
00516    /*-- save stuff --*/
00517 
00518    basis->fitim      = fitim ;
00519    basis->chol_fitim = chol_fitim ;
00520 
00521    RETURN( basis );
00522 }
00523 
00524 /*-----------------------------------------------------------------------
00525    Input:   basis  = MRI_3dalign_basis * return from setup routine above.
00526             im     = MRI_IMAGE * to align to base image
00527 
00528    Output:  Return value is aligned image;
00529             *dx, *dy, *dz, *th1, *th2, *th3 are set to estimated
00530             alignment parameters.  Note that returned image is floats.
00531 -------------------------------------------------------------------------*/
00532 
00533 MRI_IMAGE * mri_3dalign_one( MRI_3dalign_basis * basis , MRI_IMAGE * im ,
00534                              float *th1 , float *th2 , float *th3 ,
00535                              float *dx  , float *dy  , float *dz    )
00536 {
00537    MRI_IMARR * fitim ;
00538    double * chol_fitim=NULL ;
00539    float * fit , *dfit ;
00540    int iter , good , ii ;
00541    float dxt , dyt , dzt , ftop,fbot ;
00542    MRI_IMAGE * tim , * fim ;
00543 
00544 ENTRY("mri_3dalign_one") ;
00545 
00546    fitim      = basis->fitim ;
00547    chol_fitim = basis->chol_fitim ;
00548 
00549    /* use original image if possible */
00550 
00551    if( im->kind == MRI_float ) fim = im ;
00552    else                        fim = mri_to_float( im ) ;
00553 
00554    iter = 0 ;
00555 
00556    THD_rota_method( regmode ) ;
00557 
00558    /* convert displacement threshold from voxels to mm in each direction */
00559 
00560    dxt = (im->dx != 0.0) ? (fabs(im->dx) * dxy_thresh) : dxy_thresh ;
00561    dyt = (im->dy != 0.0) ? (fabs(im->dy) * dxy_thresh) : dxy_thresh ;
00562    dzt = (im->dz != 0.0) ? (fabs(im->dz) * dxy_thresh) : dxy_thresh ;
00563 
00564    if( NONZERO_INITVALS ){                                    /* 04 Sep 2000 */
00565       fit = (float *) malloc(sizeof(float)*7) ;
00566       fit[0] = 1.0 ;
00567       fit[1] = init_dth1; fit[2] = init_dth2; fit[3] = init_dth3; /* degrees */
00568       fit[4] = init_dx  ; fit[5] = init_dy  ; fit[6] = init_dz  ; /* mm      */
00569 
00570       good = 1 ;
00571    } else {
00572 
00573       /* 06 Jun 2002: do initial fit with trimmed image, if ordered */
00574 
00575       if( DOTRIM ){
00576         tim = IMTRIM(fim) ;
00577         fit = mri_delayed_lsqfit( tim , fitim , chol_fitim ) ;
00578         mri_free( tim ) ;
00579       } else {                                /* L2 fit input image */
00580         fit = mri_delayed_lsqfit( fim , fitim , chol_fitim ) ;
00581       }
00582 
00583       good = ( 10.0*fabs(fit[4]) > dxt        || 10.0*fabs(fit[5]) > dyt        ||
00584                10.0*fabs(fit[6]) > dzt        || 10.0*fabs(fit[1]) > phi_thresh ||
00585                10.0*fabs(fit[2]) > phi_thresh || 10.0*fabs(fit[3]) > phi_thresh   ) ;
00586    }
00587 
00588    if( verbose )
00589       fprintf(stderr,
00590              "\nFirst fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
00591              fit[0] , fit[1] , fit[2] , fit[3] , fit[4] , fit[5] , fit[6] ) ;
00592 
00593    /*-- iterate fit --*/
00594 
00595    while( good ){
00596       tim = THD_rota3D( fim ,
00597                         ax1,fit[1]*DFAC , ax2,fit[2]*DFAC , ax3,fit[3]*DFAC ,
00598                         dcode , fit[4],fit[5],fit[6] ) ;
00599 
00600       TRIM(tim) ; /* 06 Jun 2002: trim it if ordered to */
00601 
00602       dfit = mri_delayed_lsqfit( tim , fitim , chol_fitim ) ; /* delta angle/shift */
00603       mri_free( tim ) ;
00604 
00605       fit[1] += dfit[1] ; fit[2] += dfit[2] ; fit[3] += dfit[3] ;  /* accumulate  */
00606       fit[4] += dfit[4] ; fit[5] += dfit[5] ; fit[6] += dfit[6] ;  /* angle/shift */
00607 
00608       if( verbose ){
00609          fprintf(stderr,
00610                  "Delta fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
00611                  dfit[0], dfit[1], dfit[2], dfit[3], dfit[4], dfit[5], dfit[6] ) ;
00612          fprintf(stderr,
00613                  "Total fit: %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
00614                  dfit[0], fit[1], fit[2], fit[3], fit[4], fit[5], fit[6] ) ;
00615       }
00616 
00617       good = (++iter < max_iter) &&
00618              ( fabs(dfit[4]) > dxt        || fabs(dfit[5]) > dyt        ||
00619                fabs(dfit[6]) > dzt        || fabs(dfit[1]) > phi_thresh ||
00620                fabs(dfit[2]) > phi_thresh || fabs(dfit[3]) > phi_thresh   ) ;
00621 
00622       free(dfit) ; dfit = NULL ;
00623    } /* end while */
00624 
00625    if( verbose ) fprintf(stderr,"Iteration complete at %d steps\n",iter) ;
00626 
00627    /*-- save final alignment parameters --*/
00628 
00629    if( th1 != NULL ) *th1 = fit[1]*DFAC ;  /* convert to radians */
00630    if( th2 != NULL ) *th2 = fit[2]*DFAC ;
00631    if( th3 != NULL ) *th3 = fit[3]*DFAC ;
00632    if( dx  != NULL ) *dx  = fit[4] ;
00633    if( dy  != NULL ) *dy  = fit[5] ;
00634    if( dz  != NULL ) *dz  = fit[6] ;
00635 
00636    /*-- do the actual realignment --*/
00637 
00638    if( ! noreg ){
00639       if( final_regmode < 0 ) final_regmode = regmode ;  /* 20 Nov 1998 */
00640       THD_rota_method( final_regmode ) ;
00641       tim = THD_rota3D( fim ,
00642                         ax1,fit[1]*DFAC , ax2,fit[2]*DFAC , ax3,fit[3]*DFAC ,
00643                         dcode , fit[4],fit[5],fit[6] ) ;
00644    } else {
00645       tim = NULL ;
00646    }
00647 
00648    if( tim != NULL && clipit &&
00649        (final_regmode == MRI_QUINTIC || final_regmode==MRI_CUBIC  ||
00650         final_regmode == MRI_HEPTIC  || final_regmode==MRI_FOURIER  ) ){
00651 
00652       register int ii ;
00653       register float ftop , fbot , * tar ;
00654 
00655       ftop = mri_max( fim ); fbot = mri_min( fim );
00656       tar  = MRI_FLOAT_PTR(tim) ;
00657       for( ii=0 ; ii < tim->nvox ; ii++ ){
00658               if( tar[ii] < fbot ) tar[ii] = fbot ;
00659          else if( tar[ii] > ftop ) tar[ii] = ftop ;
00660       }
00661    }
00662 
00663    if( fim != im ) mri_free(fim) ;  /* if it was a copy, junk it */
00664 
00665    RETURN( tim );  /* 10-4, good buddy */
00666 }
00667 
00668 /*--------------------------------------------------------------------*/
00669 
00670 MRI_IMARR * mri_3dalign_many( MRI_IMAGE * im , MRI_IMAGE * imwt , MRI_IMARR * ims ,
00671                               float *th1 , float *th2 , float *th3 ,
00672                               float *dx  , float *dy  , float *dz   )
00673 {
00674    int kim ;
00675    MRI_IMAGE * tim ;
00676    MRI_IMARR * alim ;
00677    MRI_3dalign_basis * basis ;
00678 
00679 ENTRY("mri_3dalign_many") ;
00680 
00681    basis = mri_3dalign_setup( im , imwt ) ;
00682    if( basis == NULL ) RETURN( NULL );
00683 
00684    INIT_IMARR( alim ) ;
00685 
00686 #define PK(x) ( (x!=NULL) ? (x+kim) : NULL )
00687 
00688    for( kim=0 ; kim < ims->num ; kim++ ){
00689       tim = mri_3dalign_one( basis , ims->imarr[kim] ,
00690                              PK(th1), PK(th2), PK(th3),
00691                              PK(dx) , PK(dy) , PK(dz)  ) ;
00692       ADDTO_IMARR(alim,tim) ;
00693    }
00694 
00695    mri_3dalign_cleanup( basis ) ;
00696    RETURN( alim );
00697 }
00698 
00699 /*--------------------------------------------------------------------*/
00700 
00701 void mri_3dalign_cleanup( MRI_3dalign_basis * basis )
00702 {
00703 ENTRY("mri_3dalign_cleanup") ;
00704    if( basis == NULL ) EXRETURN ;
00705 
00706    if( basis->fitim      != NULL ){ DESTROY_IMARR( basis->fitim ) ; }
00707    if( basis->chol_fitim != NULL ){ free(basis->chol_fitim) ; }
00708 
00709    free(basis) ; EXRETURN ;
00710 }