SUMA_VolData.c Source File

00001 /*! Dealings with volume data */
00002 #include "SUMA_suma.h"
00003 
00004 extern SUMA_CommonFields *SUMAg_CF; 
00005 
00006 
00007 THD_fvec3 SUMA_THD_3dfind_to_3dmm( SUMA_SurfaceObject *SO, THD_fvec3 iv );
00008 THD_fvec3 SUMA_THD_3dind_to_3dmm( SUMA_SurfaceObject *SO, THD_ivec3 iv );
00009 THD_fvec3 SUMA_THD_3dmm_to_3dfind( SUMA_SurfaceObject *SO , THD_fvec3 fv );
00010 THD_ivec3 SUMA_THD_3dmm_to_3dind( SUMA_SurfaceObject *SO  , THD_fvec3 fv );
00011 THD_fvec3 SUMA_THD_3dmm_to_dicomm( SUMA_SurfaceObject *SO , THD_fvec3 imv );
00012 THD_fvec3 SUMA_THD_dicomm_to_3dmm( SUMA_SurfaceObject *SO , THD_fvec3 dicv );
00013 
00014 /*!
00015         Returns the attributes of the surface's volume parent 
00016         A volume parent is defined as:
00017                 The volume data set used to generate the surface. aka the master SPGR
00018         \ret NULL for troubles
00019 */
00020 SUMA_VOLPAR *SUMA_Alloc_VolPar (void)
00021 {
00022         static char FuncName[]={"SUMA_Alloc_VolPar"};
00023         SUMA_VOLPAR *VP;
00024 
00025         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00026 
00027         VP = (SUMA_VOLPAR *)SUMA_malloc(sizeof(SUMA_VOLPAR));
00028         if (VP == NULL) {
00029                 fprintf(SUMA_STDERR,"Error SUMA_Alloc_VolPar: Failed to allocate for VolPar\n");
00030                 SUMA_RETURN (NULL);
00031         }
00032         SUMA_RETURN(VP);
00033 }
00034 SUMA_Boolean SUMA_Free_VolPar (SUMA_VOLPAR *VP)
00035 {
00036         static char FuncName[]={"SUMA_Free_VolPar"};
00037         
00038         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00039 
00040         if (VP->prefix != NULL) SUMA_free(VP->prefix);
00041         if (VP->filecode != NULL) SUMA_free(VP->filecode);
00042         if (VP->dirname != NULL) SUMA_free(VP->dirname);
00043         if (VP->idcode_str != NULL) SUMA_free(VP->idcode_str);
00044         if (VP->idcode_date != NULL) SUMA_free(VP->idcode_date);
00045         if (VP->VOLREG_CENTER_OLD != NULL) SUMA_free(VP->VOLREG_CENTER_OLD);
00046         if (VP->VOLREG_CENTER_BASE != NULL) SUMA_free(VP->VOLREG_CENTER_BASE);
00047         if (VP->VOLREG_MATVEC != NULL) SUMA_free(VP->VOLREG_MATVEC);
00048         if (VP != NULL) SUMA_free(VP);
00049         SUMA_RETURN (YUP);
00050 }
00051 
00052 SUMA_VOLPAR *SUMA_VolPar_Attr (char *volparent_name)
00053 {
00054         ATR_float *atr;
00055         static char FuncName[]={"SUMA_VolPar_Attr"};
00056         SUMA_VOLPAR *VP;
00057         THD_3dim_dataset *dset;
00058         int ii;
00059         MCW_idcode idcode;
00060                 
00061         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00062 
00063         /* read the header of the parent volume */
00064         dset = THD_open_dataset(volparent_name);
00065         if (dset == NULL) {
00066                 fprintf (SUMA_STDERR,"Error %s: Could not read %s\n", FuncName, volparent_name);
00067                 SUMA_RETURN (NULL);
00068         }
00069         
00070         /* allocate for VP */
00071         VP = SUMA_Alloc_VolPar();
00072         if (VP == NULL) {
00073                 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Alloc_VolPar\n", FuncName);
00074                 SUMA_RETURN (NULL);
00075         }
00076         
00077         /* retain pertinent info */
00078         VP->isanat = ISANAT(dset);
00079         VP->nx = DSET_NX(dset);
00080         VP->ny = DSET_NY(dset);
00081         VP->nz = DSET_NZ(dset);
00082         VP->dx = DSET_DX(dset);
00083         VP->dy = DSET_DY(dset);
00084         VP->dz = DSET_DZ(dset);
00085         VP->xorg = DSET_XORG(dset);
00086         VP->yorg = DSET_YORG(dset);
00087         VP->zorg = DSET_ZORG(dset);
00088         ii = strlen(DSET_PREFIX(dset));
00089         VP->prefix = (char *)SUMA_malloc(ii+1);
00090         ii = strlen(DSET_FILECODE(dset));
00091         VP->filecode = (char *)SUMA_malloc(ii+1);
00092         ii = strlen(DSET_DIRNAME(dset));
00093         VP->dirname = (char *)SUMA_malloc(ii+1);
00094         ii = strlen(dset->idcode.str);
00095         VP->idcode_str = (char *)SUMA_malloc(ii+1);
00096         ii = strlen(dset->idcode.date);
00097         VP->idcode_date = (char *)SUMA_malloc(ii+1);
00098         if (VP->prefix == NULL || VP->filecode == NULL || VP->idcode_date == NULL || VP->dirname == NULL || VP->idcode_str == NULL) {
00099                 fprintf(SUMA_STDERR,"Error %s: Failed to allocate for strings. Kill me, please.\n", FuncName);
00100                 SUMA_Free_VolPar(VP);
00101                 SUMA_RETURN (NULL);
00102         }
00103         VP->prefix = strcpy(VP->prefix, DSET_PREFIX(dset));
00104         VP->filecode = strcpy(VP->filecode, DSET_FILECODE(dset));
00105         VP->dirname = strcpy(VP->dirname, DSET_DIRNAME(dset));
00106         VP->idcode_str = strcpy(VP->idcode_str, dset->idcode.str);
00107         VP->idcode_date = strcpy(VP->idcode_date, dset->idcode.date);
00108         VP->xxorient = dset->daxes->xxorient;
00109         VP->yyorient = dset->daxes->yyorient;
00110         VP->zzorient = dset->daxes->zzorient;
00111 
00112         /* Get the volreg matrix if possible*/
00113         atr = THD_find_float_atr( dset->dblk , "VOLREG_MATVEC_000000" ) ;
00114         if (atr == NULL) {
00115                 VP->VOLREG_MATVEC = NULL;
00116         }else {
00117                 VP->VOLREG_MATVEC = (float *)SUMA_calloc(12, sizeof(float));
00118                 if (VP->VOLREG_MATVEC != NULL) {
00119                         if (atr->nfl == 12) {
00120                                 for (ii=0; ii<12; ++ii) VP->VOLREG_MATVEC[ii] = atr->fl[ii];
00121                         } else {        
00122                                 fprintf(SUMA_STDERR,"Error %s: VOLREG_MATVEC_000000 does not have 12 elements.\n", FuncName);
00123                         }
00124                 } else {
00125                         fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_MATVEC\n", FuncName);
00126                 }
00127         }
00128         /* Get the center base coordinates */
00129         atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_BASE");
00130         if (atr == NULL) {
00131                 VP->VOLREG_CENTER_BASE = NULL;
00132         } else {
00133                 VP->VOLREG_CENTER_BASE = (float *)SUMA_calloc(3, sizeof(float));
00134                 if (VP->VOLREG_CENTER_BASE != NULL) {
00135                         if (atr->nfl == 3) {
00136                                 for (ii=0; ii<3; ++ii) VP->VOLREG_CENTER_BASE[ii] = atr->fl[ii];
00137                         } else {        
00138                                 fprintf(SUMA_STDERR,"Error %s: VOLREG_CENTER_BASE does not have 12 elements.\n", FuncName);
00139                         }
00140                 } else {
00141                         fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_CENTER_BASE\n", FuncName);
00142                 }
00143         }
00144         
00145         /* VOLREG_CENTER_OLD  */
00146         atr = THD_find_float_atr( dset->dblk , "VOLREG_CENTER_OLD");
00147         if (atr == NULL) {
00148                 VP->VOLREG_CENTER_OLD = NULL;
00149         } else {
00150                 VP->VOLREG_CENTER_OLD = (float *)SUMA_calloc(3, sizeof(float));
00151                 if (VP->VOLREG_CENTER_OLD != NULL) {
00152                         if (atr->nfl == 3) {
00153                                 for (ii=0; ii<3; ++ii) VP->VOLREG_CENTER_OLD[ii] = atr->fl[ii];
00154                         } else {        
00155                                 fprintf(SUMA_STDERR,"Error %s: VOLREG_CENTER_OLD does not have 12 elements.\n", FuncName);
00156                         }
00157                 } else {
00158                         fprintf(SUMA_STDERR,"Error %s: Failed to allocate for VP->VOLREG_CENTER_OLD\n", FuncName);
00159                 }
00160         }
00161 
00162    /* now free the dset pointer */
00163    THD_delete_3dim_dataset( dset , False ) ;
00164    
00165         SUMA_RETURN (VP);
00166 }
00167 
00168 /*!
00169    \brief Form a string containing the info of the volume parent
00170    
00171    \param VP (SUMA_VOLPAR *) Volume parent structure.
00172    \return s (char *) pointer to NULL terminated string containing surface info.
00173    It is your responsability to free it.
00174    
00175    \sa SUMA_Show_VolPar
00176 */
00177 
00178 char *SUMA_VolPar_Info (SUMA_VOLPAR *VP)
00179 {
00180    static char FuncName[]={"SUMA_VolPar_Info"};
00181    char stmp[1000], *s=NULL;
00182    SUMA_STRING *SS = NULL;
00183    
00184    if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00185 
00186    SS = SUMA_StringAppend (NULL, NULL);
00187    
00188    if (VP) { 
00189            sprintf (stmp,"\nVP contents:\n");
00190       SS = SUMA_StringAppend (SS, stmp);
00191            sprintf (stmp,"prefix: %s\tfilecode: %s\tdirname: %s\nId code str:%s\tID code date: %s\n", \
00192                    VP->prefix, VP->filecode, VP->dirname, VP->idcode_str, VP->idcode_date);
00193       SS = SUMA_StringAppend (SS, stmp);
00194            sprintf (stmp,"isanat: %d\n", VP->isanat);
00195       SS = SUMA_StringAppend (SS, stmp);
00196            sprintf (stmp,"Orientation: %d %d %d\n", \
00197                    VP->xxorient, VP->yyorient, VP->zzorient);
00198       SS = SUMA_StringAppend (SS, stmp);
00199            sprintf (stmp,"Origin: %f %f %f\n", \
00200                    VP->xorg, VP->yorg, VP->zorg);
00201            SS = SUMA_StringAppend (SS, stmp);
00202       sprintf (stmp,"Delta: %f %f %f\n", \
00203                    VP->dx, VP->dy, VP->dz);
00204       SS = SUMA_StringAppend (SS, stmp);
00205            sprintf (stmp,"N: %d %d %d\n",\
00206                    VP->nx, VP->ny, VP->nz);
00207       SS = SUMA_StringAppend (SS, stmp);
00208 
00209            if (VP->VOLREG_MATVEC != NULL) {
00210                    sprintf (stmp,"VP->VOLREG_MATVEC = \n\tMrot\tDelta\n");
00211          SS = SUMA_StringAppend (SS, stmp);
00212                    sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00213                    VP->VOLREG_MATVEC[0], VP->VOLREG_MATVEC[1], VP->VOLREG_MATVEC[2], VP->VOLREG_MATVEC[3]); 
00214          SS = SUMA_StringAppend (SS, stmp);
00215                    sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00216                    VP->VOLREG_MATVEC[4], VP->VOLREG_MATVEC[5], VP->VOLREG_MATVEC[6], VP->VOLREG_MATVEC[7]);
00217          SS = SUMA_StringAppend (SS, stmp);
00218                    sprintf (stmp,"|%f\t%f\t%f|\t|%f|\n", \
00219                    VP->VOLREG_MATVEC[8], VP->VOLREG_MATVEC[9], VP->VOLREG_MATVEC[10], VP->VOLREG_MATVEC[11]);
00220          SS = SUMA_StringAppend (SS, stmp);
00221            } else {
00222          sprintf (stmp,"VP->VOLREG_MATVEC = NULL\n");
00223          SS = SUMA_StringAppend (SS, stmp);
00224       }
00225 
00226            if (VP->VOLREG_CENTER_OLD != NULL) {
00227                    sprintf (stmp,"VP->VOLREG_CENTER_OLD = %f, %f, %f\n", \
00228                      VP->VOLREG_CENTER_OLD[0], VP->VOLREG_CENTER_OLD[1], VP->VOLREG_CENTER_OLD[2]); 
00229          SS = SUMA_StringAppend (SS, stmp);
00230            }else {
00231          sprintf (stmp,"VP->VOLREG_CENTER_OLD = NULL\n");
00232          SS = SUMA_StringAppend (SS, stmp);
00233       }
00234 
00235            if (VP->VOLREG_CENTER_BASE != NULL) {
00236                    sprintf (stmp,"VP->VOLREG_CENTER_BASE = %f, %f, %f\n", \
00237                    VP->VOLREG_CENTER_BASE[0], VP->VOLREG_CENTER_BASE[1], VP->VOLREG_CENTER_BASE[2]); 
00238          SS = SUMA_StringAppend (SS, stmp);
00239            } else {
00240          sprintf (stmp,"VP->VOLREG_CENTER_BASE = NULL\n");
00241          SS = SUMA_StringAppend (SS, stmp);
00242       }
00243    }else{
00244       sprintf (stmp, "NULL Volume Parent Pointer.\n");
00245       SS = SUMA_StringAppend (SS, stmp);
00246    }
00247    
00248    /* clean SS */
00249    SS = SUMA_StringAppend (SS, NULL);
00250    /* copy s pointer and free SS */
00251    s = SS->s;
00252    SUMA_free(SS); 
00253    
00254    SUMA_RETURN (s);
00255 }
00256 /*!
00257         Show the contents of SUMA_VOLPAR structure 
00258    \sa SUMA_VolPar_Info 
00259 */
00260 void SUMA_Show_VolPar(SUMA_VOLPAR *VP, FILE *Out)
00261 {
00262         static char FuncName[]={"SUMA_Show_VolPar"};
00263         char *s;
00264    
00265         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00266 
00267         if (Out == NULL) Out = SUMA_STDOUT;
00268         
00269    s =  SUMA_VolPar_Info(VP);
00270    
00271    if (s) {
00272       fprintf (Out, "%s", s);
00273       SUMA_free(s);
00274    }else {
00275       fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_VolPar_Info.\n", FuncName);
00276    }   
00277    
00278         SUMA_RETURNe;
00279                 
00280 }
00281 
00282 /*!
00283         \brief Transform the coordinates of a surface object to AFNI-DICOM convention and transform the coordinates by SO->VolPar
00284    ans = SUMA_Align_to_VolPar (SO, SF_Struct);
00285    \param SO (SUMA_SurfaceObject *)
00286    \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.
00287                               send NULL for all other surfaces.
00288    \return YUP/NOPE
00289    For SureFit and FreeSurfer surfaces, the coordinates are first set in RAI (DICOM) coordinate system before applying SO->VolPar.
00290    For other surface formats, SO->VolPar is applied to whatever coordinates are in SO->NodeList
00291 */
00292 SUMA_Boolean SUMA_Align_to_VolPar (SUMA_SurfaceObject *SO, void * S_Struct)
00293 {
00294         static char FuncName[]={"SUMA_Align_to_VolPar"};
00295         char  orcode[3];
00296         float xx, yy, zz;
00297         THD_coorder * cord_surf, *cord_RAI;
00298         int i, ND, id;
00299         SUMA_SureFit_struct *SF;
00300         SUMA_FreeSurfer_struct *FS;
00301         
00302         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00303 
00304         /* allocate for cord structure */
00305         cord_surf = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
00306         cord_RAI = (THD_coorder *)SUMA_malloc(sizeof(THD_coorder));
00307         if (cord_surf == NULL || cord_RAI == NULL) {
00308                 fprintf(SUMA_STDERR,"Error %s: failed to allocate for cord\n", FuncName);
00309                 SUMA_RETURN (NOPE);
00310         }
00311         
00312         /* do the dance to get cord for RAI */
00313         THD_coorder_fill( NULL, cord_RAI);
00314         ND = SO->NodeDim;
00315         switch (SO->FileType) {
00316                 case SUMA_INVENTOR_GENERIC:
00317                         /* Do nothing */
00318                         break;
00319                 case SUMA_FREE_SURFER:
00320                         /* For free surfer, all you need to do is 
00321                          go from LPI to RAI (DICOM)*/
00322                         sprintf(orcode,"LPI");
00323                         THD_coorder_fill(orcode , cord_surf); 
00324                         /*loop over XYZs and change them to dicom*/
00325                         for (i=0; i < SO->N_Node; ++i) {
00326                                 id = i * ND;
00327                                 THD_coorder_to_dicom (cord_surf, &(SO->NodeList[id]), &(SO->NodeList[id+1]), &(SO->NodeList[id+2])); 
00328                         }
00329                         break;
00330                 case SUMA_SUREFIT:
00331                         /* For SureFit, coordinates are actually a float version of the indices */
00332                         SF = (SUMA_SureFit_struct *)S_Struct;
00333                         {       THD_fvec3 fv, iv;
00334                                 float D[3];
00335                                 /* Calcluate Delta caused by cropping */
00336                                 for (i=0; i < 3; ++i) D[i] = SF->AC_WholeVolume[i] - SF->AC[i];
00337                                 /* fprintf (SUMA_STDERR,"%s: Shift Values: [%f, %f, %f]\n", FuncName, D[0], D[1], D[2]); */
00338                                 for (i=0; i < SO->N_Node; ++i) {
00339                                         id = i * ND;
00340                                         /* change float indices to mm coords */
00341                                         iv.xyz[0] = SO->NodeList[id] + D[0];
00342                                         iv.xyz[1] = SO->NodeList[id+1] + D[1];
00343                                         iv.xyz[2] = SO->NodeList[id+2] + D[2];
00344                                         fv = SUMA_THD_3dfind_to_3dmm( SO, iv );
00345                                         
00346                                         /* change mm to RAI coords */
00347                                         iv = SUMA_THD_3dmm_to_dicomm( SO , fv );
00348                                         SO->NodeList[id] = iv.xyz[0];
00349                                         SO->NodeList[id+1] = iv.xyz[1];
00350                                         SO->NodeList[id+2] = iv.xyz[2];
00351                                 }
00352                         }
00353                         break;
00354                 default:
00355                         fprintf(SUMA_STDERR,"Warning %s: Unknown SO->FileType. Assuming coordinates are in DICOM already.\n", FuncName);
00356                         break;
00357         }
00358         
00359    if (!SUMA_Apply_VolReg_Trans (SO)) {
00360       fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Apply_VolReg_Trans.\n", FuncName);
00361       SUMA_RETURN (NOPE);
00362    }
00363 
00364         SUMA_RETURN (YUP);
00365 }
00366 
00367 /*!
00368 \brief applies the transform in SO->VolPar to SO->NodeList. This only makes sense if 
00369 SO->NodeList is in DICOM to begin with...
00370 \param SO (SUMA_SurfaceObject *) You better have NodeList and VolPar in there...
00371 \return YUP/NOPE
00372 \sa SUMA_Align_to_VolPar
00373 NOTE: The field SO->VOLREG_APPLIED is set to NOPE if this function fails, YUP otherwise
00374 */
00375 SUMA_Boolean SUMA_Apply_VolReg_Trans (SUMA_SurfaceObject *SO)
00376 {
00377         static char FuncName[]={"SUMA_Apply_VolReg_Trans"};
00378    int i, ND, id;
00379    SUMA_Boolean Bad=YUP;
00380    
00381    if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00382 
00383    if (SO->VOLREG_APPLIED) {
00384       fprintf (SUMA_STDERR,"Error %s: Volreg already applied. Nothing done.\n", FuncName);
00385       SUMA_RETURN (NOPE);
00386    }
00387 
00388         ND = SO->NodeDim;
00389    
00390    /* perform the rotation needed to align the surface with the current experiment's data */
00391         if (SO->VolPar->VOLREG_MATVEC != NULL || SO->VolPar->VOLREG_CENTER_OLD != NULL || SO->VolPar->VOLREG_CENTER_BASE != NULL) {
00392                 Bad = NOPE;
00393                 if (SO->VolPar->VOLREG_MATVEC == NULL) {
00394                         fprintf(SUMA_STDERR,"Error %s: VP->VOLREG_MATVEC = NULL. Cannot perform alignment.\n", FuncName);
00395                         Bad = YUP;
00396                 }
00397                 if (SO->VolPar->VOLREG_CENTER_OLD == NULL) {
00398                         fprintf(SUMA_STDERR,"Error %s: VP->VOLREG_CENTER_OLD = NULL. Cannot perform alignment.\n", FuncName);
00399                         Bad = YUP;
00400                 }
00401                 if (SO->VolPar->VOLREG_CENTER_BASE == NULL) {
00402                         fprintf(SUMA_STDERR,"Error %s: VP->VOLREG_CENTER_BASE = NULL. Cannot perform alignment.\n", FuncName);
00403                         Bad = YUP;
00404                 }
00405         }
00406 
00407         if (!Bad && SO->VolPar->VOLREG_MATVEC != NULL) {
00408                 float Mrot[3][3], Delta[3], x, y, z, NetShift[3];
00409                 
00410                 /* fillerup*/
00411                 Mrot[0][0] = SO->VolPar->VOLREG_MATVEC[0];
00412                 Mrot[0][1] = SO->VolPar->VOLREG_MATVEC[1];
00413                 Mrot[0][2] = SO->VolPar->VOLREG_MATVEC[2];
00414                 Delta[0]   = SO->VolPar->VOLREG_MATVEC[3];
00415                 Mrot[1][0] = SO->VolPar->VOLREG_MATVEC[4];
00416                 Mrot[1][1] = SO->VolPar->VOLREG_MATVEC[5];
00417                 Mrot[1][2] = SO->VolPar->VOLREG_MATVEC[6];
00418                 Delta[1]   = SO->VolPar->VOLREG_MATVEC[7];      
00419                 Mrot[2][0] = SO->VolPar->VOLREG_MATVEC[8];
00420                 Mrot[2][1] = SO->VolPar->VOLREG_MATVEC[9];
00421                 Mrot[2][2] = SO->VolPar->VOLREG_MATVEC[10];
00422                 Delta[2]   = SO->VolPar->VOLREG_MATVEC[11];
00423                 
00424                 NetShift[0] = SO->VolPar->VOLREG_CENTER_BASE[0] + Delta[0];
00425                 NetShift[1] = SO->VolPar->VOLREG_CENTER_BASE[1] + Delta[1];
00426                 NetShift[2] = SO->VolPar->VOLREG_CENTER_BASE[2] + Delta[2];
00427                 
00428                 /*
00429                 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,\
00430                                         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], \
00431                                         Delta[0], Delta[1], Delta[2]);
00432                 fprintf (SUMA_STDERR,"VOLREG_CENTER_BASE = [%f %f %f]. VOLREG_CENTER_OLD = [%f %f %f]\n", \
00433                         SO->VolPar->VOLREG_CENTER_BASE[0], SO->VolPar->VOLREG_CENTER_BASE[1], SO->VolPar->VOLREG_CENTER_BASE[2], \
00434                         SO->VolPar->VOLREG_CENTER_OLD[0], SO->VolPar->VOLREG_CENTER_OLD[1], SO->VolPar->VOLREG_CENTER_OLD[2]);
00435                 */
00436                 
00437                 for (i=0; i < SO->N_Node; ++i) {
00438                         id = ND * i;
00439                         /* zero the center */ 
00440                         x = SO->NodeList[id] - SO->VolPar->VOLREG_CENTER_OLD[0];
00441                         y = SO->NodeList[id+1] - SO->VolPar->VOLREG_CENTER_OLD[1];
00442                         z = SO->NodeList[id+2] - SO->VolPar->VOLREG_CENTER_OLD[2];
00443                         
00444                         /* Apply the rotation matrix XYZn = Mrot x XYZ*/
00445                         SO->NodeList[id] = Mrot[0][0] * x + Mrot[0][1] * y + Mrot[0][2] * z;
00446                         SO->NodeList[id+1] = Mrot[1][0] * x + Mrot[1][1] * y + Mrot[1][2] * z;
00447                         SO->NodeList[id+2] = Mrot[2][0] * x + Mrot[2][1] * y + Mrot[2][2] * z;
00448                         
00449                         /*apply netshift*/
00450                         SO->NodeList[id] += NetShift[0];
00451                         SO->NodeList[id+1] += NetShift[1];
00452                         SO->NodeList[id+2] += NetShift[2];
00453                 }
00454                 SO->VOLREG_APPLIED = YUP;       
00455         } else
00456                 SO->VOLREG_APPLIED = NOPE;
00457       
00458         SUMA_RETURN (YUP);
00459 }
00460 
00461 /*! The following functions are adaptations from thd_coords.c 
00462 They are modified to avoid the use of dset data structure 
00463 which is not fully preserved in the Surface Object Data structure.
00464 
00465 Stolen Comment:
00466 ====================================================================
00467    3D coordinate conversion routines;
00468      tags for coordinate systems:
00469        fdind  = FD_brick voxel indices (ints)
00470        fdfind = FD_brick voxel indices (floats - added 30 Aug 2001)
00471        3dind  = THD_3dim_dataset voxel indices (ints)
00472        3dfind = THD_3dim_dataset floating point voxel indices
00473                   (used for subvoxel resolution)
00474        3dmm   = THD_3dim_dataset millimetric coordinates (floats)
00475        dicomm = DICOM 3.0 millimetric coordinates (floats)
00476 
00477      The 3dmm coordinate measurements are oriented the same way
00478      as the dicomm coordinates (which means that the ??del values
00479      can be negative if the voxel axes are reflected), but the 3dmm
00480      coordinates are in general a permutation of the DICOM coordinates.
00481 
00482      These routines all take as input and produce as output
00483      THD_?vec3 (i=int,f=float) structs.
00484 ======================================================================
00485 
00486 */
00487 
00488 THD_fvec3 SUMA_THD_3dfind_to_3dmm( SUMA_SurfaceObject *SO, 
00489                               THD_fvec3 iv )
00490 {
00491         static char FuncName[]={"SUMA_THD_3dfind_to_3dmm"};
00492    THD_fvec3     fv ;
00493 
00494         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00495 
00496    fv.xyz[0] = SO->VolPar->xorg + iv.xyz[0] * SO->VolPar->dx ;
00497    fv.xyz[1] = SO->VolPar->yorg + iv.xyz[1] * SO->VolPar->dy ;
00498    fv.xyz[2] = SO->VolPar->zorg + iv.xyz[2] * SO->VolPar->dz ;
00499    SUMA_RETURN(fv) ;
00500 }
00501 
00502 /*!------------------------------------------------------------------*/
00503 
00504 THD_fvec3 SUMA_THD_3dind_to_3dmm( SUMA_SurfaceObject *SO,
00505                                          THD_ivec3 iv )
00506 {
00507         static char FuncName[]={"SUMA_THD_3dind_to_3dmm"};
00508         THD_fvec3     fv ;
00509 
00510         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00511 
00512    fv.xyz[0] = SO->VolPar->xorg + iv.ijk[0] * SO->VolPar->dx ;
00513    fv.xyz[1] = SO->VolPar->yorg + iv.ijk[1] * SO->VolPar->dy ;
00514    fv.xyz[2] = SO->VolPar->zorg + iv.ijk[2] * SO->VolPar->dz ;
00515    SUMA_RETURN(fv) ;
00516 }
00517 
00518 /*!--------------------------------------------------------------------*/
00519 
00520 THD_fvec3 SUMA_THD_3dmm_to_3dfind( SUMA_SurfaceObject *SO ,
00521                               THD_fvec3 fv )
00522 {
00523         static char FuncName[]={"SUMA_THD_3dmm_to_3dfind"};
00524    THD_fvec3     iv ;
00525 
00526         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00527 
00528 
00529    iv.xyz[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx ;
00530    iv.xyz[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy ;
00531    iv.xyz[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz ;
00532 
00533         if( iv.xyz[0] < 0            ) iv.xyz[0] = 0 ;
00534    else if( iv.xyz[0] > SO->VolPar->nx-1 ) iv.xyz[0] = SO->VolPar->nx-1 ;
00535 
00536         if( iv.xyz[1] <  0           ) iv.xyz[1] = 0 ;
00537    else if( iv.xyz[1] > SO->VolPar->ny-1 ) iv.xyz[1] = SO->VolPar->ny-1 ;
00538 
00539         if( iv.xyz[2] < 0            ) iv.xyz[2] = 0 ;
00540    else if( iv.xyz[2] > SO->VolPar->nz-1 ) iv.xyz[2] = SO->VolPar->nz-1 ;
00541 
00542    SUMA_RETURN(iv) ;
00543 }
00544 
00545 /*!--------------------------------------------------------------------*/
00546 
00547 THD_ivec3 SUMA_THD_3dmm_to_3dind( SUMA_SurfaceObject *SO  ,
00548                              THD_fvec3 fv )
00549 {
00550         static char FuncName[]={"SUMA_THD_3dmm_to_3dind"};
00551    THD_ivec3     iv ;
00552 
00553         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00554 
00555    iv.ijk[0] = (fv.xyz[0] - SO->VolPar->xorg) / SO->VolPar->dx + 0.499 ;
00556    iv.ijk[1] = (fv.xyz[1] - SO->VolPar->yorg) / SO->VolPar->dy + 0.499 ;
00557    iv.ijk[2] = (fv.xyz[2] - SO->VolPar->zorg) / SO->VolPar->dz + 0.499 ;
00558 
00559         if( iv.ijk[0] < 0            ) iv.ijk[0] = 0 ;
00560    else if( iv.ijk[0] > SO->VolPar->nx-1 ) iv.ijk[0] = SO->VolPar->nx-1 ;
00561 
00562         if( iv.ijk[1] < 0            ) iv.ijk[1] = 0 ;
00563    else if( iv.ijk[1] > SO->VolPar->ny-1 ) iv.ijk[1] = SO->VolPar->ny-1 ;
00564 
00565         if( iv.ijk[2] < 0            ) iv.ijk[2] = 0 ;
00566    else if( iv.ijk[2] > SO->VolPar->nz-1 ) iv.ijk[2] = SO->VolPar->nz-1 ;
00567 
00568    SUMA_RETURN(iv) ;
00569 }
00570 
00571 /*!---------------------------------------------------------------------
00572    convert from input image oriented x,y,z to Dicom x,y,z
00573      (x axis = R->L , y axis = A->P , z axis = I->S)
00574 
00575    N.B.: image distances are oriented the same as Dicom,
00576          just in a permuted order.
00577 -----------------------------------------------------------------------*/
00578 
00579 THD_fvec3 SUMA_THD_3dmm_to_dicomm( SUMA_SurfaceObject *SO ,
00580                               THD_fvec3 imv )
00581 {
00582         static char FuncName[]={"SUMA_THD_3dmm_to_dicomm"};   
00583         THD_fvec3 dicv ;
00584    float xim,yim,zim , xdic = 0.0, ydic = 0.0, zdic = 0.0 ;
00585 
00586         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00587 
00588    xim = imv.xyz[0] ; yim = imv.xyz[1] ; zim = imv.xyz[2] ;
00589 
00590    switch( SO->VolPar->xxorient ){
00591       case ORI_R2L_TYPE:
00592       case ORI_L2R_TYPE: xdic = xim ; break ;
00593       case ORI_P2A_TYPE:
00594       case ORI_A2P_TYPE: ydic = xim ; break ;
00595       case ORI_I2S_TYPE:
00596       case ORI_S2I_TYPE: zdic = xim ; break ;
00597 
00598       default: 
00599                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00600                         exit (1);
00601    }
00602 
00603    switch( SO->VolPar->yyorient ){
00604       case ORI_R2L_TYPE:
00605       case ORI_L2R_TYPE: xdic = yim ; break ;
00606       case ORI_P2A_TYPE:
00607       case ORI_A2P_TYPE: ydic = yim ; break ;
00608       case ORI_I2S_TYPE:
00609       case ORI_S2I_TYPE: zdic = yim ; break ;
00610 
00611       default: 
00612                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00613                         exit (1);
00614    }
00615 
00616    switch( SO->VolPar->zzorient ){
00617       case ORI_R2L_TYPE:
00618       case ORI_L2R_TYPE: xdic = zim ; break ;
00619       case ORI_P2A_TYPE:
00620       case ORI_A2P_TYPE: ydic = zim ; break ;
00621       case ORI_I2S_TYPE:
00622       case ORI_S2I_TYPE: zdic = zim ; break ;
00623 
00624       default: 
00625                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00626                         exit (1);
00627    }
00628 
00629    dicv.xyz[0] = xdic ; dicv.xyz[1] = ydic ; dicv.xyz[2] = zdic ;
00630    SUMA_RETURN(dicv) ;
00631 }
00632 
00633 /*!---------------------------------------------------------------------
00634    convert to input image oriented x,y,z from Dicom x,y,z
00635 -----------------------------------------------------------------------*/
00636 
00637 THD_fvec3 SUMA_THD_dicomm_to_3dmm( SUMA_SurfaceObject *SO ,
00638                               THD_fvec3 dicv )
00639 {
00640    static char FuncName[]={"SUMA_THD_dicomm_to_3dmm"};
00641         THD_fvec3 imv ;
00642    float xim,yim,zim , xdic,ydic,zdic ;
00643 
00644         if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00645 
00646    xdic = dicv.xyz[0] ; ydic = dicv.xyz[1] ; zdic = dicv.xyz[2] ;
00647 
00648    switch( SO->VolPar->xxorient ){
00649       case ORI_R2L_TYPE:
00650       case ORI_L2R_TYPE: xim = xdic ; break ;
00651       case ORI_P2A_TYPE:
00652       case ORI_A2P_TYPE: xim = ydic ; break ;
00653       case ORI_I2S_TYPE:
00654       case ORI_S2I_TYPE: xim = zdic ; break ;
00655 
00656       default: 
00657                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00658                         exit (1);
00659    }
00660 
00661    switch( SO->VolPar->yyorient ){
00662       case ORI_R2L_TYPE:
00663       case ORI_L2R_TYPE: yim = xdic ; break ;
00664       case ORI_P2A_TYPE:
00665       case ORI_A2P_TYPE: yim = ydic ; break ;
00666       case ORI_I2S_TYPE:
00667       case ORI_S2I_TYPE: yim = zdic ; break ;
00668 
00669       default: 
00670                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00671                         exit (1);
00672    }
00673 
00674    switch( SO->VolPar->zzorient ){
00675       case ORI_R2L_TYPE:
00676       case ORI_L2R_TYPE: zim = xdic ; break ;
00677       case ORI_P2A_TYPE:
00678       case ORI_A2P_TYPE: zim = ydic ; break ;
00679       case ORI_I2S_TYPE:
00680       case ORI_S2I_TYPE: zim = zdic ; break ;
00681 
00682       default: 
00683                         fprintf(SUMA_STDERR, "SUMA_THD_3dmm_to_dicomm: illegal xxorient code.\n Exiting.") ;
00684                         exit (1);
00685    }
00686 
00687    imv.xyz[0] = xim ; imv.xyz[1] = yim ; imv.xyz[2] = zim ;
00688    SUMA_RETURN(imv) ;
00689 }
00690 
00691 /*!
00692 
00693    \brief transforms XYZ coordinates by transfrom in warp.
00694    ans = SUMA_AFNI_forward_warp_xyz( warp , XYZv,  N);
00695    
00696    \param warp (THD_warp *) afni warp structure
00697    \param XYZv (float *) pointer to coordinates vector.
00698    \param N (int) number of XYZ triplets in XYZv.
00699       The ith X,Y, Z coordinates would be XYZv[3i], XYZv[3i+1],XYZv[3i+2]   
00700    \return ans (YUP/NOPE) NOPE: NULL warp or bad warp->type
00701    
00702    - The following functions are adapted from afni.c & Vecwarp.c
00703 */
00704 SUMA_Boolean SUMA_AFNI_forward_warp_xyz( THD_warp * warp , float *xyzv, int N)
00705 {
00706    static char FuncName[]={"SUMA_AFNI_forward_warp_xyz"};
00707    THD_fvec3 new_fv, old_fv;
00708    int i, i3;
00709    
00710    if (SUMAg_CF->InOut_Notify) SUMA_DBG_IN_NOTIFY(FuncName);
00711    
00712    if( warp == NULL ) {
00713       fprintf (SUMA_STDERR, "Error %s: NULL warp.\n", FuncName);
00714       SUMA_RETURN (NOPE) ;
00715    }
00716 
00717    switch( warp->type ){
00718 
00719       default: 
00720          fprintf (SUMA_STDERR, "Error %s: bad warp->type\n", FuncName);
00721          SUMA_RETURN (NOPE); 
00722          break ;
00723 
00724       case WARP_TALAIRACH_12_TYPE:{
00725          THD_linear_mapping map ;
00726          int iw ;
00727 
00728          /* forward transform each possible case,
00729             and test if result is in bot..top of defined map */
00730 
00731          for (i=0; i < N; ++i) {
00732             i3 = 3*i;
00733             old_fv.xyz[0] = xyzv[i3];
00734             old_fv.xyz[1] = xyzv[i3+1];
00735             old_fv.xyz[2] = xyzv[i3+2];
00736             
00737             for( iw=0 ; iw < 12 ; iw++ ){
00738                map    = warp->tal_12.warp[iw] ;
00739                new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00740 
00741                if( new_fv.xyz[0] >= map.bot.xyz[0] &&
00742                    new_fv.xyz[1] >= map.bot.xyz[1] &&
00743                    new_fv.xyz[2] >= map.bot.xyz[2] &&
00744                    new_fv.xyz[0] <= map.top.xyz[0] &&
00745                    new_fv.xyz[1] <= map.top.xyz[1] &&
00746                    new_fv.xyz[2] <= map.top.xyz[2]   ) break ;  /* leave loop */
00747             }
00748             xyzv[i3] = new_fv.xyz[0];
00749             xyzv[i3+1] = new_fv.xyz[1];
00750             xyzv[i3+2] = new_fv.xyz[2];
00751             
00752          }
00753       }
00754       break ;
00755 
00756       case WARP_AFFINE_TYPE:{
00757          THD_linear_mapping map = warp->rig_bod.warp ;
00758          for (i=0; i < N; ++i) {
00759             i3 = 3*i;
00760             old_fv.xyz[0] = xyzv[i3];
00761             old_fv.xyz[1] = xyzv[i3+1];
00762             old_fv.xyz[2] = xyzv[i3+2];
00763             new_fv = MATVEC_SUB(map.mfor,old_fv,map.bvec) ;
00764             xyzv[i3] = new_fv.xyz[0];
00765             xyzv[i3+1] = new_fv.xyz[1];
00766             xyzv[i3+2] = new_fv.xyz[2];
00767          }
00768       }
00769       break ;
00770 
00771    }
00772    SUMA_RETURN(YUP);
00773 }