00644 {
00645    THD_dmat33 tod ;
00646 
00647    LOAD_ZERO_DMAT(tod) ;
00648 
00649    switch( dset->daxes->xxorient ){
00650       case ORI_R2L_TYPE: tod.mat[0][0] =  1.0 ; break ;
00651       case ORI_L2R_TYPE: tod.mat[0][0] = -1.0 ; break ;
00652       case ORI_P2A_TYPE: tod.mat[1][0] = -1.0 ; break ;
00653       case ORI_A2P_TYPE: tod.mat[1][0] =  1.0 ; break ;
00654       case ORI_I2S_TYPE: tod.mat[2][0] =  1.0 ; break ;
00655       case ORI_S2I_TYPE: tod.mat[2][0] = -1.0 ; break ;
00656 
00657       default: THD_FATAL_ERROR("illegal xxorient code") ;
00658    }
00659 
00660    switch( dset->daxes->yyorient ){
00661       case ORI_R2L_TYPE: tod.mat[0][1] =  1.0 ; break ;
00662       case ORI_L2R_TYPE: tod.mat[0][1] = -1.0 ; break ;
00663       case ORI_P2A_TYPE: tod.mat[1][1] = -1.0 ; break ;
00664       case ORI_A2P_TYPE: tod.mat[1][1] =  1.0 ; break ;
00665       case ORI_I2S_TYPE: tod.mat[2][1] =  1.0 ; break ;
00666       case ORI_S2I_TYPE: tod.mat[2][1] = -1.0 ; break ;
00667 
00668       default: THD_FATAL_ERROR("illegal yyorient code") ;
00669    }
00670 
00671    switch( dset->daxes->zzorient ){
00672       case ORI_R2L_TYPE: tod.mat[0][2] =  1.0 ; break ;
00673       case ORI_L2R_TYPE: tod.mat[0][2] = -1.0 ; break ;
00674       case ORI_P2A_TYPE: tod.mat[1][2] = -1.0 ; break ;
00675       case ORI_A2P_TYPE: tod.mat[1][2] =  1.0 ; break ;
00676       case ORI_I2S_TYPE: tod.mat[2][2] =  1.0 ; break ;
00677       case ORI_S2I_TYPE: tod.mat[2][2] = -1.0 ; break ;
00678 
00679       default: THD_FATAL_ERROR("illegal zxorient code") ;
00680    }
00681 
00682    return tod ;
00683 }

00023 {
00024    THD_dfvec3 *xx , *yy , dv ;
00025    int nvec , ii,jj, iarg ;
00026    THD_dvecmat rt , rtinv ;
00027    THD_dmat33  pp,ppt , rr ;
00028    THD_dfvec3  tt ;
00029 
00030    THD_3dim_dataset *mset=NULL , *dset=NULL ;
00031    double *ww=NULL ;
00032    int     nww=0 ;
00033    int keeptags=1 , wtval=0 , verb=0 , dummy=0 ;
00034    char * prefix = "tagalign" , *mfile=NULL ;
00035 
00036    float *fvol , cbot,ctop , dsum ;
00037    int nval , nvox , clipit , ival ;
00038 
00039    float matar[12] ;
00040 
00041    int use_3dWarp=1 , matrix_type=ROTATION ;
00042 
00043    /*--- help? ---*/
00044 
00045    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00046       printf("Usage: 3dTagalign [options] dset\n"
00047              "Rotates/translates dataset 'dset' to be aligned with the master,\n"
00048              "using the tagsets embedded in their .HEAD files.\n"
00049              "\n"
00050              "Options:\n"
00051              " -master mset  = Use dataset 'mset' as the master dataset\n"
00052              "                   [this is a nonoptional option]\n"
00053 #if 0
00054              "\n"
00055              " -wtval        = Use the numerical value attached to each tag\n"
00056              "                   in the master as weighting factor for that tag\n"
00057              " -wt1D ff      = Read the *.1D file to use as a vector of weights\n"
00058              "                   for each tag\n"
00059              "           N.B.: The default is to weight each tag the same.\n"
00060              "                   If you use weights, a larger weight means to\n"
00061              "                   count aligning that tag as more important.\n"
00062 #endif
00063              "\n"
00064              " -nokeeptags   = Don't put transformed locations of dset's tags\n"
00065              "                   into the output dataset [default = keep tags]\n"
00066              "\n"
00067              " -matvec mfile = Write the matrix+vector of the transformation to\n"
00068              "                   file 'mfile'.  This can be used as input to the\n"
00069              "                   '-matvec_out2in' option of 3dWarp, if you want\n"
00070              "                   to align other datasets in the same way (e.g.,\n"
00071              "                   functional datasets).\n"
00072 #if 0
00073              "           N.B.: The matrix+vector of the transformation is also\n"
00074              "                   saved in the .HEAD file of the output dataset\n"
00075              "                   (unless -dummy is used).  You can use the\n"
00076              "                   output dataset from 3dTagalign as the dataset\n"
00077              "                   for 3drotates's '-matvec_dset' option; this\n"
00078              "                   lets you avoid using an intermediate mfile.\n"
00079              "\n"
00080              " -3dWarp       = Use the '3dWarp' function to transform the dataset,\n"
00081              "                   rather than the '3drotate' function.  The output\n"
00082              "                   dataset will be computed on the same grid as the\n"
00083              "                   master dataset.\n"
00084              "           N.B.: This option is implied by '-affine' and '-rotscl'.\n"
00085 #endif
00086              "\n"
00087              " -rotate       = Compute the best transformation as a rotation + shift.\n"
00088              "                   This is the default.\n"
00089              "\n"
00090              " -affine       = Compute the best transformation as a general affine\n"
00091              "                   map rather than just a rotation + shift.  In all\n"
00092              "                   cases, the transformation from input to output\n"
00093              "                   coordinates is of the form\n"
00094              "                      [out] = [R] [in] + [V]\n"
00095              "                   where [R] is a 3x3 matrix and [V] is a 3-vector.\n"
00096              "                   By default, [R] is computed as a proper (det=1)\n"
00097              "                   rotation matrix (3 parameters).  The '-affine'\n"
00098              "                   option says to fit [R] as a general matrix\n"
00099              "                   (9 parameters).\n"
00100              "           N.B.: An affine transformation can rotate, rescale, and\n"
00101              "                   shear the volume.  Be sure to look at the dataset\n"
00102              "                   before and after to make sure things are OK.\n"
00103              "\n"
00104              " -rotscl       = Compute transformation as a rotation times an isotropic\n"
00105              "                   scaling; that is, [R] is an orthogonal matrix times\n"
00106              "                   a scalar.\n"
00107              "           N.B.: '-affine' and '-rotscl' do unweighted least squares.\n"
00108              "\n"
00109              " -prefix pp    = Use 'pp' as the prefix for the output dataset.\n"
00110              "                   [default = 'tagalign']\n"
00111              " -verb         = Print progress reports\n"
00112              " -dummy        = Don't actually rotate the dataset, just compute\n"
00113              "                   the transformation matrix and vector.  If\n"
00114              "                   '-matvec' is used, the mfile will be written.\n"
00115              "\n"
00116              "Nota Bene:\n"
00117              "* Cubic interpolation is used.  The transformation is carried out\n"
00118              "  using the same methods as program 3dWarp.\n"
00119              "\n"
00120              "Author: RWCox - 16 Jul 2000, etc.\n"
00121             ) ;
00122       exit(0) ;
00123    }
00124 
00125    /*- scan args -*/
00126 
00127    iarg = 1 ;
00128    while( iarg < argc && argv[iarg][0] == '-' ){
00129 
00130       /*-----*/
00131 
00132       if( strcmp(argv[iarg],"-rotate") == 0 ){   /* 22 Apr 2003 */
00133         matrix_type = ROTATION ; use_3dWarp = 1 ;
00134         iarg++ ; continue ;
00135       }
00136 
00137       /*-----*/
00138 
00139       if( strcmp(argv[iarg],"-affine") == 0 ){   /* 21 Apr 2003 */
00140         matrix_type = AFFINE ; use_3dWarp = 1 ;
00141         iarg++ ; continue ;
00142       }
00143 
00144       /*-----*/
00145 
00146       if( strcmp(argv[iarg],"-rotscl") == 0 ){   /* 22 Apr 2003 */
00147         matrix_type = ROTSCL ; use_3dWarp = 1 ;
00148         iarg++ ; continue ;
00149       }
00150 
00151 #if 0
00152       /*-----*/
00153 
00154       if( strcmp(argv[iarg],"-3dWarp") == 0 ){   /* 21 Apr 2003 */
00155         use_3dWarp = 1 ;
00156         iarg++ ; continue ;
00157       }
00158 #endif
00159 
00160       /*-----*/
00161 
00162       if( strcmp(argv[iarg],"-master") == 0 ){
00163          if( mset != NULL )                    ERREX("Can only have one -master option") ;
00164          if( ++iarg >= argc )                  ERREX("Need an argument after -master") ;
00165 
00166          mset = THD_open_dataset( argv[iarg] ) ;
00167 
00168          if( mset == NULL )                    ERREX("Can't open -master dataset") ;
00169          if( mset->tagset == NULL )            ERREX("No tags in -master dataset") ;
00170          if( TAGLIST_COUNT(mset->tagset) < 3 ) ERREX("Not enough tags in -master dataset") ;
00171 
00172          for( nvec=ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ )
00173             if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ) nvec++ ;
00174 
00175          if( nvec < 3 )                        ERREX("Not enough tags set in -master dataset") ;
00176 
00177          if( nvec < TAGLIST_COUNT(mset->tagset) )
00178             fprintf(stderr,"++ WARNING: not all tags are set in -master dataset\n") ;
00179 
00180          if( verb ) fprintf(stderr,"++ Found %d tags in -master dataset\n",nvec) ;
00181 
00182          iarg++ ; continue ;
00183       }
00184 
00185 #if 0
00186       /*-----*/
00187 
00188       if( strcmp(argv[iarg],"-wtval") == 0 ){
00189          if( ww != NULL )                      ERREX("Can't have -wtval after -wt1D") ;
00190          wtval++ ;
00191          iarg++ ; continue ;
00192       }
00193 
00194       /*-----*/
00195 
00196       if( strcmp(argv[iarg],"-wt1D") == 0 ){
00197          MRI_IMAGE * wtim ; float * wtar ;
00198 
00199          if( wtval )                           ERREX("Can't have -wt1D after -wtval") ;
00200          if( ww != NULL )                      ERREX("Can't have two -wt1D options!") ;
00201          if( ++iarg >= argc )                  ERREX("Need an argument after -wt1D") ;
00202 
00203          wtim = mri_read_1D( argv[iarg] ) ;
00204 
00205          if( wtim == NULL )                    ERREX("Can't read -wtim file") ;
00206          if( wtim->ny > 1 )                    ERREX("-wtim file has more than one columm") ;
00207 
00208          wtar = MRI_FLOAT_PTR(wtim) ;
00209          ww   = (double *) malloc(sizeof(double)*wtim->nx) ; nww = wtim->nx ;
00210          for( ii=0 ; ii < nww ; ii++ ){
00211             ww[ii] = (double) wtar[ii] ;
00212             if( ww[ii] < 0.0 )                 ERREX("Negative value found in -wt1D file") ;
00213          }
00214 
00215          mri_free(wtim) ;
00216          iarg++ ; continue ;
00217       }
00218 #endif
00219 
00220       /*-----*/
00221 
00222       if( strcmp(argv[iarg],"-nokeeptags") == 0 ){
00223          keeptags = 0 ;
00224          iarg++ ; continue ;
00225       }
00226 
00227       /*-----*/
00228 
00229       if( strncmp(argv[iarg],"-verb",5) == 0 ){
00230          verb++ ;
00231          iarg++ ; continue ;
00232       }
00233 
00234       /*-----*/
00235 
00236       if( strcmp(argv[iarg],"-dummy") == 0 ){
00237          dummy++ ;
00238          iarg++ ; continue ;
00239       }
00240 
00241       /*-----*/
00242 
00243       if( strcmp(argv[iarg],"-prefix") == 0 ){
00244          if( ++iarg >= argc )                  ERREX("Need an argument after -prefix") ;
00245          prefix = argv[iarg] ;
00246          if( !THD_filename_ok(prefix) )        ERREX("-prefix string is illegal") ;
00247          iarg++ ; continue ;
00248       }
00249 
00250       /*-----*/
00251 
00252       if( strcmp(argv[iarg],"-matvec") == 0 ){
00253          if( ++iarg >= argc )                  ERREX("Need an argument after -matvec") ;
00254          mfile = argv[iarg] ;
00255          if( !THD_filename_ok(mfile) )         ERREX("-matvec string is illegal") ;
00256          iarg++ ; continue ;
00257       }
00258 
00259 
00260       /*-----*/
00261 
00262       fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00263 
00264    } /* end of scanning command line for options */
00265 
00266    if( mset == NULL )                    ERREX("No -master option found on command line") ;
00267 
00268 #if 0
00269    if( ww != NULL && nww < nvec )        ERREX("Not enough weights found in -wt1D file") ;
00270 
00271    /*-- if -wtval, setup weights from master tag values --*/
00272 
00273    if( wtval ){
00274       ww = (double *) malloc(sizeof(double)*nvec) ; nww = nvec ;
00275       for( ii=jj=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
00276          if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
00277             ww[jj] = (double) TAG_VAL(TAGLIST_SUBTAG(mset->tagset,ii)) ;
00278 
00279             if( ww[jj] < 0.0 )           ERREX("Negative value found in -master tag values") ;
00280             jj++ ;
00281          }
00282       }
00283    }
00284 #endif
00285 
00286    /*-- read input dataset (to match to master dataset) --*/
00287 
00288    if( iarg >= argc )                    ERREX("No input dataset?") ;
00289 
00290    dset = THD_open_dataset( argv[iarg] ) ;
00291 
00292    if( dset == NULL )                    ERREX("Can't open input dataset") ;
00293    if( dset->tagset == NULL )            ERREX("No tags in input dataset") ;
00294    if( TAGLIST_COUNT(dset->tagset) !=
00295        TAGLIST_COUNT(mset->tagset)   )   ERREX("Tag counts don't match in -master and input") ;
00296 
00297    /* check if set tags match exactly */
00298 
00299    for( ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
00300       if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) !=
00301           TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii))    )
00302                                          ERREX("Set tags don't match in -master and input") ;
00303    }
00304 
00305    /*-- load vector lists: xx=master, yy=input --*/
00306 
00307    xx = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
00308    yy = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
00309    dsum = 0.0 ;
00310    for( ii=jj=0 ; ii < nvec ; ii++ ){
00311       if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
00312 
00313          LOAD_DFVEC3( xx[jj] ,                                      /* N.B.:     */
00314                      TAG_X( TAGLIST_SUBTAG(mset->tagset,ii) ) ,     /* these are */
00315                      TAG_Y( TAGLIST_SUBTAG(mset->tagset,ii) ) ,     /* in Dicom  */
00316                      TAG_Z( TAGLIST_SUBTAG(mset->tagset,ii) )  ) ;  /* order now */
00317 
00318          LOAD_DFVEC3( yy[jj] ,
00319                      TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00320                      TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00321                      TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) )  ) ;
00322 
00323          dv    = SUB_DFVEC3( xx[jj] , yy[jj] ) ;
00324          dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1] + dv.xyz[2]*dv.xyz[2] ;
00325 
00326          jj++ ;
00327       }
00328    }
00329 
00330    dsum = sqrt(dsum/nvec) ;
00331    fprintf(stderr,"++ RMS distance between tags before = %.2f mm\n" , dsum ) ;
00332 
00333    /*-- compute best transformation from mset to dset coords --*/
00334 
00335    switch( matrix_type ){
00336      default:
00337      case ROTATION:
00338        rt = DLSQ_rot_trans( nvec , yy , xx , ww ) ;  /* in thd_rot3d.c */
00339      break ;
00340 
00341      case AFFINE:
00342        rt = DLSQ_affine   ( nvec , yy , xx ) ;       /* 21 Apr 2003 */
00343      break ;
00344 
00345      case ROTSCL:
00346        rt = DLSQ_rotscl   ( nvec , yy , xx , (DSET_NZ(dset)==1) ? 2 : 3 ) ;
00347      break ;
00348    }
00349    rtinv = INV_DVECMAT(rt) ;
00350 
00351    /*-- check for floating point legality --*/
00352 
00353    nval = 0 ;
00354    for( ii=0 ; ii < 3 ; ii++ ){
00355       dsum = rt.vv.xyz[ii] ; nval += thd_floatscan(1,&dsum) ;
00356       for( jj=0 ; jj < 3 ; jj++ ){
00357          dsum = rt.mm.mat[ii][jj] ; nval += thd_floatscan(1,&dsum) ;
00358       }
00359    }
00360    if( nval > 0 ){
00361       fprintf(stderr,"** Floating point errors during calculation\n"
00362                      "** of transform matrix and translation vector\n" ) ;
00363       exit(1) ;
00364    }
00365 
00366    /*-- check for rotation matrix legality --*/
00367 
00368    dsum = DMAT_DET(rt.mm) ;
00369 
00370    if( dsum == 0.0 || (matrix_type == ROTATION && fabs(dsum-1.0) > 0.01) ){
00371      fprintf(stderr,"** Invalid transform matrix computed: tags dependent?\n"
00372                     "** computed [matrix] and [vector] follow:\n" ) ;
00373 
00374      for( ii=0 ; ii < 3 ; ii++ )
00375        fprintf(stderr,"  [ %10.5f %10.5f %10.5f ]   [ %10.5f ] \n",
00376                rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00377 
00378      exit(1) ;
00379    }
00380 
00381    /*-- print summary --*/
00382 
00383    if( verb ){
00384      fprintf(stderr,"++ Matrix & Vector [Dicom: x=R-L; y=A-P; z=I-S]\n") ;
00385      for( ii=0 ; ii < 3 ; ii++ )
00386        fprintf(stderr,"    %10.5f %10.5f %10.5f   %10.5f\n",
00387                rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00388    }
00389 
00390    if( matrix_type == ROTATION || matrix_type == ROTSCL ){
00391      double theta, costheta , dist , fac=1.0 ;
00392 
00393      if( matrix_type == ROTSCL ){
00394        fac = DMAT_DET(rt.mm); fac = fabs(fac);
00395        if( DSET_NZ(dset) == 1 ) fac = sqrt(fac) ;
00396        else                     fac = pow(fac,0.33333333);
00397      }
00398 
00399      costheta = 0.5 * sqrt(1.0 + DMAT_TRACE(rt.mm)/fac ) ;
00400      theta    = 2.0 * acos(costheta) * 180/3.14159265 ;
00401      dist     = SIZE_DFVEC3(rt.vv) ;
00402 
00403      fprintf(stderr,"++ Total rotation=%.2f degrees; translation=%.2f mm; scaling=%.2f\n",
00404              theta,dist,fac) ;
00405    }
00406 
00407    if( mfile ){
00408       FILE * mp ;
00409 
00410       if( THD_is_file(mfile) )
00411          fprintf(stderr,"++ Warning: -matvec will overwrite file %s\n",mfile) ;
00412 
00413       mp = fopen(mfile,"w") ;
00414       if( mp == NULL ){
00415          fprintf(stderr,"** Can't write to -matvec %s\n",mfile) ;
00416       } else {
00417         for( ii=0 ; ii < 3 ; ii++ )
00418           fprintf(mp,"    %10.5f %10.5f %10.5f   %10.5f\n",
00419                   rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00420         fclose(mp) ;
00421         if( verb ) fprintf(stderr,"++ Wrote matrix+vector to %s\n",mfile) ;
00422       }
00423    }
00424 
00425    if( dummy ){
00426       fprintf(stderr,"++ This was a -dummy run: no output dataset\n") ; exit(0) ;
00427    }
00428 
00429    /*-- 21 Apr 2003: transformation can be done the old way (a la 3drotate),
00430                      or the new way (a la 3dWarp).                          --*/
00431 
00432 #if 0
00433    if( !use_3dWarp ){          /**** the old way ****/
00434 
00435      /*-- now must scramble the rotation matrix and translation
00436           vector from Dicom coordinate order to dataset brick order --*/
00437 
00438      pp  = DBLE_mat_to_dicomm( dset ) ;
00439      ppt = TRANSPOSE_DMAT(pp) ;
00440      rr  = DMAT_MUL(ppt,rt.mm) ; rr = DMAT_MUL(rr,pp) ; tt = DMATVEC(ppt,rt.vv) ;
00441 
00442      /*-- now create the output dataset by screwing with the input dataset
00443           (this code is adapted from 3drotate.c)                           --*/
00444 
00445      DSET_mallocize(dset) ;
00446      DSET_load( dset ) ;
00447      if( !DSET_LOADED(dset) ) ERREX("Can't load input dataset bricks") ;
00448      dset->idcode = MCW_new_idcode() ;
00449      dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */
00450      EDIT_dset_items( dset ,
00451                          ADN_prefix , prefix ,
00452                          ADN_label1 , prefix ,
00453                       ADN_none ) ;
00454 
00455      if( THD_is_file(dset->dblk->diskptr->header_name) ){
00456         fprintf(stderr,
00457                 "** Output file %s already exists -- cannot continue!\n",
00458                 dset->dblk->diskptr->header_name ) ;
00459         exit(1) ;
00460      }
00461 
00462      tross_Make_History( "3dTagalign" , argc,argv , dset ) ;
00463 
00464      /*-- if desired, keep old tagset --*/
00465 
00466      if( keeptags ){
00467         THD_dfvec3 rv ;
00468 
00469         dsum = 0.0 ;
00470         for( jj=ii=0 ; ii < TAGLIST_COUNT(dset->tagset) ; ii++ ){
00471            if( TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii)) ){
00472               rv = DMATVEC( rt.mm , yy[jj] ) ;                     /* operating on */
00473               rv = ADD_DFVEC3( rt.vv , rv ) ;                      /* Dicom order  */
00474 
00475               dv    = SUB_DFVEC3( xx[jj] , rv ) ;
00476               dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
00477                                           + dv.xyz[2]*dv.xyz[2] ;
00478 
00479               UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00480                                   TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00481                                   TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) )  ) ;
00482 
00483               jj++ ;
00484            }
00485         }
00486         dsum = sqrt(dsum/nvec) ;
00487         fprintf(stderr,"++ RMS distance between tags after  = %.2f mm\n" , dsum ) ;
00488 
00489      } else {
00490         myXtFree(dset->tagset) ;  /* send it to the dustbin */
00491      }
00492 
00493      /*-- rotate sub-bricks --*/
00494 
00495      if( verb ) fprintf(stderr,"++ computing output BRIK") ;
00496 
00497      nvox = DSET_NVOX(dset) ;
00498      nval = DSET_NVALS(dset) ;
00499      fvol = (float *) malloc( sizeof(float) * nvox ) ;
00500 
00501      THD_rota_method( MRI_HEPTIC ) ;
00502      clipit = 1 ;
00503 
00504      for( ival=0 ; ival < nval ; ival++ ){
00505 
00506         /*- get sub-brick out of dataset -*/
00507 
00508         EDIT_coerce_type( nvox ,
00509                           DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ,
00510                           MRI_float,fvol ) ;
00511 
00512         if( clipit ){
00513            register int ii ; register float bb,tt ;
00514            bb = tt = fvol[0] ;
00515            for( ii=1 ; ii < nvox ; ii++ ){
00516                    if( fvol[ii] < bb ) bb = fvol[ii] ;
00517               else if( fvol[ii] > tt ) tt = fvol[ii] ;
00518            }
00519            cbot = bb ; ctop = tt ;
00520         }
00521 
00522         if( verb && nval < 5 ) fprintf(stderr,".") ;
00523 
00524         /*- rotate it -*/
00525 
00526         THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
00527                              fabs(DSET_DX(dset)) , fabs(DSET_DY(dset)) ,
00528                                                    fabs(DSET_DZ(dset)) ,
00529                              fvol , rr , tt ) ;
00530 
00531         if( verb ) fprintf(stderr,".") ;
00532 
00533         if( clipit ){
00534            register int ii ; register float bb,tt ;
00535            bb = cbot ; tt = ctop ;
00536            for( ii=0 ; ii < nvox ; ii++ ){
00537                    if( fvol[ii] < bb ) fvol[ii] = bb ;
00538               else if( fvol[ii] > tt ) fvol[ii] = tt ;
00539            }
00540         }
00541 
00542         if( verb && nval < 5 ) fprintf(stderr,".") ;
00543 
00544         /*- put it back into dataset -*/
00545 
00546         EDIT_coerce_type( nvox, MRI_float,fvol ,
00547                                 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) );
00548 
00549      } /* end of loop over sub-brick index */
00550 
00551      if( verb ) fprintf(stderr,":") ;
00552 
00553      /* save matrix+vector into dataset, too */
00554 
00555      UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
00556                        matar[4],matar[5],matar[6],
00557                        matar[8],matar[9],matar[10] ) ;
00558      UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
00559      THD_set_atr( dset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
00560 
00561      /* write dataset to disk */
00562 
00563      dset->dblk->master_nvals = 0 ;  /* in case this was a mastered dataset */
00564      DSET_write(dset) ;
00565 
00566      if( verb ) fprintf(stderr,"\n") ;
00567 
00568    } else
00569 #endif
00570    {   /**** the new way: use 3dWarp type transformation ****/
00571 
00572      THD_3dim_dataset *oset ;
00573      THD_vecmat tran ;
00574 
00575 #if 0
00576      DFVEC3_TO_FVEC3( rt.vv , tran.vv ) ;
00577      DMAT_TO_MAT    ( rt.mm , tran.mm ) ;
00578 #else
00579      DFVEC3_TO_FVEC3( rtinv.vv , tran.vv ) ;
00580      DMAT_TO_MAT    ( rtinv.mm , tran.mm ) ;
00581 #endif
00582 
00583      mri_warp3D_method( MRI_CUBIC ) ;
00584      oset = THD_warp3D_affine( dset, tran, mset, prefix, 0, WARP3D_NEWDSET ) ;
00585      if( oset == NULL ){
00586        fprintf(stderr,"** ERROR: THD_warp3D() fails!\n"); exit(1);
00587      }
00588 
00589      tross_Copy_History( dset , oset ) ;
00590      tross_Make_History( "3dTagalign" , argc,argv , oset ) ;
00591 
00592      UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
00593                        matar[4],matar[5],matar[6],
00594                        matar[8],matar[9],matar[10] ) ;
00595      UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
00596      THD_set_atr( oset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
00597 
00598      /*-- if desired, keep old tagset --*/
00599 
00600      if( keeptags ){
00601         THD_dfvec3 rv ;
00602 
00603         oset->tagset = myXtNew(THD_usertaglist) ;
00604         *(oset->tagset) = *(dset->tagset) ;
00605 
00606         dsum = 0.0 ;
00607         for( jj=ii=0 ; ii < TAGLIST_COUNT(oset->tagset) ; ii++ ){
00608           if( TAG_SET(TAGLIST_SUBTAG(oset->tagset,ii)) ){
00609             rv = DMATVEC( rt.mm , yy[jj] ) ;
00610             rv = ADD_DFVEC3( rt.vv , rv ) ;
00611 
00612             dv    = SUB_DFVEC3( xx[jj] , rv ) ;
00613             dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
00614                                         + dv.xyz[2]*dv.xyz[2] ;
00615 
00616             UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
00617                                 TAG_Y( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
00618                                 TAG_Z( TAGLIST_SUBTAG(oset->tagset,ii) )  ) ;
00619 
00620             jj++ ;
00621           }
00622         }
00623         dsum = sqrt(dsum/nvec) ;
00624         fprintf(stderr,"++ RMS distance between tags after  = %.2f mm\n" , dsum ) ;
00625      }
00626 
00627      DSET_write(oset) ;
00628 
00629    } /* end of 3dWarp-like work */
00630 
00631    exit(0) ;
00632 }