00001 /*! Dealings with volume data */
00002 #include "SUMA_suma.h"
00003 
00004 extern SUMA_CommonFields *SUMAg_CF; 
00005 
00006 /*! a copy of THD_handedness from ../thd_rotangles.c
00007 Dunno why the original was giving me pain linking ... */
00008 int SUMA_THD_handedness( THD_3dim_dataset * dset )
00009 {
00010    THD_dataxes * dax ;
00011    THD_mat33 q ;
00012    int col ;
00013    float val ;
00014 
00015 
00016    if( !ISVALID_DSET(dset) ) SUMA_RETURN(1) ;
00017 
00018    LOAD_ZERO_MAT(q) ;
00019    dax = dset->daxes ;
00020 
00021    col = 0 ;
00022    switch( dax->xxorient ){
00023       case 0: q.mat[0][col] =  1.0 ; break ;
00024       case 1: q.mat[0][col] = -1.0 ; break ;
00025       case 2: q.mat[1][col] = -1.0 ; break ;
00026       case 3: q.mat[1][col] =  1.0 ; break ;
00027       case 4: q.mat[2][col] =  1.0 ; break ;
00028       case 5: q.mat[2][col] = -1.0 ; break ;
00029    }
00030 
00031    col = 1 ;
00032    switch( dax->yyorient ){
00033       case 0: q.mat[0][col] =  1.0 ; break ;
00034       case 1: q.mat[0][col] = -1.0 ; break ;
00035       case 2: q.mat[1][col] = -1.0 ; break ;
00036       case 3: q.mat[1][col] =  1.0 ; break ;
00037       case 4: q.mat[2][col] =  1.0 ; break ;
00038       case 5: q.mat[2][col] = -1.0 ; break ;
00039    }
00040 
00041    col = 2 ;
00042    switch( dax->zzorient ){
00043       case 0: q.mat[0][col] =  1.0 ; break ;
00044       case 1: q.mat[0][col] = -1.0 ; break ;
00045       case 2: q.mat[1][col] = -1.0 ; break ;
00046       case 3: q.mat[1][col] =  1.0 ; break ;
00047       case 4: q.mat[2][col] =  1.0 ; break ;
00048       case 5: q.mat[2][col] = -1.0 ; break ;
00049    }
00050 
00051    val = MAT_DET(q) ;
00052    if( val > 0.0 ) SUMA_RETURN( 1) ;  /* right handed */
00053    else            SUMA_RETURN(-1) ;  /* left handed */
00054 }
00055 
00056 /*!
00057    see help for SUMA_AfniPrefix below
00058 */
00059 SUMA_Boolean SUMA_AfniExistsView(int exists, char *view)
00060 {
00061    static char FuncName[]={"SUMA_AfniExistsView"};
00062    SUMA_Boolean ans = NOPE;
00063    
00064    if (!exists) SUMA_RETURN(ans);
00065    
00066    if (strcmp(view,"+orig") == 0) {
00067       if (exists == 1 || exists == 3 || exists == 5 || exists == 7) ans = YUP; 
00068    } else if (strcmp(view,"+acpc") == 0) {
00069       if (exists == 2 || exists == 3 || exists == 6 || exists == 7) ans = YUP; 
00070    } else if (strcmp(view,"+tlrc") == 0) {
00071       if (exists == 4 || exists == 5 || exists == 6 || exists == 7) ans = YUP; 
00072    } 
00073    
00074    SUMA_RETURN(ans);
00075 
00076 }
00077 
00078 SUMA_Boolean SUMA_AfniView (char *nameorig, char *cview)
00079 {
00080    static char FuncName[]={"SUMA_AfniView"};
00081    char *tmp1 = NULL, *tmp2 = NULL;
00082    SUMA_Boolean LocalHead = NOPE;
00083    
00084    SUMA_ENTRY;
00085    
00086    if (!nameorig) SUMA_RETURN(NOPE);
00087    if (!cview) SUMA_RETURN(NOPE);
00088 
00089    tmp1 = SUMA_Extension(nameorig, ".HEAD", YUP);
00090    tmp2 = SUMA_Extension(tmp1, ".BRIK", YUP); SUMA_free(tmp1); tmp1 = NULL;
00091    /* is there a dot ?*/
00092    if (tmp2[strlen(tmp2)-1] == '.') tmp2[strlen(tmp2)-1] = '\0';
00093    
00094    if (LocalHead) fprintf(SUMA_STDERR,"%s: Searching for view of %s\n", FuncName, tmp2);
00095    
00096    /* view */
00097    if (SUMA_isExtension(tmp2, "+orig")) { 
00098       sprintf(cview, "+orig"); 
00099    } else if (SUMA_isExtension(tmp2, "+acpc")) { 
00100       sprintf(cview, "+acpc"); 
00101    } else if (SUMA_isExtension(tmp2, "+tlrc")) { 
00102       sprintf(cview, "+tlrc"); 
00103    } else {
00104       cview[0]='\0';
00105    }
00106    SUMA_free(tmp2); tmp2 = NULL;
00107 
00108    SUMA_RETURN(YUP);
00109 }
00110 
00111 SUMA_Boolean SUMA_AfniExists(char *prefix, char *c2view) 
00112 {
00113    static char FuncName[]={"SUMA_AfniExists"};
00114    char *head=NULL;
00115    SUMA_Boolean ans = NOPE;
00116    SUMA_Boolean LocalHead = NOPE;
00117    
00118    ans = NOPE;
00119 
00120    head = SUMA_append_replace_string(prefix,".HEAD", c2view,0);
00121    if (LocalHead) fprintf(SUMA_STDERR,"%s: Checking existence of %s\n", FuncName, head);
00122    if (SUMA_filexists(head)) { ans = YUP; }
00123    SUMA_free(head); head = NULL; 
00124    
00125    SUMA_RETURN(ans);
00126 }
00127 /*!
00128    \brief Return AFNI's prefix, from a file name also checks for its validity
00129    \param name (char *) dset name (can contain path)
00130    \param view (char[4]) array to return view in it
00131    \param path (char *) if any, make sure it is not duplicated in name!
00132    \param exists (int *) function checks
00133                         for all three possible views
00134                         0: No afni dset with name/prefix
00135                         1: +orig
00136                         2: +acpc
00137                            3: +orig and +acpc 
00138                         4: +tlrc
00139                            5: +orig and +tlrc
00140                            6: +acpc and +tlrc
00141                            7: +orig and +acpc and +tlrc
00142                         \sa SUMA_AfniExistsView   
00143                         
00144    \return prefix (char *) dset prefix, free it with SUMA_free
00145    \sa SUMA_AfniView
00146 */
00147 char *SUMA_AfniPrefix(char *nameorig, char *view, char *path, int *exists) 
00148 {
00149    static char FuncName[]={"SUMA_AfniPrefix"};
00150    char *tmp1 = NULL, *tmp2 = NULL, *prfx = NULL, *name=NULL;
00151    char cview[10];
00152    int iview;
00153    SUMA_Boolean LocalHead = NOPE;
00154    
00155    SUMA_ENTRY;
00156    
00157    if (!nameorig) SUMA_RETURN(prfx);
00158    if (exists) *exists = 0;
00159    
00160    if (LocalHead) fprintf(SUMA_STDERR,"%s: Working with %s\n", FuncName, nameorig);
00161    
00162    if (!path) name = SUMA_copy_string(nameorig);
00163    else {
00164       if (path[strlen(path)-1] == '/') {
00165          name = SUMA_append_replace_string(path, nameorig, "", 0);
00166       } else {
00167          name = SUMA_append_replace_string(path, nameorig, "/", 0);
00168       }
00169    }
00170    
00171    if (LocalHead) fprintf(SUMA_STDERR,"%s: name now %s\n", FuncName, name);
00172    
00173    tmp1 = SUMA_Extension(name, ".HEAD", YUP);
00174    tmp2 = SUMA_Extension(tmp1, ".BRIK", YUP); SUMA_free(tmp1); tmp1 = NULL;
00175    /* is there a dot ?*/
00176    if (tmp2[strlen(tmp2)-1] == '.') tmp2[strlen(tmp2)-1] = '\0';
00177    
00178    if (LocalHead) fprintf(SUMA_STDERR,"%s: Searching for view of %s\n", FuncName, tmp2);
00179    
00180    /* view */
00181    iview = -1;
00182    if (SUMA_isExtension(tmp2, "+orig")) { 
00183       iview = 0; sprintf(cview, "+orig"); 
00184       prfx = SUMA_Extension(tmp2, "+orig", YUP);
00185    } else if (SUMA_isExtension(tmp2, "+acpc")) { 
00186       iview = 1; sprintf(cview, "+acpc"); 
00187       prfx = SUMA_Extension(tmp2, "+acpc", YUP);
00188    } else if (SUMA_isExtension(tmp2, "+tlrc")) { 
00189       iview = 2; sprintf(cview, "+tlrc"); 
00190       prfx = SUMA_Extension(tmp2, "+tlrc", YUP);
00191    } else {
00192       prfx = SUMA_copy_string(tmp2);
00193       cview[0]='\0';
00194    }
00195    SUMA_free(tmp2); tmp2 = NULL;
00196    
00197    if (LocalHead) fprintf(SUMA_STDERR,"%s: Prefix %s\n", FuncName, prfx);
00198    
00199    /* can't tell what view is, so can't test properly quite yet */
00200    {   
00201       /* is name ok ? all I have to do is test with +orig */
00202       char *bname=SUMA_append_string(prfx,"+orig");
00203       if (LocalHead) fprintf(SUMA_STDERR,"%s: Checking name validity using +orig for view %s\n", FuncName, bname);
00204       if( !THD_filename_ok(bname) ){
00205          SUMA_SL_Err("not a proper dset name!\n") ;
00206          SUMA_free(name); name = NULL;
00207          SUMA_free(prfx); prfx = NULL;
00208          SUMA_free(bname); bname = NULL;
00209          SUMA_RETURN(prfx);
00210       }
00211    }
00212    
00213    if (exists) {
00214       { /* look for any */
00215          char *head=NULL, c2view[10];
00216          iview = 0; *exists = 0;
00217          while (iview < 3) {
00218             if (iview == 0) sprintf(c2view, "+orig"); 
00219             else if (iview == 1) sprintf(c2view, "+acpc"); 
00220             else if (iview == 2) sprintf(c2view, "+tlrc");
00221             head = SUMA_append_replace_string(prfx,".HEAD", c2view,0);
00222             if (LocalHead) fprintf(SUMA_STDERR,"%s: Checking existence of %s\n", FuncName, head);
00223             if (SUMA_filexists(head)) { *exists += (int)pow(2,iview); }
00224             SUMA_free(head); head = NULL; 
00225             ++iview;
00226          }
00227       }
00228       if (LocalHead) fprintf(SUMA_STDERR,"%s: exist number is %d\n", FuncName, *exists);
00229    }
00230    
00231    
00232    if (cview[0] != '\0') {
00233       if (LocalHead) {
00234          if (SUMA_AfniExistsView(*exists, cview)) {
00235             fprintf(SUMA_STDERR,"%s: dset with view %s does exist.\n", FuncName, cview); 
00236          } else {
00237             fprintf(SUMA_STDERR,"%s: dset with view %s does not exist.\n", FuncName, cview); 
00238          }
00239       }
00240    }
00241    
00242    if (view) {
00243       sprintf(view,"%s", cview);
00244    }
00245    if (name) SUMA_free(name); name = NULL;
00246    
00247    SUMA_RETURN(prfx);
00248 }
00249                         
00250 /*!
00251    \brief A function to find the skin of a volume
00252    \param dset (THD_3dim_dataset *) an AFNI volume
00253    \param dvec (double *) (nx * ny * nz) data vector
00254    \param thresh (double) consider only values in dvec > thresh
00255    \param N_skin (int *) number of voxels that are skin
00256    \return skin (byte *) (nx * ny * nz) vector containing 1 for skin voxels, 0 elsewhere.
00257 */
00258 byte * SUMA_isSkin(THD_3dim_dataset *dset, double *dvec, double thresh, int *N_skin)
00259 {
00260    static char FuncName[]={"SUMA_isSkin"};
00261    byte *isskin=NULL;
00262    int x, y, z, nx, ny, nz, i1D,  nxy;
00263   
00264    SUMA_ENTRY;
00265    
00266    if (!dset || !dvec) {
00267       SUMA_SL_Err("NULL input dset or dvec");
00268       SUMA_RETURN(isskin);
00269    }
00270    
00271    nx = DSET_NX(dset);
00272    ny = DSET_NY(dset);
00273    nz = DSET_NZ(dset);
00274    nxy = nx * ny;
00275    
00276    isskin = (byte *) SUMA_calloc(nxy * nz, sizeof(byte));
00277    if (!isskin) {
00278       SUMA_SL_Crit("Failed to allocate");
00279       SUMA_RETURN(NULL);
00280    }
00281    
00282    *N_skin = 0;
00283    for (z=0; z<nz; ++z) {
00284       for (y=0; y<ny; ++y) {
00285          x = 0;
00286          do { /* the upstroke */
00287             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00288             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00289             ++x; 
00290          } while (x < nx && !isskin[i1D]);
00291          x = nx - 1;
00292          do { /* the downstroke */
00293             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00294             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00295             --x; 
00296          } while (x >=0 && !isskin[i1D]);
00297       } /* y */
00298    } /* z */
00299    
00300    for (z=0; z<nz; ++z) {
00301       for (x=0; x<nx; ++x) {
00302          y = 0;
00303          do { /* the upstroke */
00304             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00305             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00306             ++y; 
00307          } while (y < ny && !isskin[i1D]);
00308          y = ny - 1;
00309          do { /* the downstroke */
00310             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00311             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00312             --y; 
00313          } while (y >=0 && !isskin[i1D]);
00314       } /* x */
00315    } /* z */
00316    
00317    for (x=0; x<nx; ++x) {
00318       for (y=0; y<ny; ++y) {
00319          z = 0;
00320          do { /* the upstroke */
00321             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00322             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00323             ++z; 
00324          } while (z < nz && !isskin[i1D]);
00325          z = nz - 1;
00326          do { /* the downstroke */
00327             i1D = SUMA_3D_2_1D_index(x, y, z, nx, nxy);
00328             if (dvec[i1D] > thresh) { isskin[i1D] = 1; ++ *N_skin; }
00329             --z; 
00330          } while (z >=0 && !isskin[i1D]);
00331       } /* y */
00332    } /* x */
00333       
00334    SUMA_RETURN(isskin);
00335 }
00336 /*!
00337    Returns the attributes of the surface's volume parent 
00338    A volume parent is defined as:
00339       The volume data set used to generate the surface. aka the master SPGR
00340    \ret NULL for troubles
00341 */
00342 SUMA_VOLPAR *SUMA_Alloc_VolPar (void)
00343 {
00344    static char FuncName[]={"SUMA_Alloc_VolPar"};
00345    SUMA_VOLPAR *VP;
00346 
00347    SUMA_ENTRY;
00348 
00349    VP = (SUMA_VOLPAR *)SUMA_malloc(sizeof(SUMA_VOLPAR));
00350    if (VP == NULL) {
00351       fprintf(SUMA_STDERR,"Error SUMA_Alloc_VolPar: Failed to allocate for VolPar\n");
00352       SUMA_RETURN (NULL);
00353    }
00354    VP->idcode_str = NULL;
00355    VP->isanat = 1;
00356    VP->nx = VP->ny = VP->nz = 0; /*!< number of voxels in the three dimensions */
00357    VP->dx = VP->dy = VP->dz = 0.0; /*!< delta x, y, z in mm */
00358    VP->xorg = VP->yorg = VP->zorg = 0.0; /*!< voxel origin in three dimensions */
00359    VP->prefix = NULL; /*!< parent volume prefix */
00360    VP->filecode = NULL; /*!< parent volume prefix + view */
00361    VP->dirname = NULL; /*!< parent volume directory name */
00362    VP->vol_idcode_str = NULL; /*!< idcode string OF parent volume*/
00363    VP->vol_idcode_date = NULL; /*!< idcode date */
00364    VP->xxorient = VP->yyorient = VP->zzorient = 0; /*!< orientation of three dimensions*/ 
00365    VP->VOLREG_CENTER_OLD = NULL; /*!< pointer to the named attribute (3x1) in the .HEAD file of the experiment-aligned Parent Volume */
00366    VP->VOLREG_CENTER_BASE = NULL; /*!< pointer to the named attribute (3x1) in the .HEAD file of the experiment-aligned Parent Volume */
00367    VP->VOLREG_MATVEC = NULL; /*!< pointer to the named attribute (12x1) in the .HEAD file of the experiment-aligned Parent Volume */
00368    VP->TAGALIGN_MATVEC = NULL; /*!< pointer to the named attribute (12x1) in the .HEAD file of the tag aligned Parent Volume */
00369    VP->ROTATE_MATVEC = NULL; /*!< pointer to the named attribute (12x1) in the .HEAD file of the 3drotated Parent Volume */
00370    VP->ROTATE_CENTER_OLD = NULL;
00371    VP->ROTATE_CENTER_BASE = NULL;
00372    VP->Hand = 1; /*!< Handedness of axis 1 RH, -1 LH*/
00373    
00374    SUMA_RETURN(VP);
00375 }
00376 SUMA_Boolean SUMA_Free_VolPar (SUMA_VOLPAR *VP)
00377 {
00378    static char FuncName[]={"SUMA_Free_VolPar"};
00379    
00380    SUMA_ENTRY;
00381 
00382    if (VP->prefix != NULL) SUMA_free(VP->prefix);
00383    if (VP->idcode_str != NULL) SUMA_free(VP->idcode_str);
00384    if (VP->filecode != NULL) SUMA_free(VP->filecode);
00385    if (VP->dirname != NULL) SUMA_free(VP->dirname);
00386    if (VP->vol_idcode_str != NULL) SUMA_free(VP->vol_idcode_str);
00387    if (VP->vol_idcode_date != NULL) SUMA_free(VP->vol_idcode_date);
00388    if (VP->VOLREG_CENTER_OLD != NULL) SUMA_free(VP->VOLREG_CENTER_OLD);
00389    if (VP->VOLREG_CENTER_BASE != NULL) SUMA_free(VP->VOLREG_CENTER_BASE);
00390    if (VP->VOLREG_MATVEC != NULL) SUMA_free(VP->VOLREG_MATVEC);
00391    if (VP->TAGALIGN_MATVEC != NULL) SUMA_free(VP->TAGALIGN_MATVEC);
00392    if (VP->ROTATE_MATVEC != NULL) SUMA_free(VP->ROTATE_MATVEC);
00393    if (VP->ROTATE_CENTER_OLD != NULL) SUMA_free(VP->ROTATE_CENTER_OLD);
00394    if (VP->ROTATE_CENTER_BASE != NULL) SUMA_free(VP->ROTATE_CENTER_BASE);
00395    if (VP != NULL) SUMA_free(VP);
00396    SUMA_RETURN (YUP);
00397 }
00398 
00399 SUMA_VOLPAR *SUMA_VolParFromDset (THD_3dim_dataset *dset)
00400 {
00401    ATR_float *atr;
00402    static char FuncName[]={"SUMA_VolParFromDset"};
00403    SUMA_VOLPAR *VP;
00404    int ii;
00405    MCW_idcode idcode;
00406    SUMA_Boolean LocalHead = NOPE;
00407       
00408    SUMA_ENTRY;
00409 
00410    /* read the header of the parent volume */
00411    if (dset == NULL) {
00412       fprintf (SUMA_STDERR,"Error %s:\nNULL dset\n", FuncName);
00413       SUMA_RETURN (NULL);
00414    }
00415    
00416    /* allocate for VP */
00417    VP = SUMA_Alloc_VolPar();
00418    if (VP == NULL) {
00419       fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Alloc_VolPar\n", FuncName);
00420       SUMA_RETURN (NULL);
00421    }
00422    
00423    /* retain pertinent info */
00424    VP->isanat = ISANAT(dset);
00425    VP->nx = DSET_NX(dset);
00426    VP->ny = DSET_NY(dset);
00427    VP->nz = DSET_NZ(dset);
00428    VP->dx = DSET_DX(dset);
00429    VP->dy = DSET_DY(dset);
00430    VP->dz = DSET_DZ(dset);
00431    VP->xorg = DSET_XORG(dset);
00432    VP->yorg = DSET_YORG(dset);
00433    VP->zorg = DSET_ZORG(dset);
00434    ii = strlen(DSET_PREFIX(dset));
00435    VP->prefix = (char *)SUMA_malloc(ii+1);
00436    ii = strlen(DSET_FILECODE(dset));
00437    VP->filecode = (char *)SUMA_malloc(ii+1);
00438    ii = strlen(DSET_DIRNAME(dset));
00439    VP->dirname = (char *)SUMA_malloc(ii+1);
00440    ii = strlen(dset->idcode.str);
00441    VP->vol_idcode_str = (char *)SUMA_malloc(ii+1);
00442    ii = strlen(dset->idcode.date);
00443    VP->vol_idcode_date = (char *)SUMA_malloc(ii+1);
00444    if (VP->prefix == NULL || VP->filecode == NULL || VP->vol_idcode_date == NULL || VP->dirname == NULL || VP->vol_idcode_str == NULL) {
00445       fprintf(SUMA_STDERR,"Error %s: Failed to allocate for strings. Kill me, please.\n", FuncName);
00446       SUMA_Free_VolPar(VP);
00447       SUMA_RETURN (NULL);
00448    }
00449    VP->prefix = strcpy(VP->prefix, DSET_PREFIX(dset));
00450    VP->filecode = strcpy(VP->filecode, DSET_FILECODE(dset));
00451    VP->dirname = strcpy(VP->dirname, DSET_DIRNAME(dset));
00452    VP->vol_idcode_str = strcpy(VP->vol_idcode_str, dset->idcode.str);
00453    VP->vol_idcode_date = strcpy(VP->vol_idcode_date, dset->idcode.date);
00454    VP->xxorient = dset->daxes->xxorient;
00455    VP->yyorient = dset->daxes->yyorient;
00456    VP->zzorient = dset->daxes->zzorient;
00457 
00458    if (LocalHead) {      
00459       fprintf (SUMA_STDERR,"%s: dset->idcode_str = %s\n", FuncName, dset->idcode.str);
00460       fprintf (SUMA_STDERR,"%s: VP->vol_idcode_str = %s\n", FuncName, VP->vol_idcode_str);
00461    }
00462    /* Get the 3drotate matrix if possible*/
00463    atr = THD_find_float_atr( dset->dblk , "ROTATE_MATVEC_000000" ) ;
00464    if (atr == NULL) {
00465       VP->ROTATE_MATVEC = NULL;
00466    }else {
00467       VP->ROTATE_MATVEC = (float *)SUMA_calloc(12, sizeof(float));
00468       if (VP->ROTATE_MATVEC != NULL) {
00469          if (atr->nfl == 12) {
00470             for (ii=0; ii<12; ++ii) VP->ROTATE_MATVEC[ii] = atr->fl[ii];
00471          } else {   
00472             fprintf(SUMA_STDERR,"Error %s: ROTATE_MATVEC does not have 12 elements.\n", FuncName);
00473          }
00474       } else {
00475          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->ROTATE_MATVEC\n", FuncName);
00476       }
00477    }
00478    
00479    /* Get the tagalign matrix if possible*/
00480    atr = THD_find_float_atr( dset->dblk , "TAGALIGN_MATVEC" ) ;
00481    if (atr == NULL) {
00482       VP->TAGALIGN_MATVEC = NULL;
00483    }else {
00484       VP->TAGALIGN_MATVEC = (float *)SUMA_calloc(12, sizeof(float));
00485       if (VP->TAGALIGN_MATVEC != NULL) {
00486          if (atr->nfl == 12) {
00487             for (ii=0; ii<12; ++ii) VP->TAGALIGN_MATVEC[ii] = atr->fl[ii];
00488          } else {   
00489             fprintf(SUMA_STDERR,"Error %s: TAGALIGN_MATVEC does not have 12 elements.\n", FuncName);
00490          }
00491       } else {
00492          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->TAGALIGN_MATVEC\n", FuncName);
00493       }
00494    }
00495    
00496    /* Get the volreg matrix if possible*/
00497    atr = THD_find_float_atr( dset->dblk , "VOLREG_MATVEC_000000" ) ;
00498    if (atr == NULL) {
00499       VP->VOLREG_MATVEC = NULL;
00500    }else {
00501       VP->VOLREG_MATVEC = (float *)SUMA_calloc(12, sizeof(float));
00502       if (VP->VOLREG_MATVEC != NULL) {
00503          if (atr->nfl == 12) {
00504             for (ii=0; ii<12; ++ii) VP->VOLREG_MATVEC[ii] = atr->fl[ii];
00505          } else {   
00506             fprintf(SUMA_STDERR,"Error %s: VOLREG_MATVEC_000000 does not have 12 elements.\n", FuncName);
00507          }
00508       } else {
00509          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_MATVEC\n", FuncName);
00510       }
00511    }
00512    /* Get the center base coordinates */
00513    atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_BASE");
00514    if (atr == NULL) {
00515       VP->VOLREG_CENTER_BASE = NULL;
00516    } else {
00517       VP->VOLREG_CENTER_BASE = (float *)SUMA_calloc(3, sizeof(float));
00518       if (VP->VOLREG_CENTER_BASE != NULL) {
00519          if (atr->nfl == 3) {
00520             for (ii=0; ii<3; ++ii) VP->VOLREG_CENTER_BASE[ii] = atr->fl[ii];
00521          } else {   
00522             fprintf(SUMA_STDERR,"Error %s: VOLREG_CENTER_BASE does not have 12 elements.\n", FuncName);
00523          }
00524       } else {
00525          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_CENTER_BASE\n", FuncName);
00526       }
00527    }
00528    /* VOLREG_CENTER_OLD  */
00529    atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_OLD");
00530    if (atr == NULL) {
00531       VP->VOLREG_CENTER_OLD = NULL;
00532    } else {
00533       VP->VOLREG_CENTER_OLD = (float *)SUMA_calloc(3, sizeof(float));
00534       if (VP->VOLREG_CENTER_OLD != NULL) {
00535          if (atr->nfl == 3) {
00536             for (ii=0; ii<3; ++ii) VP->VOLREG_CENTER_OLD[ii] = atr->fl[ii];
00537          } else {   
00538             fprintf(SUMA_STDERR,"Error %s: VOLREG_CENTER_OLD does not have 12 elements.\n", FuncName);
00539          }
00540       } else {
00541          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_CENTER_OLD\n", FuncName);
00542       }
00543    }
00544    
00545    /* Get the center base coordinates */
00546    atr = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_BASE");
00547    if (atr == NULL) {
00548       VP->ROTATE_CENTER_BASE = NULL;
00549    } else {
00550       VP->ROTATE_CENTER_BASE = (float *)SUMA_calloc(3, sizeof(float));
00551       if (VP->ROTATE_CENTER_BASE != NULL) {
00552          if (atr->nfl == 3) {
00553             for (ii=0; ii<3; ++ii) VP->ROTATE_CENTER_BASE[ii] = atr->fl[ii];
00554          } else {   
00555             fprintf(SUMA_STDERR,"Error %s: ROTATE_CENTER_BASE does not have 12 elements.\n", FuncName);
00556          }
00557       } else {
00558          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->ROTATE_CENTER_BASE\n", FuncName);
00559       }
00560    }
00561    /* ROTATE_CENTER_OLD */
00562    atr = THD_find_float_atr( dset->dblk , "ROTATE_CENTER_OLD");
00563    if (atr == NULL) {
00564       VP->ROTATE_CENTER_OLD = NULL;
00565    } else {
00566       VP->ROTATE_CENTER_OLD = (float *)SUMA_calloc(3, sizeof(float));
00567       if (VP->ROTATE_CENTER_OLD != NULL) {
00568          if (atr->nfl == 3) {
00569             for (ii=0; ii<3; ++ii) VP->ROTATE_CENTER_OLD[ii] = atr->fl[ii];
00570          } else {   
00571             fprintf(SUMA_STDERR,"Error %s: ROTATE_CENTER_OLD does not have 12 elements.\n", FuncName);
00572          }
00573       } else {
00574          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->ROTATE_CENTER_OLD\n", FuncName);
00575       }
00576    }
00577 
00578    /* handedness */
00579    VP->Hand = SUMA_THD_handedness( dset );
00580    
00581    SUMA_RETURN (VP);
00582 }
00583 
00584 SUMA_VOLPAR *SUMA_VolPar_Attr (char *volparent_name)
00585 {
00586    ATR_float *atr;
00587    static char FuncName[]={"SUMA_VolPar_Attr"};
00588    SUMA_VOLPAR *VP;
00589    THD_3dim_dataset *dset;
00590    int ii;
00591    MCW_idcode idcode;
00592    SUMA_Boolean LocalHead = NOPE;
00593       
00594    SUMA_ENTRY;
00595 
00596    /* read the header of the parent volume */
00597    dset = THD_open_dataset(volparent_name);
00598    if (dset == NULL) {
00599       fprintf (SUMA_STDERR,"Error %s: Could not read %s\n", FuncName, volparent_name);
00600       SUMA_RETURN (NULL);
00601    }
00602    
00603    VP = SUMA_VolParFromDset(dset);
00604    if (!VP) {
00605       SUMA_SL_Err("Failed in SUMA_VolParFromDset");
00606    }
00607 
00608    /* now free the dset pointer */
00609    THD_delete_3dim_dataset( dset , False ) ;
00610    
00611    SUMA_RETURN (VP);
00612 }
00613 
00614 /*!
00615    \brief Form a string containing the info of the volume parent
00616    
00617    \param VP (SUMA_VOLPAR *) Volume parent structure.
00618    \return s (char *) pointer to NULL terminated string containing surface info.
00619    It is your responsability to free it.
00620    
00621    \sa SUMA_Show_VolPar
00622 */
00623 
00624 char *SUMA_VolPar_Info (SUMA_VOLPAR *VP)
00625 {
00626    static char FuncName[]={"SUMA_VolPar_Info"};
00627    char stmp[1000], *s=NULL;
00628    SUMA_STRING *SS = NULL;
00629    
00630    SUMA_ENTRY;
00631 
00632    SS = SUMA_StringAppend (NULL, NULL);
00633    
00634    if (VP) { 
00635       sprintf (stmp,"\nVP contents:\n");
00636       SS = SUMA_StringAppend (SS, stmp);
00637       sprintf (stmp,"prefix: %s\tfilecode: %s\tdirname: %s\nId code str:%s\tID code date: %s\n", \
00638          VP->prefix, VP->filecode, VP->dirname, VP->vol_idcode_str, VP->vol_idcode_date);
00639       SS = SUMA_StringAppend (SS, stmp);
00640       if (VP->idcode_str) SS = SUMA_StringAppend (SS, "IDcode is NULL\n");
00641       else SS = SUMA_StringAppend_va (SS, "IDcode: %s\n", VP->idcode_str);
00642       
00643       sprintf (stmp,"isanat: %d\n", VP->isanat);
00644       SS = SUMA_StringAppend (SS, stmp);
00645       sprintf (stmp,"Orientation: %d %d %d\n", \
00646          VP->xxorient, VP->yyorient, VP->zzorient);
00647       if (VP->Hand == 1) SS = SUMA_StringAppend (SS, "Right Hand Coordinate System.\n");
00648       else if (VP->Hand == -1) SS = SUMA_StringAppend (SS, "Left Hand Coordinate System.\n");
00649       else SS = SUMA_StringAppend (SS, "No hand coordinate system!\n");
00650       
00651       SS = SUMA_StringAppend (SS, stmp);
00652       sprintf (stmp,"Origin: %f %f %f\n", \
00653          VP->xorg, VP->yorg, VP->zorg);
00654       SS = SUMA_StringAppend (SS, stmp);
00655       sprintf (stmp,"Delta: %f %f %f\n", \
00656          VP->dx, VP->dy, VP->dz);
00657       SS = SUMA_StringAppend (SS, stmp);
00658       sprintf (stmp,"N: %d %d %d\n",\
00659          VP->nx, VP->ny, VP->nz);
00660       SS = SUMA_StringAppend (SS, stmp);
00661 
00662       if (VP->ROTATE_MATVEC != NULL) {
00663          sprintf (stmp,"VP->ROTATE_MATVEC = \n\tMrot\tDelta\n");
00664          SS = SUMA_StringAppend (SS, stmp);
00665          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00666          VP->ROTATE_MATVEC[0], VP->ROTATE_MATVEC[1], VP->ROTATE_MATVEC[2], VP->ROTATE_MATVEC[3]); 
00667          SS = SUMA_StringAppend (SS, stmp);
00668          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00669          VP->ROTATE_MATVEC[4], VP->ROTATE_MATVEC[5], VP->ROTATE_MATVEC[6], VP->ROTATE_MATVEC[7]);
00670          SS = SUMA_StringAppend (SS, stmp);
00671          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00672          VP->ROTATE_MATVEC[8], VP->ROTATE_MATVEC[9], VP->ROTATE_MATVEC[10], VP->ROTATE_MATVEC[11]);
00673          SS = SUMA_StringAppend (SS, stmp);
00674       } else {
00675          sprintf (stmp,"VP->ROTATE_MATVEC = NULL\n");
00676          SS = SUMA_StringAppend (SS, stmp);
00677       }      
00678       if (VP->ROTATE_CENTER_OLD != NULL) {
00679          sprintf (stmp,"VP->ROTATE_CENTER_OLD = %f, %f, %f\n", \
00680            VP->ROTATE_CENTER_OLD[0], VP->ROTATE_CENTER_OLD[1], VP->ROTATE_CENTER_OLD[2]); 
00681          SS = SUMA_StringAppend (SS, stmp);
00682       }else {
00683          sprintf (stmp,"VP->ROTATE_CENTER_OLD = NULL\n");
00684          SS = SUMA_StringAppend (SS, stmp);
00685       }
00686       if (VP->ROTATE_CENTER_BASE != NULL) {
00687          sprintf (stmp,"VP->ROTATE_CENTER_BASE = %f, %f, %f\n", \
00688             VP->ROTATE_CENTER_BASE[0], VP->ROTATE_CENTER_BASE[1], VP->ROTATE_CENTER_BASE[2]); 
00689          SS = SUMA_StringAppend (SS, stmp);
00690       } else {
00691          sprintf (stmp,"VP->ROTATE_CENTER_BASE = NULL\n");
00692          SS = SUMA_StringAppend (SS, stmp);
00693       }      
00694       if (VP->TAGALIGN_MATVEC != NULL) {
00695          sprintf (stmp,"VP->TAGALIGN_MATVEC = \n\tMrot\tDelta\n");
00696          SS = SUMA_StringAppend (SS, stmp);
00697          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00698          VP->TAGALIGN_MATVEC[0], VP->TAGALIGN_MATVEC[1], VP->TAGALIGN_MATVEC[2], VP->TAGALIGN_MATVEC[3]); 
00699          SS = SUMA_StringAppend (SS, stmp);
00700          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00701          VP->TAGALIGN_MATVEC[4], VP->TAGALIGN_MATVEC[5], VP->TAGALIGN_MATVEC[6], VP->TAGALIGN_MATVEC[7]);
00702          SS = SUMA_StringAppend (SS, stmp);
00703          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00704          VP->TAGALIGN_MATVEC[8], VP->TAGALIGN_MATVEC[9], VP->TAGALIGN_MATVEC[10], VP->TAGALIGN_MATVEC[11]);
00705          SS = SUMA_StringAppend (SS, stmp);
00706       } else {
00707          sprintf (stmp,"VP->TAGALIGN_MATVEC = NULL\n");
00708          SS = SUMA_StringAppend (SS, stmp);
00709       }      
00710       
00711       if (VP->VOLREG_MATVEC != NULL) {
00712          sprintf (stmp,"VP->VOLREG_MATVEC = \n\tMrot\tDelta\n");
00713          SS = SUMA_StringAppend (SS, stmp);
00714          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00715          VP->VOLREG_MATVEC[0], VP->VOLREG_MATVEC[1], VP->VOLREG_MATVEC[2], VP->VOLREG_MATVEC[3]); 
00716          SS = SUMA_StringAppend (SS, stmp);
00717          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00718          VP->VOLREG_MATVEC[4], VP->VOLREG_MATVEC[5], VP->VOLREG_MATVEC[6], VP->VOLREG_MATVEC[7]);
00719          SS = SUMA_StringAppend (SS, stmp);
00720          sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00721          VP->VOLREG_MATVEC[8], VP->VOLREG_MATVEC[9], VP->VOLREG_MATVEC[10], VP->VOLREG_MATVEC[11]);
00722          SS = SUMA_StringAppend (SS, stmp);
00723       } else {
00724          sprintf (stmp,"VP->VOLREG_MATVEC = NULL\n");
00725          SS = SUMA_StringAppend (SS, stmp);
00726       }
00727 
00728       if (VP->VOLREG_CENTER_OLD != NULL) {
00729          sprintf (stmp,"VP->VOLREG_CENTER_OLD = %f, %f, %f\n", \
00730            VP->VOLREG_CENTER_OLD[0], VP->VOLREG_CENTER_OLD[1], VP->VOLREG_CENTER_OLD[2]); 
00731          SS = SUMA_StringAppend (SS, stmp);
00732       }else {
00733          sprintf (stmp,"VP->VOLREG_CENTER_OLD = NULL\n");
00734          SS = SUMA_StringAppend (SS, stmp);
00735       }
00736 
00737       if (VP->VOLREG_CENTER_BASE != NULL) {
00738          sprintf (stmp,"VP->VOLREG_CENTER_BASE = %f, %f, %f\n", \
00739             VP->VOLREG_CENTER_BASE[0], VP->VOLREG_CENTER_BASE[1], VP->VOLREG_CENTER_BASE[2]); 
00740          SS = SUMA_StringAppend (SS, stmp);
00741       } else {
00742          sprintf (stmp,"VP->VOLREG_CENTER_BASE = NULL\n");
00743          SS = SUMA_StringAppend (SS, stmp);
00744       }
00745    }else{
00746       sprintf (stmp, "NULL Volume Parent Pointer.\n");
00747       SS = SUMA_StringAppend (SS, stmp);
00748    }
00749    
00750    /* clean SS */
00751    SS = SUMA_StringAppend (SS, NULL);
00752    /* copy s pointer and free SS */
00753    s = SS->s;
00754    SUMA_free(SS); 
00755    
00756    SUMA_RETURN (s);
00757 }
00758 /*!
00759    Show the contents of SUMA_VOLPAR structure 
00760    \sa SUMA_VolPar_Info 
00761 */
00762 void SUMA_Show_VolPar(SUMA_VOLPAR *VP, FILE *Out)
00763 {
00764    static char FuncName[]={"SUMA_Show_VolPar"};
00765    char *s;
00766    
00767    SUMA_ENTRY;
00768 
00769    if (Out == NULL) Out = SUMA_STDOUT;
00770    
00771    s =  SUMA_VolPar_Info(VP);
00772    
00773    if (s) {
00774       fprintf (Out, "%s", s);
00775       SUMA_free(s);
00776    }else {
00777       fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_VolPar_Info.\n", FuncName);
00778    }   
00779    
00780    SUMA_RETURNe;
00781       
00782 }
00783 
00784 /*!
00785    \brief Transform the coordinates of a surface object to AFNI-DICOM convention and transform the coordinates by SO->VolPar
00786    ans = SUMA_Align_to_VolPar (SO, SF_Struct);
00787    \param SO (SUMA_SurfaceObject *)
00788    \param S_Struct (void *) That is only needed for SureFit surfaces and is nothing but a type cast of a SUMA_SureFit_struct containing information on cropping.
00789                               send NULL for all other surfaces.
00790    \return YUP/NOPE
00791    For SureFit and FreeSurfer surfaces, the coordinates are first set in RAI (DICOM) coordinate system before applying SO->VolPar.
00792    For other surface formats, SO->VolPar is applied to whatever coordinates are in SO->NodeList
00793 */
00794 SUMA_Boolean SUMA_Align_to_VolPar (SUMA_SurfaceObject *SO, void * S_Struct)
00795 {
00796    static char FuncName[]={"SUMA_Align_to_VolPar"};
00797    char  orcode[3];
00798    float xx, yy, zz;
00799    THD_coorder * cord_surf, *cord_RAI;
00800    int i, ND, id;
00801    SUMA_SureFit_struct *SF;
00802    SUMA_FreeSurfer_struct *FS;
00803    
00804    SUMA_ENTRY;
00805 
00806    /* allocate for cord structure */
00807    cord_surf = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
00808    cord_RAI = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
00809    if (cord_surf == NULL || cord_RAI == NULL) {
00810       fprintf(SUMA_STDERR,"Error %s: failed to allocate for cord\n", FuncName);
00811       SUMA_RETURN (NOPE);
00812    }
00813    
00814    /* do the dance to get cord for RAI */
00815    THD_coorder_fill( NULL, cord_RAI);
00816    ND = SO->NodeDim;
00817    switch (SO->FileType) {
00818       case SUMA_INVENTOR_GENERIC:
00819       case SUMA_OPENDX_MESH:
00820       case SUMA_PLY:
00821       case SUMA_VEC:
00822          /* Do nothing */
00823          break;
00824       case SUMA_FREE_SURFER:
00825       case SUMA_FREE_SURFER_PATCH:
00826          /* For free surfer, all you need to do is 
00827           go from LPI to RAI (DICOM)*/
00828          sprintf(orcode,"LPI");
00829          THD_coorder_fill(orcode , cord_surf); 
00830          /*loop over XYZs and change them to dicom*/
00831          for (i=0; i < SO->N_Node; ++i) {
00832             id = i * ND;
00833             THD_coorder_to_dicom (cord_surf, &(SO->NodeList[id]), &(SO->NodeList[id+1]), &(SO->NodeList[id+2])); 
00834          }
00835          break;
00836       case SUMA_SUREFIT:
00837          /* For SureFit, coordinates are actually a float version of the indices */
00838          SF = (SUMA_SureFit_struct *)S_Struct;
00839          {   THD_fvec3 fv, iv;
00840             float D[3];
00841             /* Calcluate Delta caused by cropping */
00842             for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
00843             /* fprintf (SUMA_STDERR,"%s: Shift Values: [%f, %f, %f]\n", FuncName, D[0], D[1], D[2]); */
00844             for (i=0; i < SO->N_Node; ++i) {
00845                id = i * ND;
00846                /* change float indices to mm coords */
00847                iv.xyz[0] = SO->NodeList[id] + D[0];
00848                iv.xyz[1] = SO->NodeList[id+1] + D[1];
00849                iv.xyz[2] = SO->NodeList[id+2] + D[2];
00850                fv = SUMA_THD_3dfind_to_3dmm( SO, iv );
00851                
00852                /* change mm to RAI coords */
00853                iv = SUMA_THD_3dmm_to_dicomm( SO->VolPar->xxorient, SO->VolPar->yyorient, SO->VolPar->zzorient,  fv );
00854                SO->NodeList[id] = iv.xyz[0];
00855                SO->NodeList[id+1] = iv.xyz[1];
00856                SO->NodeList[id+2] = iv.xyz[2];
00857             }
00858          }
00859          break;
00860       case SUMA_BRAIN_VOYAGER:
00861          /* For Brain Voyager, all you need to do is 
00862           go from AIR to RAI (DICOM)
00863           Note: The center of the volume is at the 1st voxel's center and that huge
00864           center shift, relative to standard AFNI dsets (centered about middle of volume)
00865           might throw off 3dVolreg. If you want to shift volume's center to be in
00866           the middle voxel, you'll need to shift the surface coordinates before transforming
00867           them to ASR*/
00868          sprintf(orcode,"ASR");
00869          THD_coorder_fill(orcode , cord_surf); 
00870          /*loop over XYZs and change them to dicom*/
00871          for (i=0; i < SO->N_Node; ++i) {
00872             id = i * ND;
00873             THD_coorder_to_dicom (cord_surf, &(SO->NodeList[id]), &(SO->NodeList[id+1]), &(SO->NodeList[id+2])); 
00874          }
00875          break;
00876       default:
00877          fprintf(SUMA_STDERR,"Warning %s: Unknown SO->FileType. Assuming coordinates are in DICOM already.\n", FuncName);
00878          break;
00879    }
00880    
00881    if (!SUMA_Apply_VolReg_Trans (SO)) {
00882       fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Apply_VolReg_Trans.\n", FuncName);
00883       SUMA_RETURN (NOPE);
00884    }
00885 
00886    SUMA_RETURN (YUP);
00887 }
00888 
00889 /*!
00890 \brief applies the transform in SO->VolPar to SO->NodeList. This only makes sense if 
00891 SO->NodeList is in DICOM to begin with...
00892 \param SO (SUMA_SurfaceObject *) You better have NodeList and VolPar in there...
00893 \return YUP/NOPE
00894 \sa SUMA_Align_to_VolPar
00895 NOTE: The field SO->VOLREG_APPLIED is set to NOPE if this function fails, YUP otherwise
00896 */
00897 SUMA_Boolean SUMA_Apply_VolReg_Trans (SUMA_SurfaceObject *SO)
00898 {
00899    static char FuncName[]={"SUMA_Apply_VolReg_Trans"};
00900    int i, ND, id;
00901    SUMA_Boolean UseVolreg, UseTagAlign, UseRotate, Bad=YUP;
00902    
00903    SUMA_ENTRY;
00904 
00905    if (SUMAg_CF->IgnoreVolreg) {
00906       SUMA_SL_Note("Ignoring any Volreg or TagAlign transforms present in Surface Volume.\n");
00907       SO->VOLREG_APPLIED = NOPE;
00908       SUMA_RETURN (YUP);
00909    }
00910    
00911    if (SO->VOLREG_APPLIED || SO->TAGALIGN_APPLIED || SO->ROTATE_APPLIED) {
00912       fprintf (SUMA_STDERR,"Error %s: Volreg (or Tagalign or rotate) already applied. Nothing done.\n", FuncName);
00913       SUMA_RETURN (NOPE);
00914    }
00915 
00916    ND = SO->NodeDim;
00917    
00918    
00919    UseVolreg = NOPE;
00920    UseTagAlign = NOPE;
00921    UseRotate = NOPE;
00922    /* perform the rotation needed to align the surface with the current experiment's data */
00923    if (SO->VolPar->VOLREG_MATVEC != NULL || SO->VolPar->VOLREG_CENTER_OLD != NULL || SO->VolPar->VOLREG_CENTER_BASE != NULL) {
00924       Bad = NOPE;
00925       if (SO->VolPar->VOLREG_MATVEC == NULL) {
00926          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->VOLREG_MATVEC = NULL. Cannot perform alignment.\n", FuncName);
00927          Bad = YUP;
00928       }
00929       if (SO->VolPar->VOLREG_CENTER_OLD == NULL) {
00930          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->VOLREG_CENTER_OLD = NULL. Cannot perform alignment.\n", FuncName);
00931          Bad = YUP;
00932       }
00933       if (SO->VolPar->VOLREG_CENTER_BASE == NULL) {
00934          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->VOLREG_CENTER_BASE = NULL. Cannot perform alignment.\n", FuncName);
00935          Bad = YUP;
00936       }
00937       if (!Bad) UseVolreg = YUP;
00938    }
00939    
00940    /* Check for Tagalign and Rotate field */
00941    if (SO->VolPar->TAGALIGN_MATVEC) UseTagAlign = YUP;
00942    if (SO->VolPar->ROTATE_MATVEC || SO->VolPar->ROTATE_CENTER_OLD || SO->VolPar->ROTATE_CENTER_BASE) {
00943       Bad = NOPE;
00944       if (SO->VolPar->ROTATE_MATVEC == NULL) {
00945          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->ROTATE_MATVEC = NULL. Cannot perform alignment.\n", FuncName);
00946          Bad = YUP;
00947       }
00948       if (SO->VolPar->ROTATE_CENTER_OLD == NULL) {
00949          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->ROTATE_CENTER_OLD = NULL. Cannot perform alignment.\n", FuncName);
00950          Bad = YUP;
00951       }
00952       if (SO->VolPar->ROTATE_CENTER_BASE == NULL) {
00953          fprintf(SUMA_STDERR,"Error %s: SO->VolPar->ROTATE_CENTER_BASE = NULL. Cannot perform alignment.\n", FuncName);
00954          Bad = YUP;
00955       }
00956       if (!Bad) UseRotate = YUP;
00957    }
00958    
00959    if (UseTagAlign && UseVolreg) {
00960       SUMA_SL_Note("Found both Volreg and TagAlign fields.\n"
00961                    "Using Volreg fields for alignment.");
00962       UseTagAlign = NOPE;
00963    }
00964    if (UseRotate && UseVolreg) {
00965       SUMA_SL_Note("Found both Volreg and Rotate fields.\n"
00966                    "Using Volreg fields for alignment.");
00967       UseRotate = NOPE;
00968    }
00969    if (UseRotate && UseTagAlign) {
00970       SUMA_SL_Note("Found both Rotate and TagAlign fields.\n"
00971                    "Using Tagalign fields for alignment.");
00972       UseRotate = NOPE;
00973    }
00974    
00975    if (UseTagAlign) {
00976       float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
00977       
00978       /* fillerup*/
00979       Mrot[0][0] = SO->VolPar->TAGALIGN_MATVEC[0];
00980       Mrot[0][1] = SO->VolPar->TAGALIGN_MATVEC[1];
00981       Mrot[0][2] = SO->VolPar->TAGALIGN_MATVEC[2];
00982       Delta[0]   = SO->VolPar->TAGALIGN_MATVEC[3];
00983       Mrot[1][0] = SO->VolPar->TAGALIGN_MATVEC[4];
00984       Mrot[1][1] = SO->VolPar->TAGALIGN_MATVEC[5];
00985       Mrot[1][2] = SO->VolPar->TAGALIGN_MATVEC[6];
00986       Delta[1]   = SO->VolPar->TAGALIGN_MATVEC[7];   
00987       Mrot[2][0] = SO->VolPar->TAGALIGN_MATVEC[8];
00988       Mrot[2][1] = SO->VolPar->TAGALIGN_MATVEC[9];
00989       Mrot[2][2] = SO->VolPar->TAGALIGN_MATVEC[10];
00990       Delta[2]   = SO->VolPar->TAGALIGN_MATVEC[11];
00991       
00992       NetShift[0] = Delta[0];
00993       NetShift[1] = Delta[1];
00994       NetShift[2] = Delta[2];
00995       
00996       /*
00997       fprintf (SUMA_STDERR,"%s: Applying Rotation.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
00998                Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
00999                Delta[0], Delta[1], Delta[2]);
01000       
01001       */
01002       
01003       for (i=0; i < SO->N_Node; ++i) {
01004          id = ND * i;
01005          /* zero the center */ 
01006          x = SO->NodeList[id] ;
01007          y = SO->NodeList[id+1] ;
01008          z = SO->NodeList[id+2] ;
01009          
01010          /* Apply the rotation matrix XYZn = Mrot x XYZ*/
01011          SO->NodeList[id] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
01012          SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
01013          SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
01014          
01015          /*apply netshift*/
01016          SO->NodeList[id] += NetShift[0];
01017          SO->NodeList[id+1] += NetShift[1];
01018          SO->NodeList[id+2] += NetShift[2];
01019       }
01020       SO->TAGALIGN_APPLIED = YUP;   
01021    } else
01022       SO->TAGALIGN_APPLIED = NOPE;
01023          
01024    if (UseVolreg) {
01025       float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
01026       
01027       /* fillerup*/
01028       Mrot[0][0] = SO->VolPar->VOLREG_MATVEC[0];
01029       Mrot[0][1] = SO->VolPar->VOLREG_MATVEC[1];
01030       Mrot[0][2] = SO->VolPar->VOLREG_MATVEC[2];
01031       Delta[0]   = SO->VolPar->VOLREG_MATVEC[3];
01032       Mrot[1][0] = SO->VolPar->VOLREG_MATVEC[4];
01033       Mrot[1][1] = SO->VolPar->VOLREG_MATVEC[5];
01034       Mrot[1][2] = SO->VolPar->VOLREG_MATVEC[6];
01035       Delta[1]   = SO->VolPar->VOLREG_MATVEC[7];   
01036       Mrot[2][0] = SO->VolPar->VOLREG_MATVEC[8];
01037       Mrot[2][1] = SO->VolPar->VOLREG_MATVEC[9];
01038       Mrot[2][2] = SO->VolPar->VOLREG_MATVEC[10];
01039       Delta[2]   = SO->VolPar->VOLREG_MATVEC[11];
01040       
01041       NetShift[0] = SO->VolPar->VOLREG_CENTER_BASE[0] + Delta[0];
01042       NetShift[1] = SO->VolPar->VOLREG_CENTER_BASE[1] + Delta[1];
01043       NetShift[2] = SO->VolPar->VOLREG_CENTER_BASE[2] + Delta[2];
01044       
01045       /*
01046       fprintf (SUMA_STDERR,"%s: Applying Rotation.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
01047                Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
01048                Delta[0], Delta[1], Delta[2]);
01049       fprintf (SUMA_STDERR,"VOLREG_CENTER_BASE = [%f %f %f]. VOLREG_CENTER_OLD = [%f %f %f]\n", \
01050          SO->VolPar->VOLREG_CENTER_BASE[0], SO->VolPar->VOLREG_CENTER_BASE[1], SO->VolPar->VOLREG_CENTER_BASE[2], \
01051          SO->VolPar->VOLREG_CENTER_OLD[0], SO->VolPar->VOLREG_CENTER_OLD[1], SO->VolPar->VOLREG_CENTER_OLD[2]);
01052       */
01053       
01054       for (i=0; i < SO->N_Node; ++i) {
01055          id = ND * i;
01056          /* zero the center */ 
01057          x = SO->NodeList[id  ] - SO->VolPar->VOLREG_CENTER_OLD[0];
01058          y = SO->NodeList[id+1] - SO->VolPar->VOLREG_CENTER_OLD[1];
01059          z = SO->NodeList[id+2] - SO->VolPar->VOLREG_CENTER_OLD[2];
01060          
01061          /* Apply the rotation matrix XYZn = Mrot x XYZ*/
01062          SO->NodeList[id  ] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
01063          SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
01064          SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
01065          
01066          /*apply netshift*/
01067          SO->NodeList[id  ] += NetShift[0];
01068          SO->NodeList[id+1] += NetShift[1];
01069          SO->NodeList[id+2] += NetShift[2];
01070       }
01071       SO->VOLREG_APPLIED = YUP;   
01072    } else
01073       SO->VOLREG_APPLIED = NOPE;
01074    
01075    if (UseRotate) {
01076       float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
01077       
01078       /* fillerup*/
01079       Mrot[0][0] = SO->VolPar->ROTATE_MATVEC[0];
01080       Mrot[0][1] = SO->VolPar->ROTATE_MATVEC[1];
01081       Mrot[0][2] = SO->VolPar->ROTATE_MATVEC[2];
01082       Delta[0]   = SO->VolPar->ROTATE_MATVEC[3];
01083       Mrot[1][0] = SO->VolPar->ROTATE_MATVEC[4];
01084       Mrot[1][1] = SO->VolPar->ROTATE_MATVEC[5];
01085       Mrot[1][2] = SO->VolPar->ROTATE_MATVEC[6];
01086       Delta[1]   = SO->VolPar->ROTATE_MATVEC[7];   
01087       Mrot[2][0] = SO->VolPar->ROTATE_MATVEC[8];
01088       Mrot[2][1] = SO->VolPar->ROTATE_MATVEC[9];
01089       Mrot[2][2] = SO->VolPar->ROTATE_MATVEC[10];
01090       Delta[2]   = SO->VolPar->ROTATE_MATVEC[11];
01091       
01092       NetShift[0] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[0];
01093       NetShift[1] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[1];
01094       NetShift[2] = SO->VolPar->ROTATE_CENTER_BASE[0] + Delta[2];
01095       
01096       /*
01097       fprintf (SUMA_STDERR,"%s: Applying Rotation.\nMrot[\t%f\t%f\t%f\n%f\t%f\t%f\n%f\t%f\t%f]\nDelta = [%f %f %f]\n", FuncName,\
01098                Mrot[0][0], Mrot[0][1], Mrot[0][2], Mrot[1][0], Mrot[1][1], Mrot[1][2], Mrot[2][0], Mrot[2][1], Mrot[2][2], \
01099                Delta[0], Delta[1], Delta[2]);
01100       fprintf (SUMA_STDERR,"ROTATE_CENTER_BASE = [%f %f %f]. ROTATE_CENTER_OLD = [%f %f %f]\n", \
01101          SO->VolPar->ROTATE_CENTER_BASE[0], SO->VolPar->ROTATE_CENTER_BASE[1], SO->VolPar->ROTATE_CENTER_BASE[2], \
01102          SO->VolPar->ROTATE_CENTER_OLD[0], SO->VolPar->ROTATE_CENTER_OLD[1], SO->VolPar->ROTATE_CENTER_OLD[2]);
01103       */
01104       
01105       for (i=0; i < SO->N_Node; ++i) {
01106          id = ND * i;
01107          /* zero the center */ 
01108          x = SO->NodeList[id  ] - SO->VolPar->ROTATE_CENTER_OLD[0];
01109          y = SO->NodeList[id+1] - SO->VolPar->ROTATE_CENTER_OLD[1];
01110          z = SO->NodeList[id+2] - SO->VolPar->ROTATE_CENTER_OLD[1];
01111          
01112          /* Apply the rotation matrix XYZn = Mrot x XYZ*/
01113          SO->NodeList[id  ] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
01114          SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
01115          SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
01116          
01117          /*apply netshift*/
01118          SO->NodeList[id  ] += NetShift[0];
01119          SO->NodeList[id+1] += NetShift[1];
01120          SO->NodeList[id+2] += NetShift[2];
01121       }
01122       SO->ROTATE_APPLIED = YUP;   
01123    } else
01124       SO->ROTATE_APPLIED = NOPE;
01125    
01126    SUMA_RETURN (YUP);
01127 }
01128 
01129 /*! The following functions are adaptations from thd_coords.c 
01130 They are modified to avoid the use of dset data structure 
01131 which is not fully preserved in the Surface Object Data structure.
01132 
01133 Stolen Comment:
01134 ====================================================================
01135    3D coordinate conversion routines;
01136      tags for coordinate systems:
01137        fdind  = FD_brick voxel indices (ints)
01138        fdfind = FD_brick voxel indices (floats - added 30 Aug 2001)
01139        3dind  = THD_3dim_dataset voxel indices (ints)
01140        3dfind = THD_3dim_dataset floating point voxel indices
01141                   (used for subvoxel resolution)
01142        3dmm   = THD_3dim_dataset millimetric coordinates (floats)
01143        dicomm = DICOM 3.0 millimetric coordinates (floats)
01144 
01145      The 3dmm coordinate measurements are oriented the same way
01146      as the dicomm coordinates (which means that the ??del values
01147      can be negative if the voxel axes are reflected), but the 3dmm
01148      coordinates are in general a permutation of the DICOM coordinates.
01149 
01150      These routines all take as input and produce as output
01151      THD_?vec3 (i=int,f=float) structs.
01152 ======================================================================
01153 
01154 */
01155 
01156 THD_fvec3 SUMA_THD_3dfind_to_3dmm( SUMA_SurfaceObject *SO, 
01157                               THD_fvec3 iv )
01158 {
01159    static char FuncName[]={"SUMA_THD_3dfind_to_3dmm"};
01160    THD_fvec3     fv ;
01161 
01162    SUMA_ENTRY;
01163 
01164    fv.xyz[0] = SO->VolPar->xorg + iv.xyz[0] * SO->VolPar->dx ;
01165    fv.xyz[1] = SO->VolPar->yorg + iv.xyz[1] * SO->VolPar->dy ;
01166    fv.xyz[2] = SO->VolPar->zorg + iv.xyz[2] * SO->VolPar->dz ;
01167    SUMA_RETURN(fv) ;
01168 }
01169 
01170 /*!------------------------------------------------------------------*/
01171 
01172 THD_fvec3 SUMA_THD_3dind_to_3dmm( SUMA_SurfaceObject *SO,
01173                                    THD_ivec3 iv )
01174 {
01175    static char FuncName[]={"SUMA_THD_3dind_to_3dmm"};
01176    THD_fvec3     fv ;
01177 
01178    SUMA_ENTRY;
01179 
01180    fv.xyz[0] = SO->VolPar->xorg + iv.ijk[0] * SO->VolPar->dx ;
01181    fv.xyz[1] = SO->VolPar->yorg + iv.ijk[1] * SO->VolPar->dy ;
01182    fv.xyz[2] = SO->VolPar->zorg + iv.ijk[2] * SO->VolPar->dz ;
01183    SUMA_RETURN(fv) ;
01184 }
01185 
01186 /*!--------------------------------------------------------------------*/
01187 
01188 THD_fvec3 SUMA_THD_3dmm_to_3dfind( SUMA_SurfaceObject *SO ,
01189                               THD_fvec3 fv )
01190 {
01191    static char FuncName[]={"SUMA_THD_3dmm_to_3dfind"};
01192    THD_fvec3     iv ;
01193 
01194    SUMA_ENTRY;
01195 
01196 
01197    iv.xyz[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx ;
01198    iv.xyz[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy ;
01199    iv.xyz[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz ;
01200 
01201         if( iv.xyz[0] < 0            ) iv.xyz[0] = 0 ;
01202    else if( iv.xyz[0] > SO->VolPar->nx-1 ) iv.xyz[0] = SO->VolPar->nx-1 ;
01203 
01204         if( iv.xyz[1] <  0           ) iv.xyz[1] = 0 ;
01205    else if( iv.xyz[1] > SO->VolPar->ny-1 ) iv.xyz[1] = SO->VolPar->ny-1 ;
01206 
01207         if( iv.xyz[2] < 0            ) iv.xyz[2] = 0 ;
01208    else if( iv.xyz[2] > SO->VolPar->nz-1 ) iv.xyz[2] = SO->VolPar->nz-1 ;
01209 
01210    SUMA_RETURN(iv) ;
01211 }
01212 
01213 /*!--------------------------------------------------------------------*/
01214 
01215 THD_ivec3 SUMA_THD_3dmm_to_3dind( SUMA_SurfaceObject *SO  ,
01216                              THD_fvec3 fv )
01217 {
01218    static char FuncName[]={"SUMA_THD_3dmm_to_3dind"};
01219    THD_ivec3     iv ;
01220 
01221    SUMA_ENTRY;
01222 
01223    iv.ijk[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx + 0.499 ;
01224    iv.ijk[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy + 0.499 ;
01225    iv.ijk[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz + 0.499 ;
01226 
01227         if( iv.ijk[0] < 0            ) iv.ijk[0] = 0 ;
01228    else if( iv.ijk[0] > SO->VolPar->nx-1 ) iv.ijk[0] = SO->VolPar->nx-1 ;
01229 
01230         if( iv.ijk[1] < 0            ) iv.ijk[1] = 0 ;
01231    else if( iv.ijk[1] > SO->VolPar->ny-1 ) iv.ijk[1] = SO->VolPar->ny-1 ;
01232 
01233         if( iv.ijk[2] < 0            ) iv.ijk[2] = 0 ;
01234    else if( iv.ijk[2] > SO->VolPar->nz-1 ) iv.ijk[2] = SO->VolPar->nz-1 ;
01235 
01236    SUMA_RETURN(iv) ;
01237 }
01238 
01239 /*! 
01240    \brief how many voxels in each of the RL AP IS directions
01241 */
01242 void SUMA_VolDims(THD_3dim_dataset *dset, int *nRL, int *nAP, int *nIS)
01243 {
01244    static char FuncName[]={"SUMA_VolDims"};
01245    
01246    SUMA_ENTRY;
01247    
01248    *nRL = *nAP = *nIS = -1;
01249    
01250    if (!dset) {
01251       SUMA_SL_Err("NULL dset");
01252       SUMA_RETURNe;
01253    }
01254    
01255    switch( dset->daxes->xxorient ){
01256       case ORI_R2L_TYPE:
01257       case ORI_L2R_TYPE: *nRL = DSET_NX(dset) ; break ;
01258       case ORI_P2A_TYPE:
01259       case ORI_A2P_TYPE: *nAP = DSET_NX(dset) ; break ;
01260       case ORI_I2S_TYPE:
01261       case ORI_S2I_TYPE: *nIS = DSET_NX(dset) ; break ;
01262    }
01263 
01264    switch( dset->daxes->yyorient ){
01265       case ORI_R2L_TYPE:
01266       case ORI_L2R_TYPE: *nRL = DSET_NY(dset) ; break ;
01267       case ORI_P2A_TYPE:
01268       case ORI_A2P_TYPE: *nAP = DSET_NY(dset) ; break ;
01269       case ORI_I2S_TYPE:
01270       case ORI_S2I_TYPE: *nIS = DSET_NY(dset) ; break ;
01271    }
01272 
01273    switch( dset->daxes->zzorient ){
01274       case ORI_R2L_TYPE:
01275       case ORI_L2R_TYPE: *nRL = DSET_NZ(dset) ; break ;
01276       case ORI_P2A_TYPE:
01277       case ORI_A2P_TYPE: *nAP = DSET_NZ(dset) ; break ;
01278       case ORI_I2S_TYPE:
01279       case ORI_S2I_TYPE: *nIS = DSET_NZ(dset) ; break ;
01280    }
01281    
01282    SUMA_RETURNe;
01283 }
01284 
01285 /*!---------------------------------------------------------------------
01286    convert from input image oriented x,y,z to Dicom x,y,z
01287      (x axis = R->L , y axis = A->P , z axis = I->S)
01288 
01289    N.B.: image distances are oriented the same as Dicom,
01290          just in a permuted order.
01291 -----------------------------------------------------------------------*/
01292 THD_fvec3 SUMA_THD_3dmm_to_dicomm( int xxorient, int yyorient, int zzorient, 
01293                               THD_fvec3 imv )
01294 {
01295    static char FuncName[]={"SUMA_THD_3dmm_to_dicomm"};   
01296    THD_fvec3 dicv ;
01297    float xim,yim,zim , xdic = 0.0, ydic = 0.0, zdic = 0.0 ;
01298 
01299    SUMA_ENTRY;
01300 
01301    xim = imv.xyz[0] ; yim = imv.xyz[1] ; zim = imv.xyz[2] ;
01302 
01303    switch( xxorient ){
01304       case ORI_R2L_TYPE:
01305       case ORI_L2R_TYPE: xdic = xim ; break ;
01306       case ORI_P2A_TYPE:
01307       case ORI_A2P_TYPE: ydic = xim ; break ;
01308       case ORI_I2S_TYPE:
01309       case ORI_S2I_TYPE: zdic = xim ; break ;
01310 
01311       default: 
01312          fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
01313          exit (1);
01314    }
01315 
01316    switch( yyorient ){
01317       case ORI_R2L_TYPE:
01318       case ORI_L2R_TYPE: xdic = yim ; break ;
01319       case ORI_P2A_TYPE:
01320       case ORI_A2P_TYPE: ydic = yim ; break ;
01321       case ORI_I2S_TYPE:
01322       case ORI_S2I_TYPE: zdic = yim ; break ;
01323 
01324       default: 
01325          fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
01326          exit (1);
01327    }
01328 
01329    switch( zzorient ){
01330       case ORI_R2L_TYPE:
01331       case ORI_L2R_TYPE: xdic = zim ; break ;
01332       case ORI_P2A_TYPE:
01333       case ORI_A2P_TYPE: ydic = zim ; break ;
01334       case ORI_I2S_TYPE:
01335       case ORI_S2I_TYPE: zdic = zim ; break ;
01336 
01337       default: 
01338          fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
01339          exit (1);
01340    }
01341 
01342    dicv.xyz[0] = xdic ; dicv.xyz[1] = ydic ; dicv.xyz[2] = zdic ;
01343    SUMA_RETURN(dicv) ;
01344 }
01345 
01346 /*!---------------------------------------------------------------------
01347    convert to input image oriented x,y,z from Dicom x,y,z
01348 -----------------------------------------------------------------------*/
01349 
01350 THD_fvec3 SUMA_THD_dicomm_to_3dmm( SUMA_SurfaceObject *SO ,
01351                               THD_fvec3 dicv )
01352 {
01353    static char FuncName[]={"SUMA_THD_dicomm_to_3dmm"};
01354    THD_fvec3 imv ;
01355    float xim,yim,zim , xdic,ydic,zdic ;
01356 
01357    SUMA_ENTRY;
01358 
01359    xdic = dicv.xyz[0] ; ydic = dicv.xyz[1] ; zdic = dicv.xyz[2] ;
01360 
01361    switch( SO->VolPar->xxorient ){
01362       case ORI_R2L_TYPE:
01363       case ORI_L2R_TYPE: xim = xdic ; break ;
01364       case ORI_P2A_TYPE:
01365       case ORI_A2P_TYPE: xim = ydic ; break ;
01366       case ORI_I2S_TYPE:
01367       case ORI_S2I_TYPE: xim = zdic ; break ;
01368 
01369       default: 
01370          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
01371          exit (1);
01372    }
01373 
01374    switch( SO->VolPar->yyorient ){
01375       case ORI_R2L_TYPE:
01376       case ORI_L2R_TYPE: yim = xdic ; break ;
01377       case ORI_P2A_TYPE:
01378       case ORI_A2P_TYPE: yim = ydic ; break ;
01379       case ORI_I2S_TYPE:
01380       case ORI_S2I_TYPE: yim = zdic ; break ;
01381 
01382       default: 
01383          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
01384          exit (1);
01385    }
01386 
01387    switch( SO->VolPar->zzorient ){
01388       case ORI_R2L_TYPE:
01389       case ORI_L2R_TYPE: zim = xdic ; break ;
01390       case ORI_P2A_TYPE:
01391       case ORI_A2P_TYPE: zim = ydic ; break ;
01392       case ORI_I2S_TYPE:
01393       case ORI_S2I_TYPE: zim = zdic ; break ;
01394 
01395       default: 
01396          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n Exiting.") ;
01397          exit (1);
01398    }
01399 
01400    imv.xyz[0] = xim ; imv.xyz[1] = yim ; imv.xyz[2] = zim ;
01401    SUMA_RETURN(imv) ;
01402 }
01403 /*!
01404    \brief changes orientation codes to string
01405    classic RAI codes: ORI_R2L_TYPE, ORI_A2P_TYPE, ORI_I2S_TYPE
01406    would result in orstr being set to: RAILPS 
01407 */
01408 void SUMA_orcode_to_orstring (int xxorient,int  yyorient,int zzorient, char *orstr)
01409 {
01410    static char FuncName[]={"SUMA_orcode_to_orstring"};
01411    
01412    SUMA_ENTRY;
01413    
01414    if (!orstr) { SUMA_SL_Err("NULL string"); SUMA_RETURNe; }
01415    
01416    orstr[0]='\0';
01417    switch( xxorient ){
01418       case ORI_R2L_TYPE: orstr[0] = 'R'; orstr[3] = 'L'; break ;
01419       case ORI_L2R_TYPE: orstr[0] = 'L'; orstr[3] = 'R'; break ;
01420       case ORI_P2A_TYPE: orstr[0] = 'P'; orstr[3] = 'A'; break ;
01421       case ORI_A2P_TYPE: orstr[0] = 'A'; orstr[3] = 'P'; break ;
01422       case ORI_I2S_TYPE: orstr[0] = 'I'; orstr[3] = 'S'; break ;
01423       case ORI_S2I_TYPE: orstr[0] = 'S'; orstr[3] = 'I'; break ;
01424 
01425       default: 
01426          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal xxorient code.\n") ;
01427          SUMA_RETURNe;
01428    }
01429 
01430    switch( yyorient ){
01431       case ORI_R2L_TYPE: orstr[1] = 'R'; orstr[4] = 'L'; break ;
01432       case ORI_L2R_TYPE: orstr[1] = 'L'; orstr[4] = 'R'; break ;
01433       case ORI_P2A_TYPE: orstr[1] = 'P'; orstr[4] = 'A'; break ;
01434       case ORI_A2P_TYPE: orstr[1] = 'A'; orstr[4] = 'P'; break ;
01435       case ORI_I2S_TYPE: orstr[1] = 'I'; orstr[4] = 'S'; break ;
01436       case ORI_S2I_TYPE: orstr[1] = 'S'; orstr[4] = 'I'; break ;
01437 
01438       default: 
01439          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal yyorient code.\n ") ;
01440          SUMA_RETURNe;
01441    }
01442 
01443    switch( zzorient ){
01444       case ORI_R2L_TYPE: orstr[2] = 'R'; orstr[5] = 'L'; break ;
01445       case ORI_L2R_TYPE: orstr[2] = 'L'; orstr[5] = 'R'; break ;
01446       case ORI_P2A_TYPE: orstr[2] = 'P'; orstr[5] = 'A'; break ;
01447       case ORI_A2P_TYPE: orstr[2] = 'A'; orstr[5] = 'P'; break ;
01448       case ORI_I2S_TYPE: orstr[2] = 'I'; orstr[5] = 'S'; break ;
01449       case ORI_S2I_TYPE: orstr[2] = 'S'; orstr[5] = 'I'; break ;
01450 
01451       default: 
01452          fprintf(SUMA_STDERR, "SUMA_THD_dicomm_to_3dmm: illegal zzorient code.\n ") ;
01453          SUMA_RETURNe;
01454    }
01455    
01456    SUMA_RETURNe;
01457 }
01458 SUMA_Boolean SUMA_orstring_to_orcode (char *orstr, int *orient)
01459 {
01460    static char FuncName[]={"SUMA_orstring_to_orcode"};
01461    int i;
01462    
01463    SUMA_ENTRY;
01464    
01465    if (!orstr) { SUMA_SL_Err("NULL string"); SUMA_RETURN(NOPE); }
01466    if (!SUMA_ok_orstring(orstr)) { SUMA_SL_Err("Bad orientation string"); SUMA_RETURN(NOPE); }
01467    for (i=0; i<3; ++i) {
01468       switch (orstr[i]) {
01469          case 'R': orient[i] = ORI_R2L_TYPE; break;
01470          case 'L': orient[i] = ORI_L2R_TYPE; break;
01471          case 'A': orient[i] = ORI_A2P_TYPE; break;
01472          case 'P': orient[i] = ORI_P2A_TYPE; break;
01473          case 'I': orient[i] = ORI_I2S_TYPE; break;
01474          case 'S': orient[i] = ORI_S2I_TYPE; break;
01475          default: fprintf(SUMA_STDERR, " SUMA_orstring_to_orcode: Bad to the bones\n"); SUMA_RETURN(NOPE); 
01476       }
01477    }
01478    
01479    SUMA_RETURN(YUP);
01480 }
01481 
01482 int SUMA_ok_orstring (char *orstr)
01483 {
01484    static char FuncName[]={"SUMA_ok_orstring"};
01485    int i, d[3];
01486    
01487    SUMA_ENTRY;
01488    
01489    if (!orstr) { SUMA_RETURN(NOPE); }
01490    d[0] = d[1] = d[2] = 0;
01491    for (i=0; i<3; ++i) {
01492       switch (orstr[i]) {
01493          case 'R': ++(d[0]); break;
01494          case 'L': ++(d[0]); break;
01495          case 'A': ++(d[1]); break;
01496          case 'P': ++(d[1]); break;
01497          case 'I': ++(d[2]); break;
01498          case 'S': ++(d[2]); break;
01499          default: SUMA_RETURN(NOPE); 
01500       }
01501    }
01502    if (d[0] != 1 || d[1] != 1 || d[2] != 1) SUMA_RETURN(NOPE);
01503     
01504    SUMA_RETURN(YUP);
01505 }
01506 int SUMA_flip_orient(int xxorient)
01507 {
01508    static char FuncName[]={"SUMA_flip_orient"};
01509    
01510    SUMA_ENTRY;
01511    
01512    switch( xxorient ){
01513       case ORI_R2L_TYPE: SUMA_RETURN(ORI_L2R_TYPE); break ;
01514       case ORI_L2R_TYPE: SUMA_RETURN(ORI_R2L_TYPE); break ;
01515       
01516       case ORI_P2A_TYPE: SUMA_RETURN(ORI_A2P_TYPE); break ;
01517       case ORI_A2P_TYPE: SUMA_RETURN(ORI_P2A_TYPE); break ;
01518       
01519       case ORI_I2S_TYPE: SUMA_RETURN(ORI_S2I_TYPE); break ;
01520       case ORI_S2I_TYPE: SUMA_RETURN(ORI_I2S_TYPE); break ;
01521 
01522       default: 
01523          fprintf(SUMA_STDERR, "SUMA_opposite_orient: illegal zzorient code.\n ") ;
01524          SUMA_RETURN(-1);
01525    }
01526    
01527 }
01528 
01529 SUMA_Boolean SUMA_CoordChange (char *orc_in, char *orc_out, float *XYZ, int N_xyz)
01530 {
01531    static char FuncName[]={"SUMA_CoordChange"};
01532    int i, or_in[3], or_out[3], j, map[3], sgn[3], i3;
01533    float xyz[3];
01534    
01535    SUMA_ENTRY;
01536    
01537    if (!SUMA_orstring_to_orcode(orc_in, or_in)) {
01538       SUMA_SL_Err("Bad in code");
01539       SUMA_RETURN(NOPE);
01540    }
01541    if (!SUMA_orstring_to_orcode(orc_out, or_out)) {
01542       SUMA_SL_Err("Bad out code");
01543       SUMA_RETURN(NOPE);
01544    }
01545    
01546    /* figure out the mapping */
01547    for (j=0; j<3; ++j) { 
01548       i = 0;
01549       while (i<3) {
01550          if (or_in[i] == or_out[j] || or_in[i] == SUMA_flip_orient(or_out[j])) {
01551             map[j] = i;
01552             if (or_in[i] == SUMA_flip_orient(or_out[j])) sgn[j] = -1;
01553             else sgn[j] = 1;
01554             i = 3; /* break */
01555          }
01556          ++i;
01557       }
01558    }
01559 
01560    for (i=0; i<N_xyz; ++i) {
01561       i3 = 3*i;
01562       xyz[0] = XYZ[i3]; xyz[1] = XYZ[i3+1]; xyz[2] = XYZ[i3+2];
01563       XYZ[i3  ] = sgn[0] * xyz[map[0]]; 
01564       XYZ[i3+1] = sgn[1] * xyz[map[1]]; 
01565       XYZ[i3+2] = sgn[2] * xyz[map[2]]; 
01566    }
01567    
01568    SUMA_RETURN(YUP);
01569 }
01570 /*!
01571    \brief takes the origin as entered to to3d and 
01572    changes them to the origin field for AFNI header
01573    Based on info in README.attributes and function T3D_save_file_CB
01574    int to3d.c
01575 */ 
01576 void SUMA_originto3d_2_originHEAD(THD_ivec3 orient, THD_fvec3 *origin)
01577 {
01578    static char FuncName[]={"SUMA_originto3d_2_originHEAD"};
01579    
01580    SUMA_ENTRY;
01581    
01582    origin->xyz[0] = (ORIENT_sign[orient.ijk[0]] == '+')
01583                   ? (-origin->xyz[0]) : ( origin->xyz[0]) ;
01584 
01585    origin->xyz[1] = (ORIENT_sign[orient.ijk[1]] == '+')
01586                   ? (-origin->xyz[1]) : ( origin->xyz[1]) ;
01587    
01588    origin->xyz[2] = (ORIENT_sign[orient.ijk[2]] == '+')
01589                   ? (-origin->xyz[2]) : ( origin->xyz[2]) ;
01590    
01591    SUMA_RETURNe;
01592    
01593 } 
01594 /*!
01595    \brief takes the size as entered to to3d and 
01596    changes them to the delta field for AFNI header
01597    Based on info in README.attributes and function T3D_save_file_CB
01598    int to3d.c
01599 */ 
01600 
01601 void SUMA_sizeto3d_2_deltaHEAD(THD_ivec3 orient, THD_fvec3 *delta)                 
01602 {
01603    static char FuncName[]={"SUMA_sizeto3d_2_deltaHEAD"};
01604    
01605    SUMA_ENTRY;
01606    
01607    delta->xyz[0] = (ORIENT_sign[orient.ijk[0]] == '+')
01608                   ? (delta->xyz[0]) : ( -delta->xyz[0]) ;
01609 
01610    delta->xyz[1] = (ORIENT_sign[orient.ijk[1]] == '+')
01611                   ? (delta->xyz[1]) : ( -delta->xyz[1]) ;
01612    
01613    delta->xyz[2] = (ORIENT_sign[orient.ijk[2]] == '+')
01614                   ? (delta->xyz[2]) : ( -delta->xyz[2]) ;
01615    
01616    SUMA_RETURNe;
01617    
01618 }    
01619 
01620 /*!
01621    \brief SUMA_Boolean SUMA_vec_3dfind_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01622    SUMA_Boolean SUMA_vec_3dmm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
01623    SUMA_Boolean SUMA_vec_dicomm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
01624    SUMA_Boolean SUMA_vec_3dfind_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
01625    SUMA_Boolean SUMA_vec_3dmm_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
01626    SUMA_Boolean SUMA_vec_dicomm_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar);
01627    a set of functions to change coordinate systems.
01628    
01629    \param NodeList (float *) vector of consecutive xyz (N_Node x 3) values
01630    \param N_Node (int) number of xyz triplets in NodeList
01631    \param VolPar (SUMA_VOLPAR *) volume parent structure (output of SUMA_VolPar_Attr )
01632    
01633    3dfind = float voxel/node index coordinate
01634    3dmm   = float voxel/node coordinate in volume space
01635    dicomm = float voxel/node coordinate in RAI space
01636    
01637    see also SUMA_THD_ functions ...
01638 */ 
01639 SUMA_Boolean SUMA_vec_3dfind_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01640 {
01641    static char FuncName[]={"SUMA_vec_3dfind_to_3dmm"};
01642    THD_fvec3 fv, iv;
01643    int i, id;
01644    SUMA_SurfaceObject SO;
01645    
01646    SUMA_ENTRY;
01647 
01648    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01649    /* create dummy struct */
01650    SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
01651    
01652    for (i=0; i < SO.N_Node; ++i) {
01653       id = i * SO.NodeDim;
01654       /* change float indices to mm coords */
01655       iv.xyz[0] = SO.NodeList[id] ;
01656       iv.xyz[1] = SO.NodeList[id+1] ;
01657       iv.xyz[2] = SO.NodeList[id+2] ;
01658       fv = SUMA_THD_3dfind_to_3dmm( &SO, iv );
01659       SO.NodeList[id] = fv.xyz[0];
01660       SO.NodeList[id+1] = fv.xyz[1];
01661       SO.NodeList[id+2] = fv.xyz[2];
01662    }
01663 
01664    SUMA_RETURN(YUP);
01665 }
01666 
01667 SUMA_Boolean SUMA_vec_3dmm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01668 {
01669    static char FuncName[]={"SUMA_vec_3dmm_to_3dfind"};
01670    THD_fvec3 fv, iv;
01671    int i, id;
01672    SUMA_SurfaceObject SO;
01673 
01674    SUMA_ENTRY;
01675 
01676    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01677    /* create dummy struct */
01678    SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
01679 
01680    for (i=0; i < SO.N_Node; ++i) {
01681       id = i * SO.NodeDim;
01682       /* change float indices to mm coords */
01683       iv.xyz[0] = SO.NodeList[id] ;
01684       iv.xyz[1] = SO.NodeList[id+1] ;
01685       iv.xyz[2] = SO.NodeList[id+2] ;
01686       fv = SUMA_THD_3dmm_to_3dfind( &SO, iv );
01687       SO.NodeList[id] = fv.xyz[0];
01688       SO.NodeList[id+1] = fv.xyz[1];
01689       SO.NodeList[id+2] = fv.xyz[2];
01690    }
01691 
01692    SUMA_RETURN(YUP);
01693 }
01694 
01695 SUMA_Boolean SUMA_vec_dicomm_to_3dfind (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01696 {
01697    static char FuncName[]={"SUMA_vec_dicomm_to_3dfind"};
01698 
01699    SUMA_ENTRY;
01700 
01701    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01702    
01703    
01704    if (!SUMA_vec_dicomm_to_3dmm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
01705    if (!SUMA_vec_3dmm_to_3dfind(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
01706 
01707    SUMA_RETURN(YUP);
01708 }
01709 
01710 SUMA_Boolean SUMA_vec_3dfind_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01711 {
01712    static char FuncName[]={"SUMA_vec_3dfind_to_dicomm"};
01713 
01714    SUMA_ENTRY;
01715 
01716    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01717 
01718    if (!SUMA_vec_3dfind_to_3dmm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
01719    if (!SUMA_vec_3dmm_to_dicomm(NodeList, N_Node, VolPar)) { SUMA_RETURN(NOPE); }
01720 
01721    SUMA_RETURN(YUP);
01722 }
01723 
01724 SUMA_Boolean SUMA_vec_3dmm_to_dicomm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01725 {
01726    static char FuncName[]={"SUMA_vec_3dmm_to_dicomm"};
01727    THD_fvec3 fv, iv;
01728    int i, id;
01729    SUMA_SurfaceObject SO;
01730    
01731    SUMA_ENTRY;
01732 
01733    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01734    /* create dummy struct */
01735    SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
01736 
01737    for (i=0; i < SO.N_Node; ++i) {
01738       id = i * SO.NodeDim;
01739       iv.xyz[0] = SO.NodeList[id] ;
01740       iv.xyz[1] = SO.NodeList[id+1] ;
01741       iv.xyz[2] = SO.NodeList[id+2] ;
01742 
01743       /* change mm to RAI coords */
01744       fv = SUMA_THD_3dmm_to_dicomm( SO.VolPar->xxorient, SO.VolPar->yyorient, SO.VolPar->zzorient,  iv );
01745       SO.NodeList[id] = fv.xyz[0];
01746       SO.NodeList[id+1] = fv.xyz[1];
01747       SO.NodeList[id+2] = fv.xyz[2];
01748    }
01749 
01750    SUMA_RETURN(YUP);
01751 }   
01752 
01753 SUMA_Boolean SUMA_vec_dicomm_to_3dmm (float *NodeList, int N_Node, SUMA_VOLPAR *VolPar)
01754 {
01755    static char FuncName[]={"SUMA_vec_dicomm_to_3dmm"};
01756    THD_fvec3 fv, iv;
01757    int i, id;
01758    SUMA_SurfaceObject SO;
01759 
01760    SUMA_ENTRY;
01761 
01762    if (!NodeList || !VolPar) { SUMA_SL_Err("Null NodeList || Null VolPar"); SUMA_RETURN(NOPE); }
01763    /* create dummy struct */
01764    SO.NodeList = NodeList; SO.N_Node = N_Node; SO.VolPar = VolPar; SO.NodeDim = 3;
01765 
01766    for (i=0; i < SO.N_Node; ++i) {
01767       id = i * SO.NodeDim;
01768 
01769       iv.xyz[0] = SO.NodeList[id] ;
01770       iv.xyz[1] = SO.NodeList[id+1] ;
01771       iv.xyz[2] = SO.NodeList[id+2] ;
01772 
01773       /* change mm to RAI coords */
01774       fv = SUMA_THD_dicomm_to_3dmm( &SO, iv );
01775       SO.NodeList[id] = fv.xyz[0];
01776       SO.NodeList[id+1] = fv.xyz[1];
01777       SO.NodeList[id+2] = fv.xyz[2];
01778    }
01779 
01780    SUMA_RETURN(YUP);
01781 }
01782 
01783 
01784 /*!
01785    \brief transforms XYZ coordinates from  AFNI'S RAI 
01786    talairach space to MNI space in in RAI or (LPI) 
01787    http://www.mrc-cbu.cam.ac.uk/Imaging/mnispace.html.
01788    see Bob's THD_tta_to_mni in thd_mnicoords.c
01789    
01790    \param NodeList (float *) vector of coordinates, 3 * N_Node long
01791    \param N_Node (int) number of nodes in vector above
01792    \param Coord (char *) one of two options
01793                   "RAI" meaning the MNI coordinates are to be in RAI
01794                   "LPI" meaning the MNI coordinates are to be in LPI
01795    \return flag (SUMA_Boolean) NOPE = failed
01796    
01797    - This function should be rewritten to call THD_tta_to_mni (although you'll have to deal 
01798    with the flipping option
01799    
01800    
01801 */
01802 SUMA_Boolean SUMA_AFNItlrc_toMNI(float *NodeList, int N_Node, char *Coord)
01803 {
01804    static char FuncName[]={"SUMA_AFNItlrc_toMNI"};
01805    SUMA_Boolean DoFlip = NOPE;
01806    int i, i3;
01807    float mx = 0.0,my = 0.0,mz  = 0.0, tx = 0.0,ty = 0.0,tz = 0.0 ;
01808    SUMA_Boolean LocalHead = NOPE;
01809    
01810    if (strcmp(Coord,"RAI") == 0) DoFlip = NOPE;
01811    else if (strcmp(Coord,"LPI") == 0) DoFlip = YUP;
01812    else {
01813       SUMA_S_Err("Can't do. Either RAI or LPI");
01814       SUMA_RETURN(NOPE);
01815    }   
01816 
01817    
01818    for (i=0; i< N_Node; ++i) {
01819       i3 = 3*i;
01820       if (DoFlip) {
01821          if (!i) SUMA_LH("Flipping to LPI");
01822          tx = - NodeList[i3]; ty = -NodeList[i3+1] ;  /* flip xy from RAI to LPI */
01823          tz = NodeList[i3+2] ;
01824       } else {
01825          if (!i) SUMA_LH("No Flipping, RAI maintained");
01826          tx =  NodeList[i3]; ty = NodeList[i3+1] ;  /* flip xy from RAI to LPI */
01827          tz =  NodeList[i3+2] ;
01828       }
01829       mx = 1.01010 * tx ;
01830       my = 1.02962 * ty - 0.05154 * tz ;
01831       mz = 0.05434 * ty + 1.08554 * tz ;
01832       if( mz < 0.0 ) mz *= 1.09523 ;
01833       NodeList[i3] = mx;
01834       NodeList[i3+1] = my;
01835       NodeList[i3+2] = mz;
01836    }
01837    
01838    
01839    SUMA_RETURN(YUP);
01840 }
01841 /*!
01842 
01843    \brief transforms XYZ coordinates by transfrom in warp.
01844    ans = SUMA_AFNI_forward_warp_xyz( warp , XYZv,  N);
01845    
01846    \param warp (THD_warp *) afni warp structure
01847    \param XYZv (float *) pointer to coordinates vector.
01848    \param N (int) number of XYZ triplets in XYZv.
01849       The ith X,Y, Z coordinates would be XYZv[3i], XYZv[3i+1],XYZv[3i+2]   
01850    \return ans (YUP/NOPE) NOPE: NULL warp or bad warp->type
01851    
01852    - The following functions are adapted from afni.c & Vecwarp.c
01853 */
01854 SUMA_Boolean SUMA_AFNI_forward_warp_xyz( THD_warp * warp , float *xyzv, int N)
01855 {
01856    static char FuncName[]={"SUMA_AFNI_forward_warp_xyz"};
01857    THD_fvec3 new_fv, old_fv;
01858    int i, i3;
01859    
01860    SUMA_ENTRY;
01861    
01862    if( warp == NULL ) {
01863       fprintf (SUMA_STDERR, "Error %s: NULL warp.\n", FuncName);
01864       SUMA_RETURN (NOPE) ;
01865    }
01866 
01867    switch( warp->type ){
01868 
01869       default: 
01870          fprintf (SUMA_STDERR, "Error %s: bad warp->type\n", FuncName);
01871          SUMA_RETURN (NOPE); 
01872          break ;
01873 
01874       case WARP_TALAIRACH_12_TYPE:{
01875          THD_linear_mapping map ;
01876          int iw ;
01877 
01878          /* forward transform each possible case,
01879             and test if result is in bot..top of defined map */
01880 
01881          for (i=0; i < N; ++i) {
01882             i3 = 3*i;
01883             old_fv.xyz[0] = xyzv[i3];
01884             old_fv.xyz[1] = xyzv[i3+1];
01885             old_fv.xyz[2] = xyzv[i3+2];
01886             
01887             for( iw=0 ; iw < 12 ; iw++ ){
01888                map    = warp->tal_12.warp[iw] ;
01889                new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
01890 
01891                if( new_fv.xyz[0] >= map.bot.xyz[0] &&
01892                    new_fv.xyz[1] >= map.bot.xyz[1] &&
01893                    new_fv.xyz[2] >= map.bot.xyz[2] &&
01894                    new_fv.xyz[0] <= map.top.xyz[0] &&
01895                    new_fv.xyz[1] <= map.top.xyz[1] &&
01896                    new_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
01897             }
01898             xyzv[i3] = new_fv.xyz[0];
01899             xyzv[i3+1] = new_fv.xyz[1];
01900             xyzv[i3+2] = new_fv.xyz[2];
01901             
01902          }
01903       }
01904       break ;
01905 
01906       case WARP_AFFINE_TYPE:{
01907          THD_linear_mapping map = warp->rig_bod.warp ;
01908          for (i=0; i < N; ++i) {
01909             i3 = 3*i;
01910             old_fv.xyz[0] = xyzv[i3];
01911             old_fv.xyz[1] = xyzv[i3+1];
01912             old_fv.xyz[2] = xyzv[i3+2];
01913             new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
01914             xyzv[i3] = new_fv.xyz[0];
01915             xyzv[i3+1] = new_fv.xyz[1];
01916             xyzv[i3+2] = new_fv.xyz[2];
01917          }
01918       }
01919       break ;
01920 
01921    }
01922    SUMA_RETURN(YUP);
01923 }