00001 #include "mrilib.h"
00002 
00003 /********************************************************************
00004  ****** Functions to create and deal with SUMA_surface structs ******
00005  ********************************************************************/
00006 
00007 #define SUMA_EXTEND_NUM 64
00008 #define SUMA_EXTEND_FAC 1.05
00009 
00010 /*------------------------------------------------------------------*/
00011 /*! Create an empty surface description.
00012 --------------------------------------------------------------------*/
00013 
00014 SUMA_surface * SUMA_create_empty_surface(void)
00015 {
00016    SUMA_surface *ag ;
00017 
00018 ENTRY("SUMA_create_empty_surface") ;
00019 
00020    ag       = (SUMA_surface *) calloc(1,sizeof(SUMA_surface)) ;
00021    ag->type = SUMA_SURFACE_TYPE ;
00022 
00023    ag->num_ixyz  = ag->num_ijk  = 0 ;
00024    ag->nall_ixyz = ag->nall_ijk = 1 ;
00025    ag->ixyz = (SUMA_ixyz *) malloc(sizeof(SUMA_ixyz)) ; /* space for */
00026    ag->ijk  = (SUMA_ijk *)  malloc(sizeof(SUMA_ijk) ) ; /* 1 of each */
00027    ag->norm = NULL ;                                  ; /* none of this */
00028 
00029    if( ag->ixyz == NULL || ag->ijk == NULL ){
00030       fprintf(stderr,"SUMA_create_empty_surface: can't malloc!\n"); EXIT(1);
00031    }
00032 
00033    ag->idcode[0] = ag->idcode_dset[0] = ag->idcode_ldp[0] =
00034       ag->label[0] = ag->label_ldp[0] = '\0' ;
00035 
00036    ag->xbot = ag->ybot = ag->zbot =  WAY_BIG ;
00037    ag->xtop = ag->ytop = ag->ztop = -WAY_BIG ;
00038    ag->xcen = ag->ycen = ag->zcen = 0.0      ;
00039 
00040    ag->seq = ag->seqbase = ag->sorted = 0 ; /* not sequential; not sorted */
00041 
00042    ag->vv = NULL ;  /* 16 Jun 2003 */
00043    ag->vn = NULL ;  /* 22 Jan 2004 */
00044 
00045    RETURN( ag ) ;
00046 }
00047 
00048 /*------------------------------------------------------------------*/
00049 /*! Throw out some trash (i.e., free the contents of a surface).
00050 --------------------------------------------------------------------*/
00051 
00052 void SUMA_destroy_surface( SUMA_surface *ag )
00053 {
00054 ENTRY("SUMA_destroy_surface") ;
00055 
00056    if( ag == NULL ) EXRETURN ;
00057    if( ag->ixyz != NULL ) free((void *)ag->ixyz) ;
00058    if( ag->ijk  != NULL ) free((void *)ag->ijk) ;
00059    if( ag->norm != NULL ) free((void *)ag->norm) ;
00060 
00061    if( ag->vv != NULL ) DESTROY_VVLIST(ag->vv) ;
00062    if( ag->vn != NULL ) SUMA_destroy_vnlist(ag->vn) ;
00063 
00064    free((void *)ag) ; EXRETURN ;
00065 }
00066 
00067 /*------------------------------------------------------------------*/
00068 /* Add a bunch of nodes to a surface.
00069 --------------------------------------------------------------------*/
00070 
00071 void SUMA_add_nodes_ixyz( SUMA_surface *ag, int nadd,
00072                           int *iadd, float *xadd, float *yadd, float *zadd )
00073 {
00074    int ii , nup ;
00075 
00076 ENTRY("SUMA_add_nodes_ixyz") ;
00077 
00078    if( ag == NULL || nadd < 1 ) EXRETURN ;
00079    if( xadd == NULL || yadd == NULL || zadd == NULL || iadd == NULL ) EXRETURN ;
00080 
00081    nup = ag->num_ixyz + nadd ;
00082 
00083    if( nup >= SUMA_MAX_NODES ){  /* 07 Sep 2001 */
00084      fprintf(stderr,
00085              "** SUMA surface can't have more than %d nodes!\n",
00086              SUMA_MAX_NODES-1 ) ;
00087      EXRETURN ;
00088    }
00089 
00090    if( nup > ag->nall_ixyz ){ /* extend length of array */
00091      ag->nall_ixyz = nup = nup*SUMA_EXTEND_FAC + SUMA_EXTEND_NUM ;
00092      ag->ixyz = (SUMA_ixyz *) realloc( (void *)ag->ixyz, sizeof(SUMA_ixyz)*nup );
00093      if( ag->ixyz == NULL ){
00094        fprintf(stderr,"SUMA_add_nodes_ixyz: can't malloc!\n"); EXIT(1);
00095      }
00096    }
00097 
00098    nup = ag->num_ixyz ;
00099 
00100    for( ii=0 ; ii < nadd ; ii++ ){
00101      ag->ixyz[ii+nup].x  = xadd[ii] ;
00102      ag->ixyz[ii+nup].y  = yadd[ii] ;
00103      ag->ixyz[ii+nup].z  = zadd[ii] ;
00104      ag->ixyz[ii+nup].id = iadd[ii] ;
00105    }
00106 
00107    ag->num_ixyz += nadd ;
00108 
00109    ag->seq = ag->sorted = 0 ; EXRETURN ;
00110 }
00111 
00112 /*--------------------------------------------------------------------
00113  * Add/replace normals on the given surface.       05 Oct 2004 [rickr]
00114  *
00115  * This function requires one normal per node, and that the
00116  * indices match.
00117  *--------------------------------------------------------------------
00118 */
00119 int SUMA_add_norms_xyz( SUMA_surface *ag, int nadd,
00120                         float *xadd, float *yadd, float *zadd )
00121 {
00122    int ii ;
00123 
00124 ENTRY("SUMA_add_norms_xyz") ;
00125 
00126    if( ag == NULL || nadd < 1 ) RETURN(-1) ;
00127    if( xadd == NULL || yadd == NULL || zadd == NULL ) RETURN(-1) ;
00128 
00129    if( nadd != ag->num_ixyz ){
00130      fprintf(stderr, "** SUMA surface has %d nodes but %d normals!\n",
00131              ag->num_ixyz, nadd ) ;
00132      RETURN(-1) ;
00133    }
00134 
00135    /* if norm is NULL, memory is needed */
00136    if( ag->norm == NULL ){
00137        ag->norm = (THD_fvec3 *)calloc(nadd, sizeof(THD_fvec3));
00138        if( ag->norm == NULL ){
00139            fprintf(stderr,"SUMA_add_norms_xyz: can't malloc!\n"); EXIT(1);
00140        }
00141    }
00142 
00143    for( ii=0 ; ii < nadd ; ii++ ){
00144      ag->norm[ii].xyz[0] = xadd[ii] ;
00145      ag->norm[ii].xyz[1] = yadd[ii] ;
00146      ag->norm[ii].xyz[2] = zadd[ii] ;
00147    }
00148 
00149    RETURN(0) ;
00150 }
00151 
00152 /*------------------------------------------------------------------*/
00153 /*! Add 1 pitiful node to a surface.
00154 --------------------------------------------------------------------*/
00155 
00156 void SUMA_add_node_ixyz( SUMA_surface *ag , int i,float x,float y,float z )
00157 {
00158    SUMA_add_nodes_ixyz( ag , 1 , &i,&x,&y,&z ) ;
00159 }
00160 
00161 /*------------------------------------------------------------------*/
00162 /*!  Add a bunch of triangles (node id triples) to a surface.
00163 --------------------------------------------------------------------*/
00164 
00165 void SUMA_add_triangles( SUMA_surface *ag, int nadd, int *it, int *jt, int *kt )
00166 {
00167    int ii , nup ;
00168 
00169 ENTRY("SUMA_add_triangles") ;
00170 
00171    if( ag == NULL || nadd < 1 ) EXRETURN ;
00172    if( it == NULL || jt == NULL || kt == NULL ) EXRETURN ;
00173 
00174    nup = ag->num_ijk + nadd ;
00175    if( nup > ag->nall_ijk ){ /* extend length of array */
00176      ag->nall_ijk = nup = nup*SUMA_EXTEND_FAC + SUMA_EXTEND_NUM ;
00177      ag->ijk = (SUMA_ijk *) realloc( (void *)ag->ijk , sizeof(SUMA_ijk)*nup ) ;
00178      if( ag->ijk == NULL ){
00179        fprintf(stderr,"SUMA_add_triangles: can't malloc!\n"); EXIT(1);
00180      }
00181    }
00182 
00183    nup = ag->num_ijk ;
00184    for( ii=0 ; ii < nadd ; ii++ ){
00185      ag->ijk[ii+nup].id = it[ii] ;
00186      ag->ijk[ii+nup].jd = jt[ii] ;
00187      ag->ijk[ii+nup].kd = kt[ii] ;
00188    }
00189 
00190    ag->num_ijk += nadd ; EXRETURN ;
00191 }
00192 
00193 /*------------------------------------------------------------------*/
00194 /*! Add 1 pitiful triangle to a surface.
00195 --------------------------------------------------------------------*/
00196 
00197 void SUMA_add_triangle( SUMA_surface *ag, int it, int jt, int kt )
00198 {
00199    SUMA_add_triangles( ag , 1 , &it,&jt,&kt ) ;
00200 }
00201 
00202 /*------------------------------------------------------------------*/
00203 /*! Truncate the memory used by the node and triangle arrays back
00204     to the minimum they need.
00205 --------------------------------------------------------------------*/
00206 
00207 void SUMA_truncate_memory( SUMA_surface *ag )
00208 {
00209    int nn ;
00210 
00211 ENTRY("SUMA_truncate_memory") ;
00212 
00213    if( ag == NULL ) EXRETURN ;
00214 
00215    if( ag->num_ixyz < ag->nall_ixyz && ag->num_ixyz > 0 ){
00216      ag->nall_ixyz = nn = ag->num_ixyz ;
00217      ag->ixyz = (SUMA_ixyz *) realloc( (void *)ag->ixyz, sizeof(SUMA_ixyz)*nn );
00218    }
00219 
00220    if( ag->num_ijk < ag->nall_ijk && ag->num_ijk > 0 ){
00221      ag->nall_ijk = nn = ag->num_ijk ;
00222      ag->ijk = (SUMA_ijk *) realloc( (void *)ag->ijk , sizeof(SUMA_ijk)*nn ) ;
00223    }
00224 
00225    EXRETURN ;
00226 }
00227 
00228 /*------------------------------------------------------------------
00229   Generate a function to sort array of SUMA_ixyz-s by their id-s
00230 --------------------------------------------------------------------*/
00231 
00232 #undef  STYPE
00233 #define STYPE     SUMA_ixyz          /* name of type to sort     */
00234 #define SLT(a,b)  ((a).id < (b).id)  /* macro to decide on order */
00235 #define SNAME     SUMA_ixyz          /* function qsort_SUMA_ixyz */
00236 #include "cs_sort_template.h"        /* generate the function    */
00237 #undef  STYPE
00238 
00239 /*** Can now use qsort_SUMA_ixyz(int,SUMA_ixyz *) ***/
00240 
00241 /*------------------------------------------------------------------*/
00242 /*!  Sort the nodes by id-s, and mark if the id-s are sequential.
00243 --------------------------------------------------------------------*/
00244 
00245 void SUMA_ixyzsort_surface( SUMA_surface *ag )
00246 {
00247    int nn , ii , ndup ;
00248    float xb,yb,zb , xt,yt,zt , xc,yc,zc ;
00249 
00250 ENTRY("SUMA_ixyzsort_surface") ;
00251 
00252    if( ag == NULL || ag->num_ixyz < 1 ) EXRETURN ;
00253 
00254    SUMA_truncate_memory( ag ) ;
00255 
00256    nn = ag->num_ixyz ;
00257 
00258    /* check if nodes are already sorted [26 Oct 2001] */
00259 
00260    for( ii=1 ; ii < nn ; ii++ )
00261      if( ag->ixyz[ii].id <= ag->ixyz[ii-1].id ) break ;
00262 
00263    /* if not in increasing order,
00264       sort them using the function generated above */
00265 
00266    if( ii < nn ){
00267      qsort_SUMA_ixyz( nn , ag->ixyz ) ;
00268    }
00269 
00270    ag->sorted = 1 ;  /* mark as sorted */
00271 
00272    /* check if node id-s are sequential */
00273 
00274    for( ii=1 ; ii < nn ; ii++ )
00275      if( ag->ixyz[ii].id != ag->ixyz[ii-1].id+1 ) break ;
00276 
00277    /* if we finished that loop all the way,
00278       mark the nodes as being sequential, and
00279       store the base of the sequence (id of node #0) */
00280 
00281    if( ii == nn ){
00282      ag->seq = 1 ; ag->seqbase = ag->ixyz[0].id ;
00283    }
00284 
00285    /* 07 Sep 2001: check for duplicate node id-s */
00286 
00287    for( ndup=0,ii=1 ; ii < nn ; ii++ )
00288      if( ag->ixyz[ii].id == ag->ixyz[ii-1].id ) ndup++ ;
00289 
00290    if( ndup > 0 )
00291      fprintf(stderr,"** SUMA WARNING: %d duplicate surface node id's found!\n",ndup);
00292 
00293    /* find bounding box of all nodes (it's useful on occasion) */
00294 
00295    xb = xt = ag->ixyz[0].x ;
00296    yb = yt = ag->ixyz[0].y ;
00297    zb = zt = ag->ixyz[0].z ;
00298    xc = yc = zc = 0.0 ;
00299    for( ii=1 ; ii < nn ; ii++ ){
00300      xc += ag->ixyz[ii].x ;
00301      yc += ag->ixyz[ii].y ;
00302      zc += ag->ixyz[ii].z ;
00303 
00304           if( ag->ixyz[ii].x < xb ) xb = ag->ixyz[ii].x ;
00305      else if( ag->ixyz[ii].x > xt ) xt = ag->ixyz[ii].x ;
00306 
00307           if( ag->ixyz[ii].y < yb ) yb = ag->ixyz[ii].y ;
00308      else if( ag->ixyz[ii].y > yt ) yt = ag->ixyz[ii].y ;
00309 
00310           if( ag->ixyz[ii].z < zb ) zb = ag->ixyz[ii].z ;
00311      else if( ag->ixyz[ii].z > zt ) zt = ag->ixyz[ii].z ;
00312    }
00313 
00314    ag->xbot = xb ; ag->xtop = xt ;
00315    ag->ybot = yb ; ag->ytop = yt ;
00316    ag->zbot = zb ; ag->ztop = zt ;
00317 
00318    ag->xcen = xc/nn ; ag->ycen = yc/nn ; ag->zcen = zc/nn ;
00319 
00320    EXRETURN ;
00321 }
00322 
00323 /*--------------------------------------------------------------------*/
00324 /*! Find a node id in a surface, and return its index into the node
00325     array; return -1 if not found.
00326 ----------------------------------------------------------------------*/
00327 
00328 int SUMA_find_node_id( SUMA_surface *ag , int target )
00329 {
00330    int nn , ii,jj,kk ;
00331 
00332    if( ag == NULL || ag->num_ixyz < 1 || target < 0 ) return( -1 );
00333 
00334    if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
00335 
00336    if( ag->seq ){  /* node id-s are sequential (the easy case) */
00337       kk = target - ag->seqbase ;
00338       if( kk >= 0 && kk < ag->num_ixyz ) return( kk );
00339       return( -1 );
00340    }
00341 
00342    /* node id-s are in increasing order, but not sequential;
00343       so, use binary search to find the node id (if present) */
00344 
00345    ii = 0 ; jj = ag->num_ixyz - 1 ;                 /* search bounds */
00346 
00347         if( target <  ag->ixyz[0].id  ) return( -1 ); /* not present */
00348    else if( target == ag->ixyz[0].id  ) return( ii ); /* at start!  */
00349 
00350         if( target >  ag->ixyz[jj].id ) return( -1 ); /* not present */
00351    else if( target == ag->ixyz[jj].id ) return( jj ); /* at end!    */
00352 
00353    while( jj - ii > 1 ){  /* while search bounds not too close */
00354 
00355       kk = (ii+jj) / 2 ;  /* midway between search bounds */
00356 
00357       nn = ag->ixyz[kk].id - target ;
00358       if( nn == 0 ) return( kk );     /* AHA! */
00359 
00360       if( nn < 0 ) ii = kk ;          /* kk before target => bottom = kk */
00361       else         jj = kk ;          /* kk after target  => top    = kk */
00362    }
00363 
00364    return( -1 );
00365 }
00366 
00367 /*-------------------------------------------------------------------------*/
00368 /*! Create the voxel-to-node list for this surface/dataset combo.
00369 ---------------------------------------------------------------------------*/
00370 
00371 SUMA_vnlist * SUMA_make_vnlist( SUMA_surface *ag , THD_3dim_dataset *dset )
00372 {
00373    int ii,jj,kk , nx,ny,nz , nxy,nxyz , nnode , pp,qq,nn,nvox  ;
00374    THD_fvec3 fv ;
00375    THD_ivec3 iv ;
00376    int *vlist , *nlist , wodsave ;
00377    SUMA_vnlist *vnlist ;
00378    float xbot,xtop , ybot,ytop , zbot,ztop ;
00379 
00380 ENTRY("SUMA_make_vnlist") ;
00381 
00382    if( ag == NULL || ag->num_ixyz < 1 || !ISVALID_DSET(dset) ) RETURN(NULL) ;
00383 
00384    if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
00385 
00386    /* setup: create arrays for voxel list and node list */
00387 
00388    nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
00389    nxy = nx*ny ; nxyz = nxy*nz ; nnode = ag->num_ixyz ;
00390    vlist = (int *) malloc(sizeof(int)*nnode) ;
00391    nlist = (int *) malloc(sizeof(int)*nnode) ;
00392    if( vlist == NULL || nlist == NULL ){
00393       fprintf(stderr,"SUMA_make_vnlist: can't malloc!\n"); EXIT(1);
00394    }
00395 
00396    /* for each node, find which voxel it is in */
00397 
00398    wodsave = dset->wod_flag ; dset->wod_flag = 0 ;
00399 
00400    xbot = DSET_XXMIN(dset) ; xtop = DSET_XXMAX(dset) ;
00401    ybot = DSET_YYMIN(dset) ; ytop = DSET_YYMAX(dset) ;
00402    zbot = DSET_ZZMIN(dset) ; ztop = DSET_ZZMAX(dset) ;
00403 
00404    for( nn=pp=0 ; pp < nnode ; pp++ ){
00405       LOAD_FVEC3( fv , ag->ixyz[pp].x, ag->ixyz[pp].y, ag->ixyz[pp].z ) ;
00406       fv = THD_dicomm_to_3dmm( dset , fv ) ; /* convert Dicom coords */
00407 
00408       if( fv.xyz[0] < xbot || fv.xyz[0] > xtop ) continue ;
00409       if( fv.xyz[1] < ybot || fv.xyz[1] > ytop ) continue ;
00410       if( fv.xyz[2] < zbot || fv.xyz[2] > ztop ) continue ;
00411 
00412       iv = THD_3dmm_to_3dind( dset , fv ) ;  /*   in surface to     */
00413       UNLOAD_IVEC3( iv , ii,jj,kk ) ;        /*   dataset indexes  */
00414 
00415       nlist[nn] = pp ;                       /* list of nodes */
00416       vlist[nn] = ii + jj*nx + kk*nxy ;      /* list of voxels */
00417       nn++ ;
00418    }
00419 
00420    nnode = nn ; /* number of nodes inside dataset volume */
00421    if( nnode == 0 ){ free(nlist); free(vlist); RETURN(NULL); }
00422 
00423    dset->wod_flag = wodsave ;
00424 
00425    /* now sort the 2 lists so that vlist is increasing
00426       (and each nlist still corresponds to its original vlist) */
00427 
00428    qsort_intint( nnode , vlist , nlist ) ;
00429 
00430    /* count how many distinct voxels we found */
00431 
00432    nvox = 1 ; ii = vlist[0] ;
00433    for( pp=1 ; pp < nnode ; pp++ ){
00434       if( vlist[pp] != ii ){ nvox++; ii = vlist[pp]; }
00435    }
00436 
00437    /* now create the output vnlist */
00438 
00439    vnlist         = (SUMA_vnlist *) malloc( sizeof(SUMA_vnlist) ) ;
00440    vnlist->nvox   = nvox ;
00441    vnlist->voxijk = (int *) malloc(sizeof(int) *nvox) ;
00442    vnlist->numnod = (int *) calloc(sizeof(int) ,nvox) ;
00443    vnlist->nlist  = (int **)malloc(sizeof(int*)*nvox);
00444    vnlist->dset   = dset ;
00445 
00446    if( vnlist->voxijk==NULL || vnlist->numnod==NULL || vnlist->nlist==NULL ){
00447      fprintf(stderr,"SUMA_make_vnlist: can't malloc!\n"); EXIT(1);
00448    }
00449 
00450    /* now count how many nodes are at each voxel in the list */
00451 
00452    ii = vlist[0] ; qq = nn = 0 ;
00453    for( pp=1 ; pp < nnode ; pp++ ){
00454      if( vlist[pp] != ii ){         /* qq..pp-1 are the same */
00455        vnlist->voxijk[nn] = ii ;
00456        vnlist->numnod[nn] = jj = pp-qq ;
00457        vnlist->nlist[nn]  = (int *) malloc(sizeof(int)*jj) ;
00458        memcpy( vnlist->nlist[nn] , nlist+qq , sizeof(int)*jj ) ;
00459        ii = vlist[pp] ; nn++ ; qq = pp ;
00460      }
00461    }
00462    vnlist->voxijk[nn] = ii ;
00463    vnlist->numnod[nn] = jj = pp-qq ;
00464    vnlist->nlist[nn]  = (int *) malloc(sizeof(int)*jj) ;
00465    memcpy( vnlist->nlist[nn] , nlist+qq , sizeof(int)*jj ) ;
00466 
00467    /* and we're done! */
00468 
00469    free(nlist) ; free(vlist) ; RETURN( vnlist ) ;
00470 }
00471 
00472 /*-------------------------------------------------------------------------*/
00473 /*! Destroy a SUMA_vnlist struct.
00474 ---------------------------------------------------------------------------*/
00475 
00476 void SUMA_destroy_vnlist( SUMA_vnlist *vnlist )
00477 {
00478    int ii ;
00479    if( vnlist == NULL ) return ;
00480    if( vnlist->voxijk != NULL ) free( vnlist->voxijk ) ;
00481    if( vnlist->numnod != NULL ) free( vnlist->numnod ) ;
00482    if( vnlist->nlist  != NULL ){
00483      for( ii=0 ; ii < vnlist->nvox ; ii++ )
00484        if( vnlist->nlist[ii] != NULL ) free( vnlist->nlist[ii] ) ;
00485      free( vnlist->nlist ) ;
00486    }
00487    free( vnlist ) ;
00488 }
00489 
00490 /*--------------------------------------------------------------------------
00491    The following routines are used to convert DICOM order coordinates
00492    (used in AFNI) to SureFit order coordinates -- 25 Oct 2001 - RWCox
00493 ----------------------------------------------------------------------------*/
00494 
00495 THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *dset , THD_fvec3 fv )
00496 {
00497    float xx,yy,zz , xbase,ybase,zbase ;
00498    THD_fvec3 vout ;
00499 
00500    xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;   /* xyz now LPI */
00501 
00502    if( dset != NULL ){
00503       THD_fvec3 v1 , v2 ;
00504       LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
00505       v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
00506       LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
00507                      DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
00508                      DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset)  ) ;
00509       v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
00510       xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;  /* Left-most */
00511       ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;  /* Posterior */
00512       zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;                   /* Inferior  */
00513    } else {
00514       xbase = ybase = zbase = 0.0 ;
00515    }
00516 
00517    vout.xyz[0] = xx - xbase ;
00518    vout.xyz[1] = yy - ybase ;
00519    vout.xyz[2] = zz - zbase ; return vout ;
00520 }
00521 
00522 /* --------------------------------------------------------------------------*/
00523 
00524 THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 fv )
00525 {
00526    float xx,yy,zz , xbase,ybase,zbase ;
00527    THD_fvec3 vout ;
00528 
00529    xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ;   /* xyz now RAI */
00530 
00531    if( dset != NULL ){
00532       THD_fvec3 v1 , v2 ;
00533       LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
00534       v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
00535       LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
00536                      DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
00537                      DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset)  ) ;
00538       v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
00539       xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
00540       ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
00541       zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
00542    } else {
00543       xbase = ybase = zbase = 0.0 ;
00544    }
00545 
00546    vout.xyz[0] = xx - xbase ;
00547    vout.xyz[1] = yy - ybase ;
00548    vout.xyz[2] = zz + zbase ; return vout ;
00549 }
00550 
00551 /****************************************************************************
00552  ********** AFNI no longer reads surface from files [22 Jan 2004] ***********
00553  ****************************************************************************/
00554 
00555 #undef ALLOW_SURFACE_FILES
00556 #ifdef ALLOW_SURFACE_FILES
00557 /*----------------------------------------------------------*/
00558 /*! Will load this many items (nodes, triangles) at a time. */
00559 
00560 #undef  NBUF
00561 #define NBUF 64
00562 
00563 /*------------------------------------------------------------------------*/
00564 /*! Read a surface description file that is associated with an AFNI
00565     dataset.  Return a surface ready to rock-n-roll.
00566     NULL is returned if something bad happens.
00567 --------------------------------------------------------------------------*/
00568 
00569 SUMA_surface * SUMA_read_surface( char *fname , THD_3dim_dataset *dset )
00570 {
00571    SUMA_surface *ag ;
00572    FILE *fp ;
00573    char lbuf[1024] , *cpt ;
00574    int  do_nod=1 , ii ;
00575    float xx[NBUF],yy[NBUF],zz[NBUF] ;
00576    int   pp[NBUF],qq[NBUF],rr[NBUF] , nn ;
00577 
00578    THD_vecmat mv ;
00579    int have_mv=0 ;
00580 
00581 ENTRY("SUMA_read_surface") ;
00582 
00583    if( fname == NULL || fname[0] == '\0' ) RETURN( NULL );
00584 
00585    /*-- open input --*/
00586 
00587    if( strcmp(fname,"-") == 0 ){   /* special case */
00588       fp = stdin ;
00589    } else {
00590       fp = fopen( fname , "r" ) ;
00591       if( fp == NULL ) RETURN( NULL );
00592    }
00593 
00594    /*-- initialize surface that we will fill up here */
00595 
00596    ag = SUMA_create_empty_surface() ;
00597 
00598    /*-- set IDCODEs of surface (from filename) and of its dataset --*/
00599 
00600    cpt = UNIQ_hashcode(fname); strcpy(ag->idcode,cpt); free(cpt);
00601 
00602    strcpy( ag->idcode_dset , dset->idcode.str ) ;
00603 
00604    MCW_strncpy( ag->label, THD_trailname(fname,0), 32 ) ;  /* 19 Aug 2002 */
00605 
00606    /*-- read data --*/
00607 
00608    nn = 0 ;
00609    while(1){
00610       cpt = fgets(lbuf,1024,fp) ;  /* read a line */
00611       if( cpt == NULL ) break ;    /* end of file */
00612 
00613       /*-- read a transformation matrix-vector? --*/
00614 
00615       if( strncmp(lbuf,"<MATVEC>",8) == 0 ){  /* 07 Sep 2001 */
00616          float a11,a12,a13 , v1 ,
00617                a21,a22,a23 , v2 ,
00618                a31,a32,a33 , v3  ;
00619 
00620          ii = sscanf( lbuf+8 , "%f%f%f%f%f%f%f%f%f%f%f%f" ,
00621                       &a11,&a12,&a13 , &v1 ,
00622                       &a21,&a22,&a23 , &v2 ,
00623                       &a31,&a32,&a33 , &v3  ) ;
00624 
00625          if( ii < 12 ){
00626             fprintf(stderr,"** SUMA: Illegal MATVEC in %s\n",fname) ;
00627             have_mv = 0 ;
00628          } else {
00629             LOAD_FVEC3(mv.vv , v1,v2,v3 ) ;
00630             LOAD_MAT  (mv.mm , a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
00631             have_mv = 1 ;
00632          }
00633          continue ; /* skip to next line */
00634       }
00635 
00636       /*-- read data from SureFit? --*/
00637 
00638       if( strncmp(lbuf,"<SureFit",8) == 0 ){ /* 26 Oct 2001 */
00639 
00640          if( nn > 0 ){   /* process existing inputs, if any */
00641             if( do_nod )
00642                SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
00643             else
00644                SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
00645             nn = 0 ;
00646          }
00647 
00648          SUMA_import_surefit( ag , lbuf , dset ) ;
00649          continue ; /* skip to next input line */
00650 
00651       } /* end of SureFit input */
00652 
00653       /*-- end of node input? --*/
00654 
00655       if( strstr(lbuf,"</NODES>") != NULL ){
00656          if( do_nod && nn > 0 ){
00657             SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
00658             nn = 0 ;
00659          }
00660 #if 1
00661          do_nod = 0 ;  /* from now on, triangles */
00662          continue ;    /* skip to next line */
00663 #else
00664          break ;                  /* don't create triangles */
00665 #endif
00666       }
00667 
00668       /*-- process line as a node? --*/
00669 
00670       if( do_nod ){               /* nodes */
00671 
00672          ii = sscanf(lbuf,"%d%f%f%f",pp+nn,xx+nn,yy+nn,zz+nn) ;
00673          if( ii < 4 ) continue ;
00674          if( have_mv ){           /* transform coords */
00675             THD_fvec3 fv , qv ;
00676             LOAD_FVEC3(fv , xx[nn],yy[nn],zz[nn] ) ;
00677             qv = VECSUB_MAT( mv.mm , fv , mv.vv ) ;
00678             UNLOAD_FVEC3( qv , xx[nn],yy[nn],zz[nn] ) ;
00679          }
00680          nn++ ;
00681          if( nn == NBUF ){        /* add nodes in chunks */
00682             SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
00683             nn = 0 ;
00684          }
00685 
00686       /*-- process line as a triangle --*/
00687 
00688       } else {                    /* triangles */
00689 
00690          ii = sscanf(lbuf,"%d%d%d",pp+nn,qq+nn,rr+nn) ;
00691          if( ii < 3 ) continue ;
00692          nn++ ;
00693          if( nn == NBUF ){        /* add triangles in chunks */
00694             SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
00695             nn = 0 ;
00696          }
00697       }
00698    } /* end of loop over input lines */
00699 
00700    /*-- finish up, eh? --*/
00701 
00702    if( fp != stdin ) fclose(fp) ;
00703    if( nn > 0 ){                  /* process leftover data lines */
00704       if( do_nod )
00705          SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
00706       else
00707          SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
00708    }
00709 
00710    /*-- hope we got something --*/
00711 
00712    if( ag->num_ixyz < 1 ){
00713       SUMA_destroy_surface(ag) ; RETURN(NULL) ;
00714    }
00715 
00716    /*-- sort the nodes (if needed), et cetera --*/
00717 
00718    SUMA_ixyzsort_surface(ag) ;
00719 
00720    /*-- done --*/
00721 
00722    RETURN( ag );
00723 }
00724 
00725 /*-----------------------------------------------------------------------------*/
00726 
00727 /*! 26 point 3x3x3 nbhd of {0,0,0}
00728     (not including the central point itself) */
00729 
00730 static int ip[26][3] = { {-1,-1,-1},{-1,-1, 0},{-1,-1, 1},
00731                          {-1, 0,-1},{-1, 0, 0},{-1, 0, 1},
00732                          {-1, 1,-1},{-1, 1, 0},{-1, 1, 1},
00733                          { 0,-1,-1},{ 0,-1, 0},{ 0,-1, 1},
00734                          { 0, 0,-1},           { 0, 0, 1},
00735                          { 0, 1,-1},{ 0, 1, 0},{ 0, 1, 1},
00736                          { 1,-1,-1},{ 1,-1, 0},{ 1,-1, 1},
00737                          { 1, 0,-1},{ 1, 0, 0},{ 1, 0, 1},
00738                          { 1, 1,-1},{ 1, 1, 0},{ 1, 1, 1} } ;
00739 
00740 /*------------------------------------------------------------------------*/
00741 /*! Given a surface and a dataset, compute a voxel-to-node map.
00742     This is an int array (the return value) the size of a dataset
00743     brick - one value per voxel.  Each entry "v" is
00744      - v < 0 means that voxel is not close to a surface node
00745      - otherwise, the closest node (index in ag->ixyz, NOT id) to the
00746         voxel is SUMA_VMAP_UNMASK(v)
00747      - the "level" of that node is SUMA_VMAP_LEVEL(v),
00748         where level is the number of nbhd expansions outward from the
00749         voxel that were needed to hit a node (0<=level<=7).
00750      - The node index (26 bit integer) is stored in bits 0..25 of v
00751      - The level (3 bit integer) is stored in bits 26..28 of v
00752      - Bits 29..30 are reserved for future purposes
00753      - Bit 31 is the sign bit, and so signifies "no node nearby"
00754 
00755     Other notes:
00756      - The input surface should have been sorted by SUMA_ixyzsort_surface().
00757         If it wasn't, it will be now.
00758      - Although this function was crafted to be efficient, it still takes
00759         a few seconds to run.
00760      - Don't add nodes to the surface after calling this, since they
00761         won't be indexed here properly.  Or you could free() the output
00762         map and re-invoke this function.
00763 ---------------------------------------------------------------------------*/
00764 
00765 int * SUMA_map_dset_to_surf( SUMA_surface *ag , THD_3dim_dataset *dset )
00766 {
00767    int *vmap , ii,jj,kk , nx,ny,nz , nxy,nxyz , pp,qq,pbest ;
00768    THD_fvec3 fv ;
00769    THD_ivec3 iv ;
00770    int ibot,jbot,kbot , itop,jtop,ktop , lev , ijk ;
00771    float xv,yv,zv , dd,dbest=0 , xp,yp,zp ;
00772    char *elev ; int ltop , ntop , lmask ;
00773 
00774 ENTRY("SUMA_map_dset_to_surf") ;
00775 
00776    if( ag == NULL || ag->num_ixyz < 1 || !ISVALID_DSET(dset) ) RETURN( NULL );
00777 
00778    if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
00779 
00780    /* setup: create the output vmap and fill it with -1 */
00781 
00782    nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
00783    nxy = nx*ny ; nxyz = nxy*nz ;
00784    vmap = (int *) malloc(sizeof(int)*nxyz) ;
00785    if( vmap == NULL ){
00786       fprintf(stderr,"SUMA_map_dset_to_surf: can't malloc!\n"); EXIT(1);
00787    }
00788    for( ii=0 ; ii < nxyz ; ii++ ) vmap[ii] = -1 ; /* not mapped yet */
00789 
00790    /* put nodes directly into voxels (level 0) */
00791 
00792 STATUS("putting nodes into voxels") ;
00793 
00794    for( pp=0 ; pp < ag->num_ixyz ; pp++ ){
00795       LOAD_FVEC3( fv , ag->ixyz[pp].x, ag->ixyz[pp].y, ag->ixyz[pp].z ) ;
00796       fv = THD_dicomm_to_3dmm( dset , fv ) ; /* convert Dicom coords */
00797       iv = THD_3dmm_to_3dind( dset , fv ) ;  /*   in surface to     */
00798       UNLOAD_IVEC3( iv , ii,jj,kk ) ;        /*   dataset indexes  */
00799       qq = vmap[ii+jj*nx+kk*nxy] ;           /* previously mapped index? */
00800       if( qq < 0 ){                          /* not mapped before */
00801          vmap[ii+jj*nx+kk*nxy] = pp ;        /* index, not id */
00802       } else {
00803          LOAD_IVEC3(iv,ii,jj,kk) ;           /* get Dicom coords of voxel */
00804          fv = THD_3dind_to_3dmm( dset , iv ) ;
00805          fv = THD_3dmm_to_dicomm( dset , fv ) ;
00806          UNLOAD_FVEC3(fv,xv,yv,zv) ;         /* voxel is at (xv,yv,zv) */
00807 
00808          /* find if old node in this voxel (#qq)
00809             is closer to voxel center than current node (#pp) */
00810 
00811          xp=xv-ag->ixyz[qq].x;
00812          yp=yv-ag->ixyz[qq].y;
00813          zp=zv-ag->ixyz[qq].z; dd=xp*xp+yp*yp+zp*zp;    /* dist^2 to old node */
00814 
00815          xp=xv-ag->ixyz[pp].x;
00816          yp=yv-ag->ixyz[pp].y;
00817          zp=zv-ag->ixyz[pp].z; dbest=xp*xp+yp*yp+zp*zp; /* dist^2 to new node */
00818 
00819          if( dbest < dd ) vmap[ii+jj*nx+kk*nxy] = pp ;   /* new is better */
00820       }
00821    }
00822 
00823    LOAD_FVEC3( fv , ag->xbot,ag->ybot,ag->zbot ) ; /* find dataset */
00824    fv = THD_dicomm_to_3dmm( dset , fv ) ;          /* indexes for */
00825    iv = THD_3dmm_to_3dind( dset , fv ) ;           /* bot point  */
00826    UNLOAD_IVEC3( iv , ibot,jbot,kbot ) ;           /* of surface */
00827 
00828    LOAD_FVEC3( fv , ag->xtop,ag->ytop,ag->ztop ) ; /* and for top */
00829    fv = THD_dicomm_to_3dmm( dset , fv ) ;
00830    iv = THD_3dmm_to_3dind( dset , fv ) ;
00831    UNLOAD_IVEC3( iv , itop,jtop,ktop ) ;
00832 
00833    /* we will only try to fill voxels inside
00834       the rectangular subvolume (ibot..itop,jbot..jtop,kbot..ktop) */
00835 
00836    if( ibot > itop ){ ii=ibot ; ibot=itop; itop=ii; }
00837    if( jbot > jtop ){ jj=jbot ; jbot=jtop; jtop=jj; }
00838    if( kbot > ktop ){ kk=kbot ; kbot=ktop; ktop=kk; }
00839    if( ibot < 1 ) ibot = 1 ; if( itop >= nx ) itop = nx-1 ;
00840    if( jbot < 1 ) jbot = 1 ; if( jtop >= ny ) jtop = ny-1 ;
00841    if( kbot < 1 ) kbot = 1 ; if( ktop >= nz ) ktop = nz-1 ;
00842 
00843    /* larger than the largest node id */
00844 
00845    ntop = ag->ixyz[ag->num_ixyz-1].id + 100 ;
00846 
00847    /* scan for voxels that are next to those already mapped,
00848       out to a given level (default max level is 4)         */
00849 
00850    elev = getenv("SUMA_NBHD_LEVEL") ;  /* find level for expanding out */
00851    if( elev != NULL ){
00852       char *cpt ;
00853       ltop = strtol( elev , &cpt , 10 ) ;
00854       if( ltop < 0 || ltop > 7 || (ltop == 0 && *cpt != '\0') ) ltop = 4 ;
00855    } else {
00856       ltop = 4 ;
00857    }
00858 
00859    /* loop over levels */
00860 
00861    for( lev=1 ; lev <= ltop ; lev++ ){  /* if ltop = 0, won't be executed */
00862 
00863     if(PRINT_TRACING){
00864      char str[256]; sprintf(str,"expansion level %d",lev); STATUS(str);
00865     }
00866 
00867     /* loop over voxel 3D indexes */
00868 
00869     for( kk=kbot ; kk < ktop ; kk++ ){
00870      for( jj=jbot ; jj < jtop ; jj++ ){
00871       for( ii=ibot ; ii < itop ; ii++ ){
00872 
00873         ijk = ii+jj*nx+kk*nxy ;         /* index into brick array */
00874         if( vmap[ijk] >= 0 ) continue ; /* already mapped */
00875 
00876         LOAD_IVEC3(iv,ii,jj,kk) ;               /* get Dicom coords */
00877         fv = THD_3dind_to_3dmm( dset , iv ) ;
00878         fv = THD_3dmm_to_dicomm( dset , fv ) ;
00879         UNLOAD_FVEC3(fv,xv,yv,zv) ;             /* coords = (xv,yv,zv) */
00880 
00881         for( pbest=-1,qq=0 ; qq < 26 ; qq++ ){ /* loop over neighbors and */
00882                                                /* find closest mapped pt */
00883                                                /* (at a lower level)    */
00884 
00885           /* pp is the volume map at the qq'th neighboring voxel */
00886 
00887           pp = vmap[(ii+ip[qq][0]) + (jj+ip[qq][1])*nx + (kk+ip[qq][2])*nxy];
00888 
00889           if( pp >= 0 ){                      /* if qq'th nbhr is mapped */
00890              pp = SUMA_VMAP_UNMASK(pp) ;      /* index of mapped node in qq */
00891              xp=xv-ag->ixyz[pp].x ;           /* coords of mapped node */
00892              yp=yv-ag->ixyz[pp].y ;
00893              zp=zv-ag->ixyz[pp].z ;
00894              dd=xp*xp+yp*yp+zp*zp ;           /* dist^2 to mapped pt */
00895 
00896              /* save pbest = index of mapped node closest to (ii,jj,kk)
00897                      dbest = dist^2 of pbest to (ii,jj,kk) voxel center */
00898 
00899              if( pbest >= 0 ){
00900                 if( dd < dbest ){ pbest = pp ; dbest = dd ; }
00901              } else {
00902                 pbest = pp ; dbest = dd ;
00903              }
00904           }
00905         } /* end of loop over 26 neighbors of (ii,jj,kk) voxel */
00906 
00907         /* save closest of the neighbors;
00908            temporarily as a large negative number,
00909            so we won't hit it again in this level of expansion */
00910 
00911         if( pbest >= 0 ) vmap[ijk] = -(pbest+ntop) ; /* closest of the neighbors */
00912 
00913     }}} /* end of loop over voxels (ii,jj,kk) */
00914 
00915     STATUS(".. masking") ;
00916 
00917     /* lmask = 3 bit int for level, shifted to bits 26..28 */
00918 
00919     lmask = SUMA_VMAP_LEVMASK(lev) ;   /* 07 Sep 2001: put on a mask */
00920                                        /* to indicate which level of */
00921                                        /* indirection this voxel was */
00922 
00923     for( kk=kbot ; kk < ktop ; kk++ ){  /* change all the nodes we found */
00924      for( jj=jbot ; jj < jtop ; jj++ ){ /* at this level to non-negative */
00925       for( ii=ibot ; ii < itop ; ii++ ){
00926         ijk = ii+jj*nx+kk*nxy ;
00927         if( vmap[ijk] < -1 ) vmap[ijk] = (-vmap[ijk] - ntop) | lmask ;
00928     }}}
00929 
00930    } /* end of loop over lev */
00931 
00932    RETURN( vmap );
00933 }
00934 
00935 /*-------------------------------------------------------------------------*/
00936 /*! Load the surface for this dataset from its file (if any).
00937 ---------------------------------------------------------------------------*/
00938 
00939 void SUMA_load( THD_3dim_dataset *dset )
00940 {
00941    int ks ;  /* 14 Aug 2002: surface index */
00942 
00943 ENTRY("SUMA_load") ;
00944 
00945    if( !ISVALID_DSET(dset)   ||
00946        dset->su_num   == 0   ||
00947        dset->su_sname == NULL  ) EXRETURN ;
00948 
00949    /* 1st time in: allocate arrays to hold surface data */
00950 
00951    if( dset->su_surf == NULL ){
00952      dset->su_surf   = (SUMA_surface **) calloc(dset->su_num,sizeof(SUMA_surface *));
00953      dset->su_vmap   = (int **)          calloc(dset->su_num,sizeof(int *)         );
00954      dset->su_vnlist = (SUMA_vnlist **)  calloc(dset->su_num,sizeof(SUMA_vnlist *) );
00955    }
00956 
00957    for( ks=0 ; ks < dset->su_num ; ks++ ){       /* loop over surfaces */
00958 
00959      if( dset->su_surf[ks] != NULL ) continue ;  /* already loaded */
00960 
00961      dset->su_surf[ks] = SUMA_read_surface( dset->su_sname[ks] , dset ) ;
00962 
00963      if( dset->su_surf[ks] == NULL ) continue ;
00964 
00965      if( dset->su_vmap[ks] != NULL ) free(dset->su_vmap[ks]) ;
00966 #if 0
00967      dset->su_vmap[ks] = SUMA_map_dset_to_surf( dset->su_surf , dset ) ;
00968 #else
00969      dset->su_vmap[ks] = NULL ;
00970 #endif
00971 
00972      if( dset->su_vnlist[ks] != NULL ){
00973         SUMA_destroy_vnlist( dset->su_vnlist[ks] ) ;
00974         dset->su_vnlist[ks] = NULL ;
00975      }
00976 #if 0
00977      dset->su_vnlist[ks] = SUMA_make_vnlist( dset->su_surf[ks] , dset ) ;
00978 #endif
00979 
00980    }
00981 
00982    EXRETURN ;
00983 }
00984 
00985 /*--------------------------------------------------------------------------*/
00986 /*! Free the surface structs for this dataset (if any).
00987 ----------------------------------------------------------------------------*/
00988 
00989 void SUMA_unload( THD_3dim_dataset *dset )
00990 {
00991    int ks ;  /* 14 Aug 2002: surface index */
00992 
00993 ENTRY("SUMA_unload") ;
00994 
00995    if( !ISVALID_DSET(dset)    ||
00996        dset->su_num   == 0    ||
00997        dset->su_sname == NULL ||
00998        dset->su_surf  == NULL   ) EXRETURN ;
00999 
01000    for( ks=0 ; ks < dset->su_num ; ks++ ){
01001 
01002      if( dset->su_sname[ks] == NULL                   ||
01003          strstr(dset->su_sname[ks],"++LOCK++") != NULL  ) continue ;
01004 
01005      if( dset->su_surf[ks] != NULL ){
01006         SUMA_destroy_surface( dset->su_surf[ks] ) ; dset->su_surf[ks] = NULL ;
01007      }
01008 
01009      if( dset->su_vmap[ks] != NULL ){
01010        free( dset->su_vmap[ks] ) ; dset->su_vmap[ks] = NULL ;
01011      }
01012 
01013      if( dset->su_vnlist[ks] != NULL ){
01014        SUMA_destroy_vnlist( dset->su_vnlist[ks] ) ; dset->su_vnlist[ks] = NULL ;
01015      }
01016 
01017    }
01018 
01019    EXRETURN ;
01020 }
01021 
01022 /*---------------------------------------------------------------------
01023   Read a <SureFit .../> line into a surface
01024     ag   = already existing surface (nodes loaded into this)
01025     lbuf = line containing "<SureFit"
01026     dset = dataset for SureFit-to-DICOM coordinate conversion
01027 -----------------------------------------------------------------------*/
01028 
01029 void SUMA_import_surefit( SUMA_surface *ag, char *lbuf, THD_3dim_dataset *dset )
01030 {
01031    float xx[NBUF],yy[NBUF],zz[NBUF] ;
01032    int   pp[NBUF] ;
01033    int nn , ii , idadd=0 ;
01034    FILE *sfp ;
01035    char sname[1024] , *cpt ;
01036    THD_fvec3 fv ;
01037 
01038 ENTRY("SUMA_import_surefit") ;
01039 
01040    /* scan input line for coord=sname, and extract into sname */
01041 
01042    cpt = strstr(lbuf,"coord=") ;
01043    if( cpt == NULL ){
01044       fprintf(stderr,"** SUMA: Illegal SureFit: no coord=\n** %s\n",lbuf) ;
01045       EXRETURN ;
01046    }
01047    cpt += 6 ;                                  /* skip coord= */
01048    if( *cpt == '\"' || *cpt == '\'' ) cpt++ ;  /* skip quote  */
01049    ii = sscanf(cpt,"%s",sname) ;               /* get sname   */
01050    if( ii == 0 ){
01051       fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
01052       EXRETURN ;
01053    }
01054    ii = strlen(sname) ;
01055    if( ii == 0 ){
01056       fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
01057       EXRETURN ;
01058    }
01059    if( sname[ii-1] == '\'' || sname[ii-1] == '\"' ) sname[ii-1] = '\0' ;
01060    if( strlen(sname) == 0 ){
01061       fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
01062       EXRETURN ;
01063    }
01064 
01065    /* add dataset directory name to start of sname? */
01066 
01067    if( sname[0] != '/' ){
01068       char buf[1024] ;
01069       sprintf(buf,"%s%s",DSET_DIRNAME(dset),sname) ;
01070       strcpy(sname,buf) ;
01071    }
01072 
01073    /* scan line for IDadd=value, and extract into idadd */
01074 
01075    cpt = strstr(lbuf,"IDadd=") ;
01076    if( cpt != NULL ){
01077       cpt += 6 ;
01078       if( *cpt == '\"' || *cpt == '\'' ) cpt++ ;
01079       ii = sscanf(cpt,"%d",&idadd) ;
01080       if( ii == 0 || idadd < 0 ){
01081          fprintf(stderr,"** SUMA: Illegal SureFit: bad IDadd=\n** %s\n",lbuf) ;
01082          EXRETURN ;
01083       }
01084    }
01085 
01086    /* open sname */
01087 
01088    sfp = fopen( sname , "r" ) ;
01089    if( sfp == NULL ){
01090       fprintf(stderr,"** SUMA: Illegal SureFit: can't open file %s\n** %s\n",sname,lbuf) ;
01091       EXRETURN ;
01092    }
01093 
01094    nn = 0 ;
01095 
01096    while(1){
01097       cpt = fgets(sname,1024,sfp) ;  /* read a line */
01098       if( cpt == NULL ) break ;      /* end of file */
01099 
01100       if( strstr(sname,"BeginHeader") != NULL ){  /* skip SureFit header */
01101          do{
01102             cpt = fgets(sname,1024,sfp) ;                /* get next line */
01103             if( cpt == NULL ){ fclose(sfp); EXRETURN; }  /* bad */
01104          } while( strstr(sname,"EndHeader") == NULL ) ;
01105          cpt = fgets(sname,1024,sfp) ;                   /* 1 more line */
01106          if( cpt == NULL ){ fclose(sfp); EXRETURN; }     /* bad */
01107          continue ;                                      /* start over */
01108       }
01109 
01110       ii = sscanf(sname,"%d%f%f%f",pp+nn,xx+nn,yy+nn,zz+nn) ;
01111       if( ii < 4 ) continue ;   /* skip this line; it's bad */
01112 
01113       /* process value to AFNI-ize it from SureFit */
01114 
01115       pp[nn] += idadd ;
01116       LOAD_FVEC3(fv,xx[nn],yy[nn],zz[nn]) ;
01117       fv = THD_surefit_to_dicomm( dset , fv ) ;
01118       UNLOAD_FVEC3(fv,xx[nn],yy[nn],zz[nn]) ;
01119 
01120       nn++ ;
01121       if( nn == NBUF ){
01122          SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
01123          nn = 0 ;
01124       }
01125    } /* end of loop over input lines */
01126 
01127    fclose(sfp) ;
01128 
01129    if( nn > 0 ){
01130       SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
01131    }
01132 
01133    EXRETURN ;
01134 }
01135 
01136 /*-------------------------------------------------------------------------*/
01137 /*! Load the surface name into the dataset struct (derived by replacing
01138     .HEAD with .SURF).
01139 ---------------------------------------------------------------------------*/
01140 
01141 void SUMA_get_surfname( THD_3dim_dataset *dset )
01142 {
01143    char *snam ;
01144    int ii , ks ;
01145 
01146 ENTRY("THD_get_surfname") ;
01147 
01148    if( !ISVALID_DSET(dset) || dset->su_num > 0 ) EXRETURN ;
01149 
01150    snam = strdup( DSET_HEADNAME(dset) ) ;
01151    ii = strlen(snam) ;
01152    if( ii > 5 ){
01153       strcpy(snam+ii-4,"SURF") ;
01154       if( THD_is_file(snam) ){
01155         dset->su_num      = 1 ;
01156         dset->su_sname    = (char **) malloc(sizeof(char *)) ;
01157         dset->su_sname[0] = snam;
01158         EXRETURN;
01159       }
01160    }
01161    free(snam) ; EXRETURN ;  /* .SURF file does not exist */
01162 }
01163 #endif  /* ALLOW_SURFACE_FILES */