do{ int p=abas.nparam ;                                        \
     abas.param = (MRI_warp3D_param_def *) realloc(             \
                      (void *)abas.param ,                      \
                      sizeof(MRI_warp3D_param_def)*(p+1) ) ;    \
     abas.param[p].min   = (bb) ; abas.param[p].max   = (tt) ;  \
     abas.param[p].delta = (dd) ; abas.param[p].toler = (ll) ;  \
     abas.param[p].ident = abas.param[p].val_init = (id) ;      \
     strcpy( abas.param[p].name , (nm) ) ;                      \
     abas.param[p].fixed = 0 ;                                  \
     abas.nparam = p+1 ;                                        \
 } while(0)

00037 {
00038    THD_fvec3 xxx , yyy ;
00039 
00040    LOAD_FVEC3( xxx , ii,jj,kk ) ;
00041    yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
00042    warp_for(  yyy.xyz[0] ,   yyy.xyz[1] ,   yyy.xyz[2] ,
00043             &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
00044    yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
00045    *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
00046 }

00052 {
00053    THD_fvec3 xxx , yyy ;
00054 
00055    LOAD_FVEC3( xxx , ii,jj,kk ) ;
00056    yyy = VECMAT_VEC( ijk_to_xyz , xxx ) ;
00057    warp_inv(  yyy.xyz[0] ,   yyy.xyz[1] ,   yyy.xyz[2] ,
00058             &(xxx.xyz[0]), &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
00059    yyy = VECMAT_VEC( xyz_to_ijk , xxx ) ;
00060    *pp = yyy.xyz[0] ; *qq = yyy.xyz[1] ; *rr = yyy.xyz[2] ;
00061 }

00021                                      {
00022   memcpy( parvec , pv , sizeof(float)*np ) ;
00023   if( warp_parset != NULL ) warp_parset() ;
00024 }

00360 {
00361    THD_3dim_dataset *inset=NULL, *outset=NULL, *baset=NULL, *wtset=NULL ;
00362    MRI_warp3D_align_basis abas ;
00363    char *prefix="warpdriven" ;
00364    int warpdrive_code=-1 , nerr , nx,ny,nz , nopt=1 ;
00365    float dx,dy,dz , vp ;
00366    int kim , nvals , kpar , np , nfree ;
00367    MRI_IMAGE *qim , *tim , *fim ;
00368    float clip_baset=0.0f , clip_inset=0.0f ;
00369    char *W_1Dfile=NULL ;                      /* 04 Jan 2005 */
00370    float **parsave=NULL ;
00371    int output_float=0 ;                      /* 06 Jul 2005 */
00372    char *base_idc=NULL , *wt_idc=NULL ;
00373    int ctstart = NI_clock_time() ;
00374 
00375    /*-- help? --*/
00376 
00377    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00378      printf("\n"
00379             "Usage: 3dWarpDrive [options] dataset\n"
00380             "Warp a dataset to match another one (the base).\n"
00381             "\n"
00382             "This program is a generalization of 3dvolreg.  It tries to find\n"
00383             "a spatial transformation that warps a given dataset to match an\n"
00384             "input dataset (given by the -base option).  It will be slow.\n"
00385             "\n"
00386             "--------------------------\n"
00387             "Transform Defining Options: [exactly one of these must be used]\n"
00388             "--------------------------\n"
00389             "  -shift_only         =  3 parameters (shifts)\n"
00390             "  -shift_rotate       =  6 parameters (shifts + angles)\n"
00391             "  -shift_rotate_scale =  9 parameters (shifts + angles + scale factors)\n"
00392             "  -affine_general     = 12 parameters (3 shifts + 3x3 matrix)\n"
00393             "  -bilinear_general   = 39 parameters (3 + 3x3 + 3x3x3)\n"
00394             "\n"
00395             "  N.B.: At this time, the image intensity is NOT \n"
00396             "         adjusted for the Jacobian of the transformation.\n"
00397             "  N.B.: -bilinear_general is not yet implemented.\n"
00398             "\n"
00399             "-------------\n"
00400             "Other Options:\n"
00401             "-------------\n"
00402             "  -linear   }\n"
00403             "  -cubic    } = Chooses spatial interpolation method.\n"
00404             "  -NN       } =   [default = linear; inaccurate but fast]\n"
00405             "  -quintic  }     [for accuracy, try '-cubic -final quintic']\n"
00406             "\n"
00407             "  -base bbb   = Load dataset 'bbb' as the base to which the\n"
00408             "                  input dataset will be matched.\n"
00409             "                  [This is a mandatory option]\n"
00410             "\n"
00411             "  -verb       = Print out lots of information along the way.\n"
00412             "  -prefix ppp = Sets the prefix of the output dataset.\n"
00413             "  -input ddd  = You can put the input dataset anywhere in the\n"
00414             "                  command line option list by using the '-input'\n"
00415             "                  option, instead of always putting it last.\n"
00416             "\n"
00417             "-----------------\n"
00418             "Technical Options:\n"
00419             "-----------------\n"
00420             "  -maxite    m  = Allow up to 'm' iterations for convergence.\n"
00421             "  -delta     d  = Distance, in voxel size, used to compute\n"
00422             "                   image derivatives using finite differences.\n"
00423             "                   [Default=1.0]\n"
00424             "  -weight  wset = Set the weighting applied to each voxel\n"
00425             "                   proportional to the brick specified here.\n"
00426             "                   [Default=computed by program from base]\n"
00427             "  -thresh    t  = Set the convergence parameter to be RMS 't' voxels\n"
00428             "                   movement between iterations.  [Default=0.03]\n"
00429             "  -twopass      = Do the parameter estimation in two passes,\n"
00430             "                   coarse-but-fast first, then fine-but-slow second\n"
00431             "                   (much like the same option in program 3dvolreg).\n"
00432             "                   This is useful if large-ish warping is needed to\n"
00433             "                   align the volumes.\n"
00434             "  -final 'mode' = Set the final warp to be interpolated using 'mode'\n"
00435             "                   instead of the spatial interpolation method used\n"
00436             "                   to find the warp parameters.\n"
00437             "  -parfix n v   = Fix the n'th parameter of the warp model to\n"
00438             "                   the value 'v'.  More than one -parfix option\n"
00439             "                   can be used, to fix multiple parameters.\n"
00440             "  -1Dfile ename = Write out the warping parameters to the file\n"
00441             "                   named 'ename'.  Each sub-brick of the input\n"
00442             "                   dataset gets one line in this file.  Each\n"
00443             "                   parameter in the model gets one column.\n"
00444             "  -float        = Write output dataset in float format, even if\n"
00445             "                   input dataset is short or byte.\n"
00446             "\n"
00447             "----------------------\n"
00448             "AFFINE TRANSFORMATIONS:\n"
00449             "----------------------\n"
00450             "The options below control how the affine tranformations\n"
00451             "(-shift_rotate, -shift_rotate_scale, -affine_general)\n"
00452             "are structured in terms of 3x3 matrices:\n"
00453             "\n"
00454             "  -SDU or -SUD }= Set the order of the matrix multiplication\n"
00455             "  -DSU or -DUS }= for the affine transformations:\n"
00456             "  -USD or -UDS }=   S = triangular shear (params #10-12)\n"
00457             "                    D = diagonal scaling matrix (params #7-9)\n"
00458             "                    U = rotation matrix (params #4-6)\n"
00459             "                  Default order is '-SDU', which means that\n"
00460             "                  the U matrix is applied first, then the\n"
00461             "                  D matrix, then the S matrix.\n"
00462             "\n"
00463             "  -Supper      }= Set the S matrix to be upper or lower\n"
00464             "  -Slower      }= triangular [Default=lower triangular]\n"
00465             "\n"
00466             "  -ashift OR   }= Apply the shift parameters (#1-3) after OR\n"
00467             "  -bshift      }= before the matrix transformation. [Default=after]\n"
00468             "\n"
00469             "The matrices are specified in DICOM-ordered (x=-R+L,y=-A+P,z=-I+S)\n"
00470             "coordinates as:\n"
00471             "\n"
00472             "  [U] = [Rotate_y(param#6)] [Rotate_x(param#5)] [Rotate_z(param #4)]\n"
00473             "        (angles are in degrees)\n"
00474             "\n"
00475             "  [D] = diag( param#7 , param#8 , param#9 )\n"
00476             "\n"
00477             "        [    1        0     0 ]        [ 1 param#10 param#11 ]\n"
00478             "  [S] = [ param#10    1     0 ]   OR   [ 0    1     param#12 ]\n"
00479             "        [ param#11 param#12 1 ]        [ 0    0        1     ]\n"
00480             "\n"
00481             " For example, the default (-SDU/-ashift/-Slower) has the warp\n"
00482             " specified as [x]_warped = [S] [D] [U] [x]_in + [shift].\n"
00483             " The shift vector comprises parameters #1, #2, and #3.\n"
00484             "\n"
00485             " The goal of the program is to find the warp parameters such that\n"
00486             "   I([x]_warped) = s * J([x]_in)\n"
00487             " as closely as possible in a weighted least squares sense, where\n"
00488             " 's' is a scaling factor (an extra, invisible, parameter), J(x)\n"
00489             " is the base image, I(x) is the input image, and the weight image\n"
00490             " is a blurred copy of J(x).\n"
00491             "\n"
00492             " Using '-parfix', you can specify that some of these parameters\n"
00493             " are fixed.  For example, '-shift_rotate_scale' is equivalent\n"
00494             " '-affine_general -parfix 10 0 -parfix 11 0 -parfix 12 0'.\n"
00495             " Don't attempt to use the '-parfix' option unless you understand\n"
00496             " this example!\n"
00497             "\n"
00498             "-------------------------\n"
00499             "  RWCox - November 2004\n"
00500             "-------------------------\n"
00501            ) ;
00502      exit(0) ;
00503    }
00504 
00505    /*-- startup mechanics --*/
00506 
00507    mainENTRY("3dWarpDrive main"); machdep(); AFNI_logger("3dWarpDrive",argc,argv);
00508    PRINT_VERSION("3dWarpDrive") ;
00509 
00510    /* initialize parameters of the alignment basis struct */
00511 
00512    abas.nparam     = 0 ;
00513    abas.param      = NULL ;
00514    abas.scale_init = 1.0f ;
00515    abas.delfac     = 1.0f ;
00516    abas.tolfac     = 0.03f ;
00517    abas.twoblur    = 0.0f ;
00518    abas.regmode    = MRI_LINEAR ;
00519    abas.regfinal   = -1 ;
00520    abas.verb       = 0 ;
00521    abas.max_iter   = 0 ;
00522    abas.wtproc     = 1 ;
00523    abas.imbase     = NULL ;
00524    abas.imwt       = NULL ;
00525    abas.vwfor      = ijk_warp_for ;
00526    abas.vwinv      = ijk_warp_inv ;
00527    abas.vwset      = load_parvec ;
00528    abas.vwdet      = NULL ;
00529 
00530    abas.xedge = abas.yedge = abas.zedge = -1 ;
00531    abas.imww  = abas.imap  = abas.imps  = abas.imsk = NULL ;
00532    abas.imps_blur = NULL ;
00533 
00534    nparfix = 0 ;
00535 
00536    /*-- command line options --*/
00537 
00538    while( nopt < argc && argv[nopt][0] == '-' ){
00539 
00540      /*-----*/
00541 
00542      if( strcmp(argv[nopt],"-float") == 0 ){   /* 06 Jul 2005 */
00543        output_float = 1 ; nopt++ ; continue ;
00544      }
00545 
00546      /*-----*/
00547 
00548      if( strcmp(argv[nopt],"-1Dfile") == 0 ){
00549        if( ++nopt >= argc )
00550          ERROR_exit("Need 1 parameter afer -1Dfile!\n");
00551        W_1Dfile = strdup( argv[nopt] ) ;
00552        if( !THD_filename_ok(W_1Dfile) )
00553          ERROR_exit("Name after -1Dfile has bad characters!\n") ;
00554        nopt++ ; continue ;
00555      }
00556 
00557      /*-----*/
00558 
00559      if( strcmp(argv[nopt],"-twopass") == 0 ){
00560        float bbb = AFNI_numenv("AFNI_WARPDRIVE_TWOBLUR") ;
00561        abas.twoblur = (bbb==0.0f) ? 3.0f : bbb ;
00562        nopt++ ; continue ;
00563      }
00564 
00565      /*-----*/
00566 
00567      if( strcmp(argv[nopt],"-SDU") == 0 ){
00568        matorder = MATORDER_SDU ; nopt++ ; continue ;
00569      }
00570      if( strcmp(argv[nopt],"-SUD") == 0 ){
00571        matorder = MATORDER_SUD ; nopt++ ; continue ;
00572      }
00573      if( strcmp(argv[nopt],"-DSU") == 0 ){
00574        matorder = MATORDER_DSU ; nopt++ ; continue ;
00575      }
00576      if( strcmp(argv[nopt],"-DUS") == 0 ){
00577        matorder = MATORDER_DUS ; nopt++ ; continue ;
00578      }
00579      if( strcmp(argv[nopt],"-USD") == 0 ){
00580        matorder = MATORDER_USD ; nopt++ ; continue ;
00581      }
00582      if( strcmp(argv[nopt],"-UDS") == 0 ){
00583        matorder = MATORDER_UDS ; nopt++ ; continue ;
00584      }
00585      if( strcmp(argv[nopt],"-ashift") == 0 ){
00586        dcode = DELTA_AFTER     ; nopt++ ; continue ;
00587      }
00588      if( strcmp(argv[nopt],"-bshift") == 0 ){
00589        dcode = DELTA_BEFORE    ; nopt++ ; continue ;
00590      }
00591      if( strcmp(argv[nopt],"-Slower") == 0 ){
00592        smat  = SMAT_LOWER      ; nopt++ ; continue ;
00593      }
00594      if( strcmp(argv[nopt],"-Supper") == 0 ){
00595        smat  = SMAT_UPPER      ; nopt++ ; continue ;
00596      }
00597 
00598      /*-----*/
00599 
00600      if( strcmp(argv[nopt],"-final") == 0 ){
00601        char *str ;
00602 
00603        if( ++nopt >= argc )
00604          ERROR_exit("Need 1 parameter afer -final!\n") ;
00605        str = argv[nopt] ; if( *str == '-' ) str++ ;
00606 
00607             if( strcmp(str,"cubic")   == 0 ) abas.regfinal = MRI_CUBIC ;
00608        else if( strcmp(str,"quintic") == 0 ) abas.regfinal = MRI_QUINTIC ;
00609        else if( strcmp(str,"linear") == 0  ) abas.regfinal = MRI_LINEAR ;
00610        else if( strcmp(str,"NN")      == 0 ) abas.regfinal = MRI_NN ;
00611        else
00612          ERROR_exit("Illegal mode after -final\n");
00613        nopt++ ; continue ;
00614      }
00615 
00616      /*-----*/
00617 
00618      if( strcmp(argv[nopt],"-parfix") == 0 ){
00619        int np , ip ; float vp ;
00620        if( ++nopt >= argc-1 )
00621          ERROR_exit("Need 2 parameters afer -parfix!\n");
00622        np = strtol( argv[nopt] , NULL , 10 ) ; nopt++ ;
00623        vp = strtod( argv[nopt] , NULL ) ;
00624        if( np <= 0 || np > MAXPAR )
00625          ERROR_exit("Param #%d after -parfix is illegal!\n",np) ;
00626        for( ip=0 ; ip < nparfix ; ip++ ){
00627          if( parfix[ip].np == np ){
00628            WARNING_message("Ignoring later -parfix option for param #%d\n",ip);
00629            break ;
00630          }
00631        }
00632        if( ip == nparfix ) nparfix++ ;
00633        parfix[ip].np = np ; parfix[ip].vp = vp ;
00634        nopt++ ; continue ;
00635      }
00636 
00637      /*-----*/
00638 
00639      if( strcmp(argv[nopt],"-shift_only") == 0 ){
00640        warpdrive_code = WARPDRIVE_SHIFT ; nopt++ ; continue ;
00641      }
00642 
00643      if( strcmp(argv[nopt],"-shift_rotate") == 0 ){
00644        warpdrive_code = WARPDRIVE_ROTATE ; nopt++ ; continue ;
00645      }
00646 
00647      if( strcmp(argv[nopt],"-shift_rotate_scale") == 0 ){
00648        warpdrive_code = WARPDRIVE_SCALE ; nopt++ ; continue ;
00649      }
00650 
00651      if( strcmp(argv[nopt],"-affine_general") == 0 ){
00652        warpdrive_code = WARPDRIVE_AFFINE ; nopt++ ; continue ;
00653      }
00654 
00655      if( strcmp(argv[nopt],"-bilinear_general") == 0 ){  /* not implemented */
00656 #if 0
00657        ERROR_exit("3dWarpDrive -bilinear_general NOT IMPLEMENTED!\n");
00658 #else
00659        WARNING_message("-bilinear_general transformations are experimental!") ;
00660 #endif
00661        warpdrive_code = WARPDRIVE_BILINEAR ; nopt++ ; continue ;
00662      }
00663 
00664      /*-----*/
00665 
00666      if( strcmp(argv[nopt],"-NN")     == 0 ){
00667        abas.regmode = MRI_NN      ; nopt++ ; continue ;
00668      }
00669      if( strcmp(argv[nopt],"-linear") == 0 ){
00670        abas.regmode = MRI_LINEAR  ; nopt++ ; continue ;
00671      }
00672      if( strcmp(argv[nopt],"-cubic")  == 0 ){
00673        abas.regmode = MRI_CUBIC   ; nopt++ ; continue ;
00674      }
00675      if( strcmp(argv[nopt],"-quintic") == 0 ){
00676        abas.regmode = MRI_QUINTIC ; nopt++ ; continue ;
00677      }
00678 
00679      /*-----*/
00680 
00681      if( strcmp(argv[nopt],"-prefix") == 0 ){
00682        if( ++nopt >= argc )
00683          ERROR_exit("Need an argument after -prefix!\n");
00684        if( !THD_filename_ok(argv[nopt]) )
00685          ERROR_exit("-prefix argument is invalid!\n");
00686        prefix = argv[nopt] ; nopt++ ; continue ;
00687      }
00688 
00689      /*-----*/
00690 
00691      if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00692        abas.verb++ ; nopt++ ; continue ;
00693      }
00694 
00695      /*-----*/
00696 
00697      if( strcmp(argv[nopt],"-base") == 0 ){
00698        if( ++nopt >= argc )
00699          ERROR_exit("Need an argument after -base!\n");
00700        baset = THD_open_dataset( argv[nopt] ) ;
00701        if( baset == NULL )
00702          ERROR_exit("Can't open -base dataset %s\n",argv[nopt]);
00703        if( DSET_NVALS(baset) > 1 )
00704          WARNING_message(
00705            "-base dataset %s has %d sub-bricks; will only use #0\n",
00706            argv[nopt],DSET_NVALS(baset) ) ;
00707        nopt++ ; continue ;
00708      }
00709 
00710      /*-----*/
00711 
00712      if( strcmp(argv[nopt],"-weight") == 0 ){
00713        if( ++nopt >= argc )
00714          ERROR_exit("Need an argument after -weight!\n");
00715        wtset = THD_open_dataset( argv[nopt] ) ;
00716        if( wtset == NULL )
00717          ERROR_exit("Can't open -weight dataset %s\n",argv[nopt]);
00718        if( DSET_NVALS(wtset) > 1 )
00719          WARNING_message(
00720            "-weight dataset %s has %d sub-bricks; will only use #0\n",
00721            argv[nopt],DSET_NVALS(wtset) ) ;
00722        nopt++ ; continue ;
00723      }
00724 
00725      /*-----*/
00726 
00727      if( strcmp(argv[nopt],"-input") == 0 ){
00728        if( ++nopt >= argc )
00729          ERROR_exit("Need an argument after -input!\n");
00730        inset = THD_open_dataset( argv[nopt] ) ;
00731        if( inset == NULL )
00732          ERROR_exit("Can't open -input dataset %s\n",argv[nopt]);
00733        nopt++ ; continue ;
00734      }
00735 
00736      /*-----*/
00737 
00738      if( strcmp(argv[nopt],"-maxite") == 0 ){
00739        int ival ;
00740        if( ++nopt >= argc )
00741          ERROR_exit("Need an argument after -maxite!\n");
00742        ival = strtol( argv[nopt] , NULL , 10 ) ;
00743        if( ival > 1 ) abas.max_iter = ival ;
00744        nopt++ ; continue ;
00745      }
00746 
00747      /*-----*/
00748 
00749      if( strcmp(argv[nopt],"-delta") == 0 ){
00750        float val ;
00751        if( ++nopt >= argc )
00752          ERROR_exit("Need an argument after -delta!\n");
00753        val = strtod( argv[nopt] , NULL ) ;
00754        if( val > 0.0499 && val < 49.99 ) abas.delfac = val ;
00755        nopt++ ; continue ;
00756      }
00757 
00758      /*-----*/
00759 
00760      if( strcmp(argv[nopt],"-thresh") == 0 ){
00761        float val ;
00762        if( ++nopt >= argc )
00763          ERROR_exit("Need an argument after -thresh!\n");
00764        val = strtod( argv[nopt] , NULL ) ;
00765        if( val > 0.001 && val < 3.01 ) abas.tolfac = val ;
00766        nopt++ ; continue ;
00767      }
00768 
00769      /*-----*/
00770 
00771      ERROR_exit("Unknown option %s\n",argv[nopt]) ;
00772 
00773    } /*--- end of loop over command line options ---*/
00774 
00775    /*-- 1 remaining argument should be a dataset --*/
00776 
00777    if( inset == NULL && nopt != argc-1 )
00778      ERROR_exit("Command line should have exactly 1 dataset!\n"
00779                 "**         Whereas there seems to be %d of them!\n",
00780              argc-nopt ) ;
00781 
00782    /*-- input dataset header --*/
00783 
00784    if( inset == NULL ){
00785      inset = THD_open_dataset( argv[nopt] ) ;
00786      if( !ISVALID_DSET(inset) )
00787        ERROR_exit("Can't open dataset %s\n",argv[nopt]);
00788    }
00789 
00790    nx = DSET_NX(inset) ; ny = DSET_NY(inset) ; nz = DSET_NZ(inset) ;
00791    dx = DSET_DX(inset) ; dy = DSET_DY(inset) ; dz = DSET_DZ(inset) ;
00792 
00793    if( abas.verb ) INFO_message("Checking inputs\n") ;
00794 
00795    /*-- parameterize the warp model --*/
00796 
00797    /*! Macro to add a parameter to the warp3D model.
00798         - nm = name of parameter
00799         - bb = min value allowed
00800         - tt = max value allowed
00801         - id = value for identity warp
00802         - dd = delta to use for stepsize
00803         - ll = tolerance for convergence test */
00804 
00805 #define ADDPAR(nm,bb,tt,id,dd,ll)                               \
00806  do{ int p=abas.nparam ;                                        \
00807      abas.param = (MRI_warp3D_param_def *) realloc(             \
00808                       (void *)abas.param ,                      \
00809                       sizeof(MRI_warp3D_param_def)*(p+1) ) ;    \
00810      abas.param[p].min   = (bb) ; abas.param[p].max   = (tt) ;  \
00811      abas.param[p].delta = (dd) ; abas.param[p].toler = (ll) ;  \
00812      abas.param[p].ident = abas.param[p].val_init = (id) ;      \
00813      strcpy( abas.param[p].name , (nm) ) ;                      \
00814      abas.param[p].fixed = 0 ;                                  \
00815      abas.nparam = p+1 ;                                        \
00816  } while(0)
00817 
00818    nerr = 0 ;
00819    if( warpdrive_code <= 0 ){
00820      ERROR_message("Need a transform-specifying option!\n");
00821      nerr++ ;
00822    } else if( WARPDRIVE_IS_AFFINE(warpdrive_code) ) {
00823 
00824        char *lab09, *lab10, *lab11 ;
00825 
00826        /* add all 12 parameters (may ignore some, later) */
00827 
00828        ADDPAR( "x-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00829        ADDPAR( "y-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00830        ADDPAR( "z-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00831 
00832        ADDPAR( "z-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;  /* degrees */
00833        ADDPAR( "x-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00834        ADDPAR( "y-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00835 
00836        ADDPAR( "x-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;  /* identity */
00837        ADDPAR( "y-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;  /*  == 1.0 */
00838        ADDPAR( "z-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00839 
00840        switch( smat ){
00841          default:
00842          case SMAT_LOWER:
00843            lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
00844          break ;
00845 
00846          case SMAT_UPPER:
00847            lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
00848          break ;
00849        }
00850        ADDPAR( lab09 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00851        ADDPAR( lab10 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00852        ADDPAR( lab11 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00853 
00854        /* initialize transform parameter vector */
00855 
00856        for( kpar=0 ; kpar < 12 ; kpar++ )
00857          parvec[kpar] = abas.param[kpar].ident ;
00858 
00859        /* initialize transformation function pointers */
00860 
00861        if( warpdrive_code == WARPDRIVE_SHIFT ){
00862          warp_parset = parset_shift ;
00863          warp_for    = warper_shift_for ;
00864          warp_inv    = warper_shift_inv ;
00865        } else {
00866          warp_parset = parset_affine ;
00867          warp_for    = warper_affine_for ;
00868          warp_inv    = warper_affine_inv ;
00869        }
00870 
00871        /* how many parameters to actually pay attention to */
00872 
00873        switch( warpdrive_code ){
00874          case WARPDRIVE_SHIFT:  abas.nparam =  3 ; break ;
00875          case WARPDRIVE_ROTATE: abas.nparam =  6 ; break ;
00876          case WARPDRIVE_SCALE:  abas.nparam =  9 ; break ;
00877          case WARPDRIVE_AFFINE: abas.nparam = 12 ; break ;
00878        }
00879 
00880    } else if( warpdrive_code == WARPDRIVE_BILINEAR ){
00881 
00882        char *lab09, *lab10, *lab11 , labxx[16] ;
00883        float xr,yr,zr,rr ;
00884 
00885        /* add all 39 parameters (may ignore some, later) */
00886 
00887        ADDPAR( "x-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00888        ADDPAR( "y-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00889        ADDPAR( "z-shift" , -100.0 , 100.0 , 0.0 , 0.0 , 0.0 ) ;
00890 
00891        ADDPAR( "z-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00892        ADDPAR( "x-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00893        ADDPAR( "y-angle" , -180.0 , 180.0 , 0.0 , 0.0 , 0.0 ) ;
00894 
00895        ADDPAR( "x-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00896        ADDPAR( "y-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00897        ADDPAR( "z-scale" , 0.618  , 1.618 , 1.0 , 0.0 , 0.0 ) ;
00898 
00899        switch( smat ){
00900          default:
00901          case SMAT_LOWER:
00902            lab09 = "y/x-shear" ; lab10 = "z/x-shear" ; lab11 = "z/y-shear" ;
00903          break ;
00904 
00905          case SMAT_UPPER:
00906            lab09 = "x/y-shear" ; lab10 = "x/z-shear" ; lab11 = "y/z-shear" ;
00907          break ;
00908        }
00909        ADDPAR( lab09 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00910        ADDPAR( lab10 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00911        ADDPAR( lab11 , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00912 
00913        xr = 0.5f*fabs(dx)*nx ; yr = 0.5f*fabs(dy)*ny ; zr = 0.5f*fabs(dz)*nz ;
00914        rr = MAX(xr,yr)       ; rr = MAX(rr,zr)       ; dd_fac = 1.0f / rr ;
00915        for( kpar=12 ; kpar < 39 ; kpar++ ){
00916          sprintf(labxx,"blin_%02d",kpar+1) ;
00917          ADDPAR( labxx , -0.2222 , 0.2222 , 0.0 , 0.0 , 0.0 ) ;
00918        }
00919 
00920        /* initialize transform parameter vector */
00921 
00922        for( kpar=0 ; kpar < 39 ; kpar++ )
00923          parvec[kpar] = abas.param[kpar].ident ;
00924 
00925        abas.vwdet = warper_bilinear_det ;
00926 
00927        /* initialize transformation function pointers */
00928 
00929        warp_parset = parset_bilinear ;
00930        warp_for    = warper_bilinear_for ;
00931        warp_inv    = warper_bilinear_inv ;
00932 
00933        /* how many parameters to actually pay attention to */
00934 
00935        abas.nparam = 39 ;
00936 
00937    } else {
00938      ERROR_message("Unimplemented transform model!\n") ;
00939      nerr++ ;
00940    }
00941 
00942    /* Deal with -parfix options; nfree will be number of free parameters */
00943 
00944    nfree = abas.nparam ;
00945    for( kpar=0 ; kpar < nparfix ; kpar++ ){
00946      np = parfix[kpar].np - 1 ; vp = parfix[kpar].vp ;
00947      if( np >= 0 && np < abas.nparam ){
00948        if( vp >= abas.param[np].min && vp <= abas.param[np].max ){
00949          abas.param[np].fixed     = 1  ;
00950          abas.param[np].val_fixed = vp ;
00951          nfree -- ;
00952        } else {        /* bad value */
00953          ERROR_message("-parfix for param #%d has illegal value!\n",np+1) ;
00954          nerr++ ;
00955        }
00956      } else {          /* bad index */
00957        WARNING_message("-parfix for param #%d is out of range 1..%d\n",
00958                         np+1 , abas.nparam+1 ) ;
00959      }
00960    }
00961    if( nfree <= 0 ){
00962      ERROR_message("No free parameters in transform model!\n") ;
00963      nerr++ ;
00964    }
00965 
00966    /* default number of iterations allowed */
00967 
00968    if( abas.max_iter <= 0 ) abas.max_iter = 11*nfree+5 ;
00969 
00970    /*-- other checks for good set of inputs --*/
00971 
00972    if( baset == NULL ){
00973      ERROR_message("Need to specify a base dataset!\n") ;
00974      nerr++ ;
00975    }
00976 
00977    if( nerr ) exit(1) ;  /** bad user!! **/
00978 
00979    if( abas.verb ) INFO_message("Creating empty output dataset\n") ;
00980 
00981    outset = EDIT_empty_copy( inset ) ;
00982 
00983    EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;
00984 
00985    if( output_float ){
00986      EDIT_dset_items( outset , ADN_datum_all , MRI_float , ADN_none ) ;
00987      for( kim=0 ; kim < nvals ; kim++ )
00988        EDIT_BRICK_FACTOR( outset , kim , 0.0 ) ;
00989    }
00990 
00991    if( THD_is_file( DSET_HEADNAME(outset) ) )
00992      ERROR_exit("Output file %s already exists!\n",DSET_HEADNAME(outset) ) ;
00993 
00994    tross_Copy_History( inset , outset ) ;
00995    tross_Make_History( "3dWarpDrive" , argc,argv , outset ) ;
00996 
00997    outset->daxes->xxorg = inset->daxes->xxorg ;
00998    outset->daxes->yyorg = inset->daxes->yyorg ;
00999    outset->daxes->zzorg = inset->daxes->zzorg ;
01000 
01001    /*-- more checks --*/
01002 
01003    if( DSET_NX(baset) != nx || DSET_NY(baset) != ny || DSET_NZ(baset) != nz ){
01004      ERROR_message("base and input grid dimensions don't match!\n") ;
01005      ERROR_message("base  is %d X %d X %d voxels\n",
01006                     DSET_NX(baset),DSET_NY(baset),DSET_NZ(baset) ) ;
01007      ERROR_exit   ("input is %d X %d X %d voxels\n",nx,ny,nz) ;
01008    }
01009 
01010    /* the following aren't fatal errors, but merit a warning slap in the face */
01011 
01012    if( FLDIF(DSET_DX(baset),dx) ||
01013        FLDIF(DSET_DY(baset),dy) || FLDIF(DSET_DZ(baset),dz) ){
01014      WARNING_message("base and input grid spacings don't match!\n") ;
01015      WARNING_message("base  grid = %.5f X %.5f X %.5f mm\n",
01016                      DSET_DX(baset),DSET_DY(baset),DSET_DZ(baset) ) ;
01017      WARNING_message("input grid = %.5f X %.5f X %.5f mm\n",
01018                      DSET_DX(inset),DSET_DY(inset),DSET_DZ(inset) ) ;
01019    }
01020 
01021    if( FLDIF(DSET_XORG(baset),DSET_XORG(inset)) ||
01022        FLDIF(DSET_YORG(baset),DSET_YORG(inset)) ||
01023        FLDIF(DSET_ZORG(baset),DSET_ZORG(inset))   ){
01024      WARNING_message("base and input grid offsets don't match!\n") ;
01025      WARNING_message("base  offsets = %.5f X %.5f X %.5f mm\n",
01026                      DSET_XORG(baset),DSET_YORG(baset),DSET_ZORG(baset) ) ;
01027      WARNING_message("input offsets = %.5f X %.5f X %.5f mm\n",
01028                      DSET_XORG(inset),DSET_YORG(inset),DSET_ZORG(inset) ) ;
01029    }
01030 
01031    if( baset->daxes->xxorient != inset->daxes->xxorient ||
01032        baset->daxes->yyorient != inset->daxes->yyorient ||
01033        baset->daxes->zzorient != inset->daxes->zzorient   ){
01034      WARNING_message("base and input orientations don't match!\n") ;
01035      WARNING_message("base  = %s X %s X %s\n",
01036              ORIENT_shortstr[baset->daxes->xxorient] ,
01037              ORIENT_shortstr[baset->daxes->yyorient] ,
01038              ORIENT_shortstr[baset->daxes->zzorient]  ) ;
01039      WARNING_message("input = %s X %s X %s\n",
01040              ORIENT_shortstr[inset->daxes->xxorient] ,
01041              ORIENT_shortstr[inset->daxes->yyorient] ,
01042              ORIENT_shortstr[inset->daxes->zzorient]  ) ;
01043    }
01044 
01045    /* however, this is a fatal error */
01046 
01047    if( wtset != NULL &&
01048       (DSET_NX(wtset) != nx || DSET_NY(wtset) != ny || DSET_NZ(wtset) != nz) ){
01049      ERROR_message("weight and input grid dimensions don't match!\n") ;
01050      ERROR_message("weight is %d X %d X %d voxels\n",
01051                    DSET_NX(wtset),DSET_NY(wtset),DSET_NZ(wtset)    ) ;
01052      ERROR_exit   ("input  is %d X %d X %d voxels\n",nx,ny,nz) ;
01053    }
01054 
01055    /* load datasets from disk; if can't do so, fatal errors all around */
01056 
01057    if( abas.verb ) INFO_message("Loading datasets\n") ;
01058 
01059    DSET_load(inset) ;
01060    if( !DSET_LOADED(inset) ){
01061      ERROR_exit("Can't load input dataset into memory!\n") ;
01062    } else {
01063      nvals = DSET_NVALS(inset) ;
01064      if( nvals == 1 ){
01065        clip_inset = THD_cliplevel( DSET_BRICK(inset,0) , 0.0 ) ;
01066        if( DSET_BRICK_FACTOR(inset,0) > 0.0f )
01067          clip_inset *= DSET_BRICK_FACTOR(inset,0) ;
01068      } else {
01069        qim = THD_median_brick( inset ) ;
01070        clip_inset = THD_cliplevel( qim , 0.0 ) ;
01071        mri_free(qim) ;
01072      }
01073    }
01074 
01075    DSET_load(baset) ;
01076    if( !DSET_LOADED(baset) ){
01077      ERROR_exit("Can't load base dataset into memory!\n") ;
01078    } else {
01079      clip_baset  = THD_cliplevel( DSET_BRICK(baset,0) , 0.0 ) ;
01080      abas.imbase = mri_to_float( DSET_BRICK(baset,0) ) ;
01081      base_idc = strdup(baset->idcode.str) ;
01082      DSET_delete(baset) ; baset = NULL ;
01083    }
01084 
01085    if( wtset != NULL ){
01086      DSET_load(wtset) ;
01087      if( !DSET_LOADED(wtset) ){
01088        ERROR_exit("Can't load weight dataset into memory!\n") ;
01089      } else {
01090        abas.imwt = mri_to_float( DSET_BRICK(wtset,0) ) ;
01091        wt_idc = strdup(wtset->idcode.str) ;
01092        DSET_delete(wtset) ; wtset = NULL ;
01093      }
01094    }
01095 
01096    /*-- set up (x,y,z) <-> (i,j,k) transformations ---*/
01097 
01098    { THD_vecmat ijk_to_inset_xyz , xyz_to_dicom ;
01099 
01100      LOAD_DIAG_MAT( ijk_to_inset_xyz.mm ,
01101                     inset->daxes->xxdel ,
01102                     inset->daxes->yydel , inset->daxes->zzdel );
01103 
01104      if( warpdrive_code == WARPDRIVE_BILINEAR ){
01105        /* define (x,y,z)=(0,0,0) at mid-point of dataset 3D array */
01106 
01107        LOAD_FVEC3   ( ijk_to_inset_xyz.vv ,
01108                       -0.5*(nx-1)*inset->daxes->xxdel ,
01109                       -0.5*(ny-1)*inset->daxes->yydel ,
01110                       -0.5*(nz-1)*inset->daxes->zzdel  ) ;
01111 
01112      } else {
01113        /* define (x,y,z) based strictly on dataset coords */
01114 
01115        LOAD_FVEC3   ( ijk_to_inset_xyz.vv ,
01116                       DSET_XORG(inset) , DSET_YORG(inset), DSET_ZORG(inset) ) ;
01117      }
01118 
01119      xyz_to_dicom.mm = inset->daxes->to_dicomm ;
01120      LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
01121 
01122      ijk_to_xyz = MUL_VECMAT( xyz_to_dicom , ijk_to_inset_xyz ) ;
01123      xyz_to_ijk = INV_VECMAT( ijk_to_xyz ) ;
01124 
01125 #if 0
01126      if( abas.verb ){
01127        DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
01128        DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
01129      }
01130 #endif
01131    }
01132 
01133    /*-- make the shell of the new dataset --*/
01134 
01135    if( clip_baset > 0.0f && clip_inset > 0.0f ){
01136      float fac = clip_inset / clip_baset ;
01137      if( fac > 0.005 && fac < 2000.0 ){ /* zss: have encountered OK data needing very large scale corrections...*/
01138        abas.scale_init = fac ;
01139        if( abas.verb ) INFO_message("Scale factor set to %.2f/%.2f=%.2g\n",
01140                                clip_baset , clip_inset , fac ) ;
01141        if (fac > 200) {
01142          fprintf(stderr,"++ Warning: Scale init = %g is large. Check output.\n", fac) ;
01143        }
01144      } else {
01145       fprintf(stderr,"-- Large scale difference between datasets.\n"
01146                      "   Scale init = %g\n"
01147                      "   3dWarpDrive might not converge.",fac) ;  
01148      }
01149    }
01150 
01151    /*===== do the hard work =====*/
01152 
01153    if( abas.verb ) INFO_message("Beginning alignment setup\n") ;
01154 
01155    /* 04 Jan 2005: set up to save the computed parameters */
01156 
01157    parsave = (float **)malloc( sizeof(float *) * abas.nparam ) ;
01158    for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01159      parsave[kpar] = (float *)calloc( sizeof(float) , nvals ) ;
01160 
01161    mri_warp3D_align_setup( &abas ) ;
01162 
01163    if( abas.verb ) INFO_message("Beginning alignment loop\n") ;
01164 
01165    /** loop over input sub-bricks **/
01166 
01167    for( kim=0 ; kim < nvals ; kim++ ){
01168 
01169      for( kpar=0 ; kpar < abas.nparam ; kpar++ ){  /** init params **/
01170        if( abas.param[kpar].fixed )
01171          abas.param[kpar].val_init = abas.param[kpar].val_fixed ;
01172        else
01173          abas.param[kpar].val_init = abas.param[kpar].ident ;
01174      }
01175 
01176      /** create copy of input brick into qim
01177          then warp-align it, with result into tim **/
01178 
01179      qim = mri_scale_to_float( DSET_BRICK_FACTOR(inset,kim) ,
01180                                DSET_BRICK(inset,kim)         ) ;
01181      tim = mri_warp3d_align_one( &abas , qim ) ;
01182      mri_free( qim ) ; DSET_unload_one( inset , kim ) ;
01183 
01184      if( abas.verb ){ DUMP_VECMAT( "end mv_for" , mv_for ) ; }
01185 
01186      /** save output parameters for later **/
01187 
01188      for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01189        parsave[kpar][kim] = abas.param[kpar].val_out ;  /* 04 Jan 2005 */
01190 
01191      /** convert output image from float to whatever **/
01192 
01193      switch( DSET_BRICK_TYPE(outset,kim) ){
01194 
01195          default:
01196            ERROR_message("Can't store bricks of type %s\n",
01197                          MRI_TYPE_name[DSET_BRICK_TYPE(inset,kim)] ) ;
01198            /* fall thru on purpose */
01199 
01200          case MRI_float:
01201            EDIT_substitute_brick( outset, kim, MRI_float, MRI_FLOAT_PTR(tim) );
01202            mri_fix_data_pointer( NULL , tim ) ; mri_free( tim ) ;
01203          break ;
01204 
01205          case MRI_short:{
01206            double fac = DSET_BRICK_FACTOR(inset,kim) ;
01207            fac = (fac >  0.0) ? 1.0/fac : 1.0 ;
01208            fim = mri_to_short(fac,tim) ; mri_free( tim ) ;
01209            EDIT_substitute_brick( outset, kim, MRI_short, MRI_SHORT_PTR(fim) );
01210            fac = (fac != 1.0) ? 1.0/fac : 0.0 ;
01211            EDIT_BRICK_FACTOR( outset , kim , fac ) ;
01212            mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
01213          }
01214          break ;
01215 
01216          case MRI_byte:
01217            vp = mri_min(tim) ;
01218            if( vp < 0.0f ){
01219              WARNING_message(
01220               "output sub-brick #%d is byte, but has negative values\n",
01221               kim ) ;
01222            }
01223            fim = mri_to_byte(tim) ; mri_free( tim ) ;
01224            EDIT_substitute_brick( outset, kim, MRI_byte, MRI_BYTE_PTR(fim) ) ;
01225            mri_fix_data_pointer( NULL , fim ) ; mri_free( fim ) ;
01226          break ;
01227       }
01228    }
01229    DSET_unload( inset ) ;
01230 
01231    /*===== hard work is done =====*/
01232 
01233    mri_warp3D_align_cleanup( &abas ) ;
01234 
01235    /*-- 06 Jul 2005:
01236         write the affine transform matrices into the output header --*/
01237 
01238    THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_IDCODE" ,
01239                                       inset->idcode.str         ) ;
01240    THD_set_string_atr( outset->dblk , "WARPDRIVE_INPUT_NAME" ,
01241                                       DSET_HEADNAME(inset)      ) ;
01242    if( base_idc != NULL )
01243      THD_set_string_atr( outset->dblk , "WARPDRIVE_BASE_IDCODE" , base_idc ) ;
01244    if( wt_idc != NULL )
01245      THD_set_string_atr( outset->dblk , "WARPDRIVE_WEIGHT_IDCODE" , wt_idc ) ;
01246 
01247    if( WARPDRIVE_IS_AFFINE(warpdrive_code) ){  /* can't do bilinear here */
01248      float matar[12] ; char anam[64] ;
01249 
01250      for( kpar=0 ; kpar < 12 ; kpar++ ) parvec[kpar] = 0.0 ;
01251 
01252      for( kim=0 ; kim < nvals ; kim++ ){
01253        for( kpar=0 ; kpar < abas.nparam ; kpar++ )  /* load params */
01254          parvec[kpar] = parsave[kpar][kim] ;
01255        parset_affine() ;                            /* compute matrices */
01256 
01257        UNLOAD_MAT(mv_for.mm,matar[0],matar[1],matar[2],
01258                             matar[4],matar[5],matar[6],
01259                             matar[8],matar[9],matar[10] ) ;
01260        UNLOAD_FVEC3(mv_for.vv,matar[3],matar[7],matar[11]) ;
01261        sprintf(anam,"WARPDRIVE_MATVEC_FOR_%06d",kim) ;
01262        THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
01263 
01264        UNLOAD_MAT(mv_inv.mm,matar[0],matar[1],matar[2],
01265                             matar[4],matar[5],matar[6],
01266                             matar[8],matar[9],matar[10] ) ;
01267        UNLOAD_FVEC3(mv_inv.vv,matar[3],matar[7],matar[11]) ;
01268        sprintf(anam,"WARPDRIVE_MATVEC_INV_%06d",kim) ;
01269        THD_set_float_atr( outset->dblk , anam , 12 , matar ) ;
01270      }
01271    }
01272 
01273    /*-- write the results to disk for all of history to see --*/
01274 
01275    DSET_write( outset ) ;  DSET_unload( outset ) ;
01276    if( abas.verb ) WROTE_DSET(outset) ;
01277 
01278    if( W_1Dfile != NULL ){
01279      FILE *fp ;
01280      if( abas.verb ) INFO_message("Writing 1Dfile: %s\n",W_1Dfile) ;
01281      if( THD_is_file(W_1Dfile) )
01282        WARNING_message("Overwriting file %s\n",W_1Dfile) ;
01283 
01284      fp = fopen( W_1Dfile , "w" ) ;
01285      if( fp != NULL ){
01286 
01287        fprintf(fp,"#") ;
01288        for( kim=0 ; kim < argc ; kim++ ) fprintf(fp," %s",argv[kim]) ;
01289        fprintf(fp,"\n") ;
01290 
01291        fprintf(fp,"#") ;
01292        for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01293          fprintf(fp," %-13.13s",abas.param[kpar].name) ;
01294        fprintf(fp,"\n") ;
01295 
01296        fprintf(fp,"#") ;
01297        for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01298          fprintf(fp," -------------") ;
01299        fprintf(fp,"\n") ;
01300 
01301        for( kim=0 ; kim < nvals ; kim++ ){
01302          for( kpar=0 ; kpar < abas.nparam ; kpar++ )
01303            fprintf(fp," %13.6g",parsave[kpar][kim]) ;
01304          fprintf(fp,"\n") ;
01305        }
01306        fclose(fp) ;
01307      }
01308 
01309    }
01310 
01311    if( abas.verb ){
01312      double tt = (NI_clock_time()-ctstart)*0.001 ;
01313      INFO_message("Total elapsed time = %.2f s\n",tt) ;
01314    }
01315 
01316    exit(0) ;
01317 }

00152 {
00153    THD_mat33 ss,dd,uu,aa,bb ;
00154    THD_fvec3 vv ;
00155 
00156 #if 0
00157 { int ii;fprintf(stderr,"\nparset:");
00158   for(ii=0;ii<12;ii++)fprintf(stderr," %g",parvec[ii]); fprintf(stderr,"\n");}
00159 #endif
00160 
00161    /* rotation */
00162 
00163    uu = rot_matrix( 2, D2R*parvec[3] , 0, D2R*parvec[4] , 1, D2R*parvec[5] ) ;
00164 
00165    /* scaling */
00166 
00167    LOAD_DIAG_MAT( dd , parvec[6] , parvec[7] , parvec[8] ) ;
00168 
00169    /* shear */
00170 
00171    switch( smat ){
00172      default:
00173      case SMAT_LOWER:
00174        LOAD_MAT( ss , 1.0        , 0.0        , 0.0 ,
00175                       parvec[9]  , 1.0        , 0.0 ,
00176                       parvec[10] , parvec[11] , 1.0  ) ;
00177      break ;
00178 
00179      case SMAT_UPPER:
00180        LOAD_MAT( ss , 1.0 , parvec[9] , parvec[10] ,
00181                       0.0 , 1.0       , parvec[11] ,
00182                       0.0 , 0.0       , 1.0         ) ;
00183      break ;
00184    }
00185 
00186    /* multiply them, as ordered */
00187 
00188    switch( matorder ){
00189      case MATORDER_SDU:  aa = MAT_MUL(ss,dd) ; bb = uu ; break ;
00190      case MATORDER_SUD:  aa = MAT_MUL(ss,uu) ; bb = dd ; break ;
00191      case MATORDER_DSU:  aa = MAT_MUL(dd,ss) ; bb = uu ; break ;
00192      case MATORDER_DUS:  aa = MAT_MUL(dd,uu) ; bb = ss ; break ;
00193      case MATORDER_USD:  aa = MAT_MUL(uu,ss) ; bb = dd ; break ;
00194      case MATORDER_UDS:  aa = MAT_MUL(uu,dd) ; bb = ss ; break ;
00195    }
00196    mv_for.mm = MAT_MUL(aa,bb) ;
00197 
00198    LOAD_FVEC3( vv , parvec[0] , parvec[1] , parvec[2] ) ;
00199 
00200    switch( dcode ){
00201      case DELTA_AFTER:  mv_for.vv = vv ; break ;
00202      case DELTA_BEFORE: mv_for.vv = MATVEC( mv_for.mm , vv ) ; break ;
00203    }
00204 
00205    mv_inv = INV_VECMAT( mv_for ) ;
00206 
00207 #if 0
00208 DUMP_VECMAT("mv_for",mv_for) ; DUMP_VECMAT("mv_inv",mv_inv) ;
00209 #endif
00210 }

00307 {
00308    THD_mat33 ai ; THD_fvec3 df,di ; int i,j,k ;
00309 
00310    parset_affine() ;  /* sets up numerator matrices: mv_for and mv_inv */
00311 
00312    /* load forward denominator 3-tensor */
00313 
00314    dd_for[0][0][0] = parvec[12]; dd_for[0][0][1] = parvec[13]; dd_for[0][0][2] = parvec[14];
00315    dd_for[0][1][0] = parvec[15]; dd_for[0][1][1] = parvec[16]; dd_for[0][1][2] = parvec[17];
00316    dd_for[0][2][0] = parvec[18]; dd_for[0][2][1] = parvec[19]; dd_for[0][2][2] = parvec[20];
00317    dd_for[1][0][0] = parvec[21]; dd_for[1][0][1] = parvec[22]; dd_for[1][0][2] = parvec[23];
00318    dd_for[1][1][0] = parvec[24]; dd_for[1][1][1] = parvec[25]; dd_for[1][1][2] = parvec[26];
00319    dd_for[1][2][0] = parvec[27]; dd_for[1][2][1] = parvec[28]; dd_for[1][2][2] = parvec[29];
00320    dd_for[2][0][0] = parvec[30]; dd_for[2][0][1] = parvec[31]; dd_for[2][0][2] = parvec[32];
00321    dd_for[2][1][0] = parvec[33]; dd_for[2][1][1] = parvec[34]; dd_for[2][1][2] = parvec[35];
00322    dd_for[2][2][0] = parvec[36]; dd_for[2][2][1] = parvec[37]; dd_for[2][2][2] = parvec[38];
00323 
00324    for( i=0 ; i < 3 ; i++ )
00325     for( j=0 ; j < 3 ; j++ )
00326      for( k=0 ; k < 3 ; k++ ) dd_for[i][j][k] *= dd_fac ;  /* 18 Jul 2005 */
00327 
00328    /* computer inverse denominator 3-tensor */
00329 
00330    ai = mv_inv.mm ;
00331    for( k=0 ; k < 3 ; k++ ){
00332      for( j=0 ; j < 3 ; j++ ){
00333        LOAD_FVEC3( df , -dd_for[0][k][j] , -dd_for[1][k][j] , -dd_for[2][k][j] ) ;
00334        di = MATVEC( ai , df ) ;
00335        UNLOAD_FVEC3( di , dd_inv[0][j][k] , dd_inv[1][j][k] , dd_inv[2][j][k] ) ;
00336      }
00337    }
00338 
00339 #if 0
00340    fprintf(stderr,"++++++++++++ parset_bilinear: dd_fac = %g\n",dd_fac) ;
00341    for( k=0 ; k < 3 ; k++ ){
00342      fprintf(stderr," +-------+ dd_for[][][%d]                       | dd_inv[][][%d]\n"
00343                     "            %11.4g %11.4g %11.4g |  %11.4g %11.4g %11.4g\n"
00344                     "            %11.4g %11.4g %11.4g |  %11.4g %11.4g %11.4g\n"
00345                     "            %11.4g %11.4g %11.4g |  %11.4g %11.4g %11.4g\n" ,
00346          k , k ,
00347          dd_for[0][0][k] , dd_for[0][1][k] , dd_for[0][2][k] ,
00348            dd_inv[0][0][k] , dd_inv[0][1][k] , dd_inv[0][2][k] ,
00349          dd_for[1][0][k] , dd_for[1][1][k] , dd_for[1][2][k] ,
00350            dd_inv[1][0][k] , dd_inv[1][1][k] , dd_inv[1][2][k] ,
00351          dd_for[2][0][k] , dd_for[2][1][k] , dd_for[2][2][k] ,
00352            dd_inv[2][0][k] , dd_inv[2][1][k] , dd_inv[2][2][k]  ) ;
00353    }
00354 #endif
00355 }

00080 {
00081    xsh = parvec[0] ; ysh = parvec[1] ; zsh = parvec[2] ;
00082 }

00126 {
00127    THD_mat33 q , p ;
00128    LOAD_ROT_MAT( q , th1 , ax1 ) ;
00129    LOAD_ROT_MAT( p , th2 , ax2 ) ; q = MAT_MUL( p , q ) ;
00130    LOAD_ROT_MAT( p , th3 , ax3 ) ; q = MAT_MUL( p , q ) ;
00131    return q ;
00132 }

00103 {
00104    THD_fvec3 v , w ;
00105    LOAD_FVEC3( v , aa,bb,cc ) ;
00106    w = VECMAT_VEC( mv_for , v ) ;
00107    *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
00108 }

00112 {
00113    THD_fvec3 v , w ;
00114    LOAD_FVEC3( v , aa,bb,cc ) ;
00115    w = VECMAT_VEC( mv_inv , v ) ;
00116    *p = w.xyz[0] ; *q = w.xyz[1] ; *r = w.xyz[2] ;
00117 }

00267 {
00268    THD_fvec3 x,v,w ; THD_mat33 dd,ee ; int i,j ; float edet,adet,ddet , aa,bb,cc ;
00269    static int first=1 ;
00270 
00271    LOAD_FVEC3( x , ii,jj,kk )   ; v = VECMAT_VEC( ijk_to_xyz , x ) ;
00272    UNLOAD_FVEC3( v , aa,bb,cc ) ; w = VECMAT_VEC( mv_for , v ) ;
00273 
00274    dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
00275    dd.mat[0][1] =        dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
00276    dd.mat[0][2] =        dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
00277    dd.mat[1][0] =        dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
00278    dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
00279    dd.mat[1][2] =        dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
00280    dd.mat[2][0] =        dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
00281    dd.mat[2][1] =        dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
00282    dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
00283 
00284    ddet = MAT_DET(dd) ;
00285    if( first && fabs(ddet) < 0.001f ){
00286      fprintf(stderr,"******* ddet=%g  ii,jj,kk=%g %g %g  aa,bb,cc=%g %g %g\n",
00287                     ddet,ii,jj,kk, aa,bb,cc ) ;
00288      DUMP_MAT33("dd",dd) ;
00289      DUMP_VECMAT("ijk_to_xyz",ijk_to_xyz) ;
00290      DUMP_VECMAT("xyz_to_ijk",xyz_to_ijk) ;
00291      first = 0 ;
00292    }
00293 
00294    ee = MAT_INV(dd) ; edet = MAT_DET(ee) ; v = MATVEC(ee,w) ;
00295 
00296    for( i=0 ; i < 3 ; i++ )
00297     for( j=0 ; j < 3 ; j++ )
00298      dd.mat[i][j] = mv_for.mm.mat[i][j] - dd_for[i][0][j]*v.xyz[0]
00299                                         - dd_for[i][1][j]*v.xyz[1]
00300                                         - dd_for[i][2][j]*v.xyz[2] ;
00301 
00302    adet = MAT_DET(dd) ;
00303    return (adet*edet) ;
00304 }

00220 {
00221    THD_fvec3 v,w ; THD_mat33 dd,ee ;
00222 
00223    LOAD_FVEC3( v , aa,bb,cc ) ;
00224    w = VECMAT_VEC( mv_for , v ) ;
00225 
00226    dd.mat[0][0] = 1.0f + dd_for[0][0][0]*aa + dd_for[0][0][1]*bb + dd_for[0][0][2]*cc ;
00227    dd.mat[0][1] =        dd_for[0][1][0]*aa + dd_for[0][1][1]*bb + dd_for[0][1][2]*cc ;
00228    dd.mat[0][2] =        dd_for[0][2][0]*aa + dd_for[0][2][1]*bb + dd_for[0][2][2]*cc ;
00229    dd.mat[1][0] =        dd_for[1][0][0]*aa + dd_for[1][0][1]*bb + dd_for[1][0][2]*cc ;
00230    dd.mat[1][1] = 1.0f + dd_for[1][1][0]*aa + dd_for[1][1][1]*bb + dd_for[1][1][2]*cc ;
00231    dd.mat[1][2] =        dd_for[1][2][0]*aa + dd_for[1][2][1]*bb + dd_for[1][2][2]*cc ;
00232    dd.mat[2][0] =        dd_for[2][0][0]*aa + dd_for[2][0][1]*bb + dd_for[2][0][2]*cc ;
00233    dd.mat[2][1] =        dd_for[2][1][0]*aa + dd_for[2][1][1]*bb + dd_for[2][1][2]*cc ;
00234    dd.mat[2][2] = 1.0f + dd_for[2][2][0]*aa + dd_for[2][2][1]*bb + dd_for[2][2][2]*cc ;
00235 
00236    ee = MAT_INV(dd) ;
00237    v  = MATVEC(ee,w) ;
00238 
00239    *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
00240 }

00244 {
00245    THD_fvec3 v,w ; THD_mat33 dd,ee ;
00246 
00247    LOAD_FVEC3( v , aa,bb,cc ) ;
00248    w = VECMAT_VEC( mv_inv , v ) ;
00249 
00250    dd.mat[0][0] = 1.0f + dd_inv[0][0][0]*aa + dd_inv[0][0][1]*bb + dd_inv[0][0][2]*cc ;
00251    dd.mat[0][1] =        dd_inv[0][1][0]*aa + dd_inv[0][1][1]*bb + dd_inv[0][1][2]*cc ;
00252    dd.mat[0][2] =        dd_inv[0][2][0]*aa + dd_inv[0][2][1]*bb + dd_inv[0][2][2]*cc ;
00253    dd.mat[1][0] =        dd_inv[1][0][0]*aa + dd_inv[1][0][1]*bb + dd_inv[1][0][2]*cc ;
00254    dd.mat[1][1] = 1.0f + dd_inv[1][1][0]*aa + dd_inv[1][1][1]*bb + dd_inv[1][1][2]*cc ;
00255    dd.mat[1][2] =        dd_inv[1][2][0]*aa + dd_inv[1][2][1]*bb + dd_inv[1][2][2]*cc ;
00256    dd.mat[2][0] =        dd_inv[2][0][0]*aa + dd_inv[2][0][1]*bb + dd_inv[2][0][2]*cc ;
00257    dd.mat[2][1] =        dd_inv[2][1][0]*aa + dd_inv[2][1][1]*bb + dd_inv[2][1][2]*cc ;
00258    dd.mat[2][2] = 1.0f + dd_inv[2][2][0]*aa + dd_inv[2][2][1]*bb + dd_inv[2][2][2]*cc ;
00259 
00260    ee = MAT_INV(dd) ;
00261    v  = MATVEC(ee,w) ;
00262 
00263    *p = v.xyz[0] ; *q = v.xyz[1] ; *r = v.xyz[2] ;
00264 }

00086 {
00087    *p = aa+xsh ; *q = bb+ysh ; *r = cc+zsh ;
00088 }

00092 {
00093    *p = aa-xsh ; *q = bb-ysh ; *r = cc-zsh ;
00094 }