00833 {
00834    static char FuncName[]={"SUMA_alloc_problem"};
00835    SUMA_ENTRY;
00836 
00837    printf ("\n\n\aError in memory allocation\n");
00838    printf ("Error origin : %s\n\n",s1);
00839    printf ("Exiting Program ..\n\n");
00840    exit (0);
00841 }

00868 {
00869    int i;
00870    char **A;
00871    static char FuncName[]={"SUMA_allocate2D"};
00872 
00873    SUMA_ENTRY;
00874    
00875    #ifdef USE_SUMA_MALLOC
00876       SUMA_SL_Err("NO LONGER SUPPORTED");
00877       SUMA_RETURN(NULL);
00878    #else
00879       pause_mcw_malloc();
00880    #endif
00881    
00882    /* try to allocate the request */
00883    switch(element_size) {
00884      case sizeof(short): {    /* integer matrix */
00885          short **int_matrix;
00886          int_matrix = (short **)calloc(rows,sizeof(short *));
00887          if(!int_matrix) {
00888              printf("\nError making pointers in %dx%d int matrix\n"
00889                          ,rows,cols);
00890              exit(1);
00891          }
00892          for(i = 0 ; i < rows ; i++) {
00893              int_matrix[i] = (short *)calloc(cols,sizeof(short));
00894              if(!int_matrix[i]) {
00895                  printf("\nError making row %d in %dx%d int matrix\n"
00896                          ,i,rows,cols);
00897                  exit(1);
00898              }
00899          }
00900          A = (char **)int_matrix;
00901          break;
00902      }
00903      case sizeof(float): {    /* float matrix */
00904          float **float_matrix;
00905          float_matrix = (float **)calloc(rows,sizeof(float *));
00906          if(!float_matrix) {
00907              printf("\nError making pointers in %dx%d float matrix\n"
00908                          ,rows,cols);
00909              exit(1);
00910          }
00911          for(i = 0 ; i < rows ; i++) {
00912              float_matrix[i] = (float *)calloc(cols,sizeof(float));
00913              if(!float_matrix[i]) {
00914                  printf("\nError making row %d in %dx%d float matrix\n"
00915                          ,i,rows,cols);
00916                  exit(1);
00917              }
00918          }
00919          A = (char **)float_matrix;
00920          break;
00921      }
00922      case sizeof(double): {   /* double matrix */
00923          double **double_matrix;
00924          double_matrix = (double **)calloc(rows,sizeof(double *));
00925          if(!double_matrix) {
00926              printf("\nError making pointers in %dx%d double matrix\n"
00927                          ,rows,cols);
00928              exit(1);
00929          }
00930          for(i = 0 ; i < rows ; i++) {
00931              double_matrix[i] = (double *)calloc(cols,sizeof(double));
00932              if(!double_matrix[i]) {
00933                  printf("\nError making row %d in %dx%d double matrix\n"
00934                          ,i,rows,cols);
00935                  exit(1);
00936              }
00937          }
00938          A = (char **)double_matrix;
00939          break;
00940      }
00941      default:
00942          printf("\nERROR in matrix_allocate: unsupported type\n");
00943          exit(1);
00944    }
00945    
00946    #ifdef USE_SUMA_MALLOC
00947    SUMA_SL_Err("NO LONGER SUPPORTED");
00948    SUMA_RETURN(NULL);
00949 
00950    #if SUMA_MEMTRACE_FLAG
00951    if (SUMAg_CF->MemTrace) {
00952       ++SUMAg_CF->Mem->N_alloc;
00953       if (SUMAg_CF->Mem->N_MaxPointers <= SUMAg_CF->Mem->N_alloc) {
00954          /* must reallocate */
00955          /* SUMA_ShowMemTrace (SUMAg_CF->Mem, NULL);*/
00956          SUMAg_CF->Mem->N_MaxPointers += SUMA_MEMTRACE_BLOCK;
00957 
00958          SUMAg_CF->Mem->Pointers = (void **)realloc (SUMAg_CF->Mem->Pointers, sizeof(void*) * SUMAg_CF->Mem->N_MaxPointers);
00959          SUMAg_CF->Mem->Size  = (int *)realloc ((void *)SUMAg_CF->Mem->Size, sizeof(int) * SUMAg_CF->Mem->N_MaxPointers);
00960          if (!SUMAg_CF->Mem->Pointers || !SUMAg_CF->Mem->Pointers) {
00961             fprintf (SUMA_STDERR, "Error %s: Failed to reallocate.\nTurning off memory tracing.\n", \
00962                FuncName);
00963             /* free up allocated space, clean up pointers, turn off memory tracing DO NOT USE SUMA_free here*/
00964             if (SUMAg_CF->Mem->Pointers) free(SUMAg_CF->Mem->Pointers); SUMAg_CF->Mem->Pointers = NULL;
00965             if (SUMAg_CF->Mem->Size) free(SUMAg_CF->Mem->Size); SUMAg_CF->Mem->Size = NULL;
00966             SUMAg_CF->MemTrace = 0;
00967             SUMAg_CF->Mem->N_alloc = 0;
00968             SUMAg_CF->Mem->N_MaxPointers =0;
00969          }
00970       }
00971       SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1] = A;
00972       SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1] = rows * cols * element_size;
00973    }
00974    #endif
00975    
00976    #else
00977       resume_mcw_malloc();
00978    #endif
00979    
00980    SUMA_RETURN(A);
00981 }

04829 {
04830    static char FuncName[]={"SUMA_FACESET_FIRST_EDGE_NEIGHB"};
04831    SUMA_FACESET_FIRST_EDGE_NEIGHB *SFFN;
04832    
04833    SUMA_ENTRY;
04834 
04835    SFFN = SUMA_malloc(sizeof(SUMA_FACESET_FIRST_EDGE_NEIGHB));
04836    if (SFFN == NULL) {
04837       fprintf (SUMA_STDERR, "Error %s: Could not allocate for SFFN.\n", FuncName);
04838       SUMA_RETURN (NULL);
04839    }
04840    
04841    SFFN->FirstNeighb = (int **) SUMA_allocate2D(N_FaceSet, SUMA_MAX_FACESET_EDGE_NEIGHB, sizeof(int));
04842    SFFN->N_Neighb = (int *) SUMA_calloc (N_FaceSet, sizeof(int));
04843    if (SFFN->FirstNeighb == NULL || SFFN->N_Neighb == NULL) {
04844       fprintf (SUMA_STDERR, "Error %s: Could not allocate for FirstNeighb or N_Neighb.\n", FuncName);
04845       SUMA_RETURN (NULL);
04846    } 
04847    
04848    SFFN->N_Neighb_max = -1; /* ridiculously low */
04849    SFFN->N_FaceSet = N_FaceSet;
04850    SFFN->N_Neighb_min = 100; /* ridiculously high */
04851    SUMA_RETURN (SFFN);
04852 }

05328 {
05329    static char FuncName[]={"SUMA_Build_FirstNeighb"};
05330    int i, j, n1, n2,  **FirstNeighb, N_ELm1, jj, tmp, TessErr_Cnt=0, IOtrace = 0;
05331    SUMA_Boolean skp;
05332    SUMA_NODE_FIRST_NEIGHB *FN;
05333    
05334    SUMA_ENTRY;
05335 
05336    if (DBG_trace > 1) IOtrace = 1;
05337    
05338    if (el == NULL || N_Node == 0) {
05339       fprintf(SUMA_STDERR, "Error %s: el == NULL or N_Node == 0, nothing to do.\n", FuncName);
05340       SUMA_RETURN (NULL);
05341    }   
05342    
05343    FN = (SUMA_NODE_FIRST_NEIGHB *)SUMA_malloc(sizeof(SUMA_NODE_FIRST_NEIGHB));
05344    if (FN == NULL) {
05345       fprintf(SUMA_STDERR, "Error %s: Could not allocate space for FN\n", FuncName);
05346       SUMA_RETURN (NULL);
05347    }
05348    
05349    FN->N_links = 0;
05350    if (ownerid) sprintf(FN->owner_id, "%s", ownerid);
05351    else FN->owner_id[0] = '\0';
05352    FN->LinkedPtrType = SUMA_LINKED_ND_FRST_NEI_TYPE;
05353    
05354    FN->idcode_str = NULL;
05355    SUMA_NEW_ID(FN->idcode_str, NULL);
05356    
05357    /* allocate space for FN's matrices */
05358    FN->N_Node = N_Node;
05359    FN->N_Neighb_max = 0;
05360    
05361    FN->FirstNeighb = (int **) SUMA_allocate2D(FN->N_Node, SUMA_MAX_NUMBER_NODE_NEIGHB+1, sizeof (int));
05362    FN->N_Neighb = (int *) SUMA_calloc (FN->N_Node, sizeof(int));
05363    FN->NodeId = (int *) SUMA_calloc (FN->N_Node, sizeof(int));
05364    
05365    if (FN->FirstNeighb == NULL || FN->N_Neighb == NULL || FN->NodeId == NULL ){
05366       fprintf(SUMA_STDERR, "Error %s: Could not allocate space forFN->FirstNeighb &/| FN->N_Neighb &/| FN->NodeId.\n", FuncName);
05367       SUMA_RETURN (NULL);
05368    } 
05369    
05370    /*fprintf(SUMA_STDOUT, "%s: Creating list ...\n", FuncName);*/
05371    
05372    FN->N_Neighb_max = 0;
05373    N_ELm1 = el->N_EL-1;
05374    j=0;
05375    while (j < el->N_EL) 
05376    {
05377       n1 = el->EL[j][0];
05378       n2 = el->EL[j][1];
05379       
05380       if (FN->N_Neighb[n1] > SUMA_MAX_NUMBER_NODE_NEIGHB || FN->N_Neighb[n2] > SUMA_MAX_NUMBER_NODE_NEIGHB) {
05381          fprintf(SUMA_STDERR, "Critical Error %s\a: Maximum number of node neighbors for node %d or node %d exceeds %d (SUMA_MAX_NUMBER_NODE_NEIGHB)\n SUMA will try to launch but some functions may not work properly.\n", FuncName, n1, n2, SUMA_MAX_NUMBER_NODE_NEIGHB);
05382       }else {
05383          /*register the neighbors for both nodes*/
05384          FN->NodeId[n1] = n1; /* this field may come in handy when operations need to be performed on subsets of the nodes making up the surface */
05385          FN->NodeId[n2] = n2;
05386          FN->FirstNeighb[n1][FN->N_Neighb[n1]] = n2;
05387          FN->FirstNeighb[n2][FN->N_Neighb[n2]] = n1;
05388 
05389          /* increment neighbor count for nodes in edge */
05390          FN->N_Neighb[n1] += 1;
05391          FN->N_Neighb[n2] += 1;
05392 
05393          if (FN->N_Neighb[n1] > FN->N_Neighb_max) FN->N_Neighb_max = FN->N_Neighb[n1];
05394          if (FN->N_Neighb[n2] > FN->N_Neighb_max) FN->N_Neighb_max = FN->N_Neighb[n2];
05395 
05396          /* skip duplicate edges */
05397          if (j < N_ELm1) {
05398             skp = NOPE;
05399             do {
05400                if (el->EL[j+1][0] == el->EL[j][0] && el->EL[j+1][1] == el->EL[j][1]) {
05401                   ++j;
05402                } else {
05403                   skp = YUP;
05404                }
05405             } while (!skp && j < N_ELm1);
05406          }
05407       }
05408       
05409       ++j;
05410    }/* for j */
05411 
05412    /* now SUMA_reallocate for final FirstNeighb */
05413    FirstNeighb = (int **) SUMA_allocate2D(FN->N_Node, FN->N_Neighb_max, sizeof (int));
05414    if (FirstNeighb == NULL){
05415       fprintf(SUMA_STDERR, "Error %s: Could not allocate space for FirstNeighb\n", FuncName);
05416       SUMA_Free_FirstNeighb (FN);
05417       SUMA_RETURN (NULL);
05418    } 
05419 
05420    /* crop left over allocated space and rearrange neighboring nodes in order */
05421    for (i=0; i < N_Node; ++i) {
05422       #ifdef NoOrder
05423       for (j=0; j < FN->N_Neighb[i]; ++j) {
05424           FirstNeighb[i][j] = FN->FirstNeighb[i][j];
05425       }
05426       #else /* ordered nodes, Tue Jan  7 13:21:57 EST 2003 */
05427         /* copy first node */
05428         FirstNeighb[i][0] = FN->FirstNeighb[i][0];
05429         j = 1;
05430         jj = 1;
05431         while (j < FN->N_Neighb[i]) {
05432             if (SUMA_whichTri (el, i, FirstNeighb[i][jj-1], FN->FirstNeighb[i][j], IOtrace) >= 0) {
05433                FirstNeighb[i][jj] = FN->FirstNeighb[i][j];
05434                /* now swap in FN->FirstNeighb[i] the positions of jj and j */
05435                tmp =  FN->FirstNeighb[i][jj];
05436                FN->FirstNeighb[i][jj] = FN->FirstNeighb[i][j];
05437                FN->FirstNeighb[i][j] = tmp;
05438                ++jj;
05439                j = jj;
05440             } else {
05441                ++j;
05442             }
05443         }
05444         if (jj != FN->N_Neighb[i]) {
05445             if (!TessErr_Cnt) {
05446                fprintf (SUMA_STDERR,"Error %s:\n"
05447                                     " Failed in copying neighbor list! jj=%d, FN->N_Neighb[%d]=%d\n"
05448                                     " If this is a closed surface, the problem is likely due to a\n"
05449                                     " tessellation error. One or more edges may not be part of 2 \n"
05450                                     " and only 2 triangles. Neighbor list for node %d will not be\n"
05451                                     " ordered as connected vertices.\n"
05452                                     " Further occurences of this error will not be reported.\n"
05453                                     ,  FuncName, jj, i, FN->N_Neighb[i], i);
05454             }
05455             ++TessErr_Cnt;
05456             while (jj < FN->N_Neighb[i]) {
05457                FirstNeighb[i][jj] = FN->FirstNeighb[i][jj];
05458                ++jj;
05459             }
05460         }    
05461       #endif
05462    }
05463    if (TessErr_Cnt) {
05464       fprintf (SUMA_STDERR, " %d similar occurences of the error above were found in this mesh.\n", TessErr_Cnt);
05465    }
05466    SUMA_free2D((char **)FN->FirstNeighb, N_Node);
05467    FN->FirstNeighb = FirstNeighb;
05468    /* SUMA_disp_dmat (FN->FirstNeighb, N_Node, FN->N_Neighb_max, 0); */
05469    SUMA_RETURN (FN);
05470 }

00218    {
00219       return (calloc(nmemb, size));
00220    }

00320 {
00321    static char FuncName[]={"SUMA_Cart2Sph"};
00322    float v[3], *f;
00323    int i, i3;
00324    float *sph=NULL;
00325    
00326    SUMA_ENTRY;
00327    
00328    if (Nval <= 0) {
00329       SUMA_RETURN(NULL);
00330    }
00331    
00332    sph = (float *)SUMA_malloc(Nval*sizeof(float)*3);
00333    if (!sph) {
00334       SUMA_SL_Crit("Failed to allocate");
00335       SUMA_RETURN(NULL);
00336    }
00337    
00338    for (i=0; i<Nval; ++i) {
00339       i3 = 3*i;
00340       if (center) {
00341          v[0] = coord[i3+0] - center[0]; 
00342          v[1] = coord[i3+1] - center[1]; 
00343          v[2] = coord[i3+2] - center[2]; 
00344       } else {
00345          v[0] = coord[i3+0]; 
00346          v[1] = coord[i3+1]; 
00347          v[2] = coord[i3+2]; 
00348       }
00349       f = &(sph[i3]);
00350       SUMA_CART_2_SPH(v,f);
00351    }
00352    
00353    SUMA_RETURN(sph);
00354 }

02619 {/*SUMA_compare_int*/
02620     if (*a < *b)
02621       return (-1);
02622    else if (*a == *b)
02623       return (0);
02624    else
02625       return (1);
02626    
02627 }/*SUMA_compare_int*/

06218 {
06219    static char FuncName[]={"SUMA_Convexity"};
06220    float *C=NULL;
06221    
06222    SUMA_ENTRY;
06223    
06224    C = SUMA_Convexity_Engine (NL, N_N, NNL, FN, NULL);
06225    
06226    SUMA_RETURN(C);
06227    
06228 }

06253 {
06254    static char FuncName[]={"SUMA_Convexity_Engine"};
06255    float *C, d, D, dij;
06256    int i, j, jj, in, id, ind, ND;
06257    FILE *fid = NULL;
06258    
06259    SUMA_ENTRY;
06260 
06261    C = NULL;
06262    
06263    /* allocate for C */
06264    C = (float *)SUMA_calloc (N_N, sizeof(float));
06265    
06266    if (C == NULL) {
06267       fprintf (SUMA_STDERR,"Error %s: Could not allocate for C.\n", FuncName);
06268       SUMA_RETURN (C);
06269    }
06270    
06271 
06272    if (DetailFile) {
06273       fprintf (SUMA_STDERR,"%s:\nSaving convexity Info to %s.\n", FuncName, DetailFile);
06274       fid = fopen(DetailFile,"w");
06275    }
06276    
06277    ND = 3;
06278    for (i=0; i < N_N; ++i) {
06279       id = ND * i;
06280       /* the plane at node i, having normal [NNL(id), NNL(id+1), NNL(id+2)] (id = 3*i) has the equation A X + B Y + C Z + D = 0
06281       NNL[id] NL[id]  + NNL[id+1] NNL[id+1]  + NNL[id+2] NL[id+2] + D = 0 */
06282       
06283       D = -NNL[id] * NL[id] - NNL[id+1] * NL[id+1] - NNL[id+2] * NL[id+2];
06284       
06285       if (DetailFile) fprintf(fid,"%d   %d   ", i, FN->N_Neighb[i]);
06286       
06287       for (j=0; j < FN->N_Neighb[i]; ++j) {
06288          /* find the distance between the neighboring node j and the tangent plane at i 
06289          d = (A X + B Y + C Z + D ) / (sqrt(A*A + B*B + C*C))
06290          denominator is norm of Normals which should be 1
06291          */
06292          in = FN->FirstNeighb[i][j];
06293          ind = in * ND;
06294          d = NNL[id] * NL[ind] + NNL[id+1] * NL[ind+1] + NNL[id+2] * NL[ind+2] + D ;
06295          
06296          /* calculate the distance between node i and it's neighbor */
06297          dij = sqrt( (NL[ind] - NL[id]) * (NL[ind] - NL[id]) + (NL[ind+1] - NL[id+1]) * (NL[ind+1] - NL[id+1]) + (NL[ind+2] - NL[id+2]) * (NL[ind+2] - NL[id+2]));
06298          
06299          /* assuming normals are normalized d is the cosine of the angle between the two vectors */
06300          /* if d > 0, then angle is > 0..90 degrees */
06301          
06302          /* as a measure of curvature, compute the sum of signed distances of negihbors to tangent plane at i.
06303          use distances normalized by length of segment ij to account for differences in segment length */
06304           
06305          C[i] -= d/dij; /* used to be C[i] += d/dij; prior to May 06 04 */
06306          
06307          if (DetailFile) fprintf(fid,"%f\t%f\t%f\t", d, dij, d/dij);
06308          
06309       }
06310       
06311       if (DetailFile) {
06312          /* fill with -1 until you reach FN->N_Neighb_max */
06313          for (jj=FN->N_Neighb[i]; jj < FN->N_Neighb_max; ++jj) fprintf(fid,"-1\t-1\t-1\t");
06314          fprintf(fid,"\n");
06315       }
06316    
06317    }
06318    
06319    if (DetailFile) fclose (fid);  /* close previous file */
06320    
06321    #if 0
06322    {
06323       /* Now write the results to disk just for debugging */
06324       fprintf(SUMA_STDOUT,"%s: Writing convexity to Conv.txt ...", FuncName);
06325       fid = fopen("Conv.txt","w");
06326       for (i=0; i < N_N; ++i) {
06327          fprintf(fid,"%f\n", C[i]);
06328       }
06329       fclose (fid);
06330       
06331       fprintf(SUMA_STDOUT,"%s: Done.\n", FuncName);
06332    }
06333    #endif
06334    
06335    SUMA_RETURN (C);
06336 } 

00636 {
00637    SUMA_IRGB *irgb=NULL;
00638    static char FuncName[]={"SUMA_Create_IRGB"};
00639    
00640    SUMA_ENTRY;
00641    
00642    irgb = (SUMA_IRGB *)SUMA_malloc(sizeof(SUMA_IRGB));
00643    
00644    
00645    irgb->i = (int *)SUMA_calloc(n_el, sizeof(int));
00646    irgb->r = (float*)SUMA_calloc(n_el, sizeof(float));
00647    irgb->g = (float*)SUMA_calloc(n_el, sizeof(float));
00648    irgb->b = (float*)SUMA_calloc(n_el, sizeof(float));
00649    irgb->N = n_el;
00650    if (!irgb->i || !irgb->r || !irgb->g || !irgb->b) {
00651       SUMA_S_Crit ("Failed to allocate for i, r, g and/or b.");
00652       if (irgb) SUMA_free(irgb);
00653       SUMA_RETURN (NULL);
00654    }
00655    
00656    SUMA_RETURN(irgb);
00657 }

00227                                                    {
00228    static char FuncName[]={"SUMA_Create_MemTrace"};
00229    SUMA_MEMTRACE_STRUCT *Mem;
00230  
00231    #ifdef USE_SUMA_MALLOC
00232    SUMA_SL_Err("NO LONGER SUPPORTED");
00233    SUMA_RETURN(NULL);
00234    /* you cannot use SUMAg_CF here because the function that allocates for SUMAg_CF calls that one */
00235    
00236    /* DO NOT USE SUMA_malloc function here ! */
00237    Mem = malloc (sizeof(SUMA_MEMTRACE_STRUCT));
00238    /* allocate for the Pointers and Size vectors */
00239    Mem->Pointers = (void **)calloc(SUMA_MEMTRACE_BLOCK, sizeof(void *));
00240    
00241    Mem->Size = (int *)calloc(SUMA_MEMTRACE_BLOCK, sizeof(int));
00242    Mem->N_MaxPointers = SUMA_MEMTRACE_BLOCK;
00243    Mem->N_alloc = 0;
00244    
00245    if (!Mem->Pointers || !Mem->Size) {
00246       fprintf (SUMA_STDERR, "Error %s: Failed to allocate.\n", FuncName);
00247       return (NULL);
00248    }
00249    return(Mem);
00250    #else
00251    return(NULL);
00252    #endif
00253 }

01225 {/*SUMA_disp_dmat*/
01226    char spc [40]; 
01227    int i,j;
01228    static char FuncName[]={"SUMA_disp_dmat"};
01229 
01230    SUMA_ENTRY;
01231 
01232    if (!SpcOpt)
01233       sprintf(spc," ");
01234    else if (SpcOpt == 1)
01235       sprintf(spc,"\t");
01236    else
01237       sprintf(spc," , ");
01238 
01239    fprintf (SUMA_STDOUT,"\n");
01240    for (i=0; i < nr; ++i)
01241       {
01242          for (j=0; j < nc; ++j)
01243                fprintf (SUMA_STDOUT,"%d%s",v[i][j],spc);
01244          fprintf (SUMA_STDOUT,"\n");
01245       }
01246    SUMA_RETURNe;
01247 }/*SUMA_disp_dmat*/

01553 {   int i;
01554    static char FuncName[]={"SUMA_disp_dvect"};
01555    
01556    SUMA_ENTRY;
01557 
01558    fprintf (SUMA_STDOUT,"\n");
01559    if ((l-1) == 0)
01560       fprintf (SUMA_STDOUT, "%d\n",*v);
01561    else 
01562    {
01563    for (i=0;i<l;++i)
01564       fprintf (SUMA_STDOUT,"%d\t",v[i]);
01565 
01566    fprintf (SUMA_STDOUT,"\n");
01567    }
01568    SUMA_RETURNe;
01569 }

01273 {/*SUMA_disp_mat*/
01274    char spc [40]; 
01275     int i,j;
01276    static char FuncName[]={"SUMA_disp_mat"};
01277       
01278    SUMA_ENTRY;
01279 
01280    if (!SpcOpt)
01281       sprintf(spc," ");
01282    else if (SpcOpt == 1)
01283       sprintf(spc,"\t");
01284    else
01285       sprintf(spc," , ");
01286    
01287    fprintf (SUMA_STDOUT,"\n");
01288    for (i=0; i < nr; ++i)
01289       {
01290          for (j=0; j < nc; ++j)
01291                fprintf (SUMA_STDOUT, "%4.2f%s",v[i][j],spc);
01292          fprintf (SUMA_STDOUT,"\n");
01293       }
01294    SUMA_RETURNe;
01295 }/*SUMA_disp_mat*/

01409 {/*SUMA_disp_vecdmat*/
01410    char spc [40]; 
01411    int i,j;
01412    FILE *foutp;
01413    static char FuncName[]={"SUMA_disp_vectdmat"};
01414       
01415    SUMA_ENTRY;
01416 
01417    if (!fout) foutp = stdout;
01418    else foutp = fout;
01419    
01420    if (!SpcOpt)
01421       sprintf(spc," ");
01422    else if (SpcOpt == 1)
01423       sprintf(spc,"\t");
01424    else
01425       sprintf(spc," , ");
01426    
01427    if (!fout) fprintf (SUMA_STDOUT,"\n"); /* a blank 1st line when writing to screen */
01428    switch (d_order) {
01429       case SUMA_ROW_MAJOR:
01430          for (i=0; i < nr; ++i) {
01431             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01432             for (j=0; j < nc; ++j) fprintf (foutp, "%d%s",v[i*nc+j],spc);
01433             fprintf (foutp,"\n");
01434          }
01435          break;
01436       case SUMA_COLUMN_MAJOR:
01437          for (i=0; i < nr; ++i) {
01438             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01439             for (j=0; j < nc; ++j) fprintf (foutp, "%d%s",v[i+j*nr],spc);
01440             fprintf (foutp,"\n");
01441          }
01442          break;
01443       default:
01444          SUMA_SL_Err("Bad order.\n");
01445          SUMA_RETURNe;
01446          break;
01447    }
01448    SUMA_RETURNe;
01449 }/*SUMA_disp_vecdmat*/

01329 {/*SUMA_disp_vecmat*/
01330    char spc [40]; 
01331    int i,j;
01332    FILE *foutp;
01333    static char FuncName[]={"SUMA_disp_vecmat"};
01334    SUMA_Boolean LocalHead = NOPE;
01335       
01336    SUMA_ENTRY;
01337 
01338    if (!fout) foutp = stdout;
01339    else foutp = fout;
01340    
01341    if (LocalHead) fprintf(SUMA_STDERR,"%s:\nExpecting to write %d rows/%d columns\n", FuncName, nr, nc);
01342    
01343    if (!SpcOpt)
01344       sprintf(spc," ");
01345    else if (SpcOpt == 1)
01346       sprintf(spc,"\t");
01347    else
01348       sprintf(spc," , ");
01349    
01350    if (!fout) fprintf (SUMA_STDOUT,"\n"); /* a blank 1st line when writing to screen */
01351    switch (d_order) {
01352       case SUMA_ROW_MAJOR:
01353          for (i=0; i < nr; ++i) {
01354             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01355             for (j=0; j < nc; ++j) fprintf (foutp, "%f%s",v[i*nc+j],spc);
01356             fprintf (foutp,"\n");
01357          }
01358          break;
01359       case SUMA_COLUMN_MAJOR:
01360          for (i=0; i < nr; ++i) {
01361             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01362             for (j=0; j < nc; ++j) fprintf (foutp, "%f%s",v[i+j*nr],spc);
01363             fprintf (foutp,"\n");
01364          }
01365          break;
01366       default:
01367          SUMA_SL_Err("Bad order.\n");
01368          SUMA_RETURNe;
01369          break;
01370    }
01371 
01372    SUMA_RETURNe;
01373 }/*SUMA_disp_vecmat*/

01517 { int i;
01518    static char FuncName[]={"SUMA_disp_vect"};
01519    
01520    SUMA_ENTRY;
01521 
01522    fprintf (SUMA_STDOUT,"\n");
01523    if ((l-1) == 0)
01524       fprintf (SUMA_STDOUT,"%f\n",*v);
01525    else 
01526    {
01527    for (i=0;i<l;++i)
01528                  fprintf (SUMA_STDOUT,"%f\t",v[i]);
01529    fprintf (SUMA_STDOUT,"\n");
01530    }
01531    SUMA_RETURNe;
01532 }

01456 {/*SUMA_disp_vecucmat*/
01457    char spc [40]; 
01458    int i,j;
01459    FILE *foutp;
01460    static char FuncName[]={"SUMA_disp_vecucmat"};
01461       
01462    SUMA_ENTRY;
01463 
01464    if (!fout) foutp = stdout;
01465    else foutp = fout;
01466    
01467    if (!SpcOpt)
01468       sprintf(spc," ");
01469    else if (SpcOpt == 1)
01470       sprintf(spc,"\t");
01471    else
01472       sprintf(spc," , ");
01473    
01474    if (!fout) fprintf (SUMA_STDOUT,"\n"); /* a blank 1st line when writing to screen */
01475    switch (d_order) {
01476       case SUMA_ROW_MAJOR:
01477          for (i=0; i < nr; ++i) {
01478             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01479             for (j=0; j < nc; ++j) fprintf (foutp, "%d%s",v[i*nc+j],spc);
01480             fprintf (foutp,"\n");
01481          }
01482          break;
01483       case SUMA_COLUMN_MAJOR:
01484          for (i=0; i < nr; ++i) {
01485             if (AddRowInd) fprintf (foutp, "%d%s", i, spc);
01486             for (j=0; j < nc; ++j) fprintf (foutp, "%d%s",v[i+j*nr],spc);
01487             fprintf (foutp,"\n");
01488          }
01489          break;
01490       default:
01491          SUMA_SL_Err("Bad order.\n");
01492          SUMA_RETURNe;
01493          break;
01494    }
01495    SUMA_RETURNe;
01496 }/*SUMA_disp_vecucmat*/

02239 {/*SUMA_dPoint_At_Distance*/
02240    static char FuncName[]={"SUMA_dPoint_At_Distance"}; 
02241    double bf, **P2, P1orig[3], Uorig[3];
02242    double m, n, p, q, D, A, B, C, epsi = 0.000001;
02243    int flip, i;
02244    SUMA_Boolean LocalHead = NOPE;
02245    
02246    SUMA_ENTRY;
02247   
02248    SUMA_SL_Warn ("useless piece of junk, use SUMA_POINT_AT_DISTANCE instead!");
02249 
02250    if (d == 0) {
02251       fprintf(SUMA_STDERR,"Error %s: d is 0. Not good, Not good at all.\n", FuncName);
02252       SUMA_RETURN (NULL);
02253    }
02254    
02255    if (LocalHead) {
02256       fprintf (SUMA_STDOUT,"%s: U %f, %f, %f, P1 %f %f %f, d %f\n", FuncName,\
02257          U[0], U[1], U[2], P1[0], P1[1], P1[2], d);
02258    }
02259          
02260    /* store initial values */
02261    P1orig[0] = P1[0];    
02262    P1orig[1] = P1[1]; 
02263    P1orig[2] = P1[2]; 
02264 
02265    Uorig[0] = U[0];
02266    Uorig[1] = U[1];
02267    Uorig[2] = U[2];
02268    
02269    /* normalize U such that U(0) = 1 */
02270    flip = 0;
02271    if (fabs(U[0]) < epsi) { /* must flip X with some other coordinate */
02272       if (fabs(U[1]) > epsi) {/*U[1] != 0; */
02273          U[0] = U[1]; U[1] = 0;
02274          bf = P1[0]; P1[0] = P1[1]; P1[1] = bf;
02275          flip = 1;
02276       } else {   /*U[1] = 0; */
02277          if (fabs(U[2]) > epsi) { /* U[2] != 0 */
02278             U[0] = U[2]; U[2] = 0;
02279             bf = P1[0]; P1[0] = P1[2]; P1[2] = bf;
02280             flip = 2;
02281          } else { /* U[2] = 0 */
02282             fprintf(SUMA_STDERR, "Error %s: 0 direction vector.\n", FuncName);
02283             SUMA_RETURN (NULL);
02284          }
02285       }/*U[1] = 0; */
02286    }/*U[0] = 0; */
02287 
02288    if (LocalHead) fprintf (SUMA_STDERR, "%s: flip = %d\n", FuncName, flip);
02289       
02290    if (LocalHead) fprintf (SUMA_STDERR, "%s: U original: %f, %f, %f\n", FuncName, U[0], U[1], U[2]);
02291    U[1] /= U[0];
02292    U[2] /= U[0];
02293    U[0] = 1.0; 
02294    if (LocalHead) fprintf (SUMA_STDERR, "%s: U normalized: %f, %f, %f\n", FuncName, U[0], U[1], U[2]);
02295    
02296    /* Now U is clean, calculate P2 */   
02297    m = U[1];
02298    n = U[2];
02299 
02300    q = P1[1] - m*P1[0];
02301    p = P1[2] - n*P1[0];
02302 
02303    if (LocalHead) fprintf (SUMA_STDERR, "%s: m=%f n=%f, p=%f, q=%f\n", FuncName, m, n, p, q);
02304 
02305    /* Now find P2 */
02306    A = (1 + n*n + m*m);
02307    B = -2 * P1[0] + 2 * m * (q - P1[1]) + 2 * n * (p - P1[2]);
02308    C = P1[0]*P1[0] + (q - P1[1])*(q - P1[1]) + (p - P1[2])*(p - P1[2]) - d*d;
02309 
02310    D = B*B - 4*A*C;
02311    
02312    if (LocalHead) fprintf (SUMA_STDERR, "%s: A=%f B=%f, C=%f, D=%f\n", FuncName, A, B, C, D);
02313    
02314    if (D < 0) {
02315       fprintf(SUMA_STDERR, "Error %s: Negative Delta: %f.\n"
02316                            "Input values were: \n"
02317                            "U :[%f %f %f]\n"
02318                            "P1:[%f %f %f]\n"
02319                            "d :[%f]\n"
02320                            , FuncName, D, Uorig[0], Uorig[1], Uorig[2], 
02321                            P1orig[0], P1orig[1], P1orig[2], d);
02322       SUMA_RETURN(NULL);
02323    }
02324 
02325    P2 = (double **)SUMA_allocate2D(2,3, sizeof(double));
02326    if (P2 == NULL) {
02327       fprintf(SUMA_STDERR, "Error %s: Could not allocate for 6 floats! What is this? What is the matter with you?!\n", FuncName);
02328       SUMA_RETURN (NULL);
02329    }
02330 
02331    P2[0][0] = (-B + sqrt(D)) / (2 *A);
02332    P2[1][0] = (-B - sqrt(D)) / (2 *A);
02333 
02334    P2[0][1] = m * P2[0][0] + q;
02335    P2[1][1] = m * P2[1][0] + q;
02336 
02337    P2[0][2] = n * P2[0][0] + p;
02338    P2[1][2] = n * P2[1][0] + p;
02339 
02340 
02341    /* if flipping was performed, undo it */
02342    if (flip == 1) {
02343     for (i=0; i < 2; ++i) {
02344        bf = P2[i][1];
02345        P2[i][1] = P2[i][0];
02346        P2[i][0] = bf;
02347       }
02348    } else if (flip == 2){
02349     for (i=0; i < 2; ++i) {
02350        bf = P2[i][2]; 
02351        P2[i][2] = P2[i][0];
02352        P2[i][0] = bf;
02353       }
02354    }   
02355 
02356    for (i=0; i < 3; ++i) {
02357       P1[i] = P1orig[i];
02358       U[i] = Uorig[i];
02359    }
02360 
02361    if (LocalHead) {
02362       fprintf(SUMA_STDOUT,"%s: P1 = %f, %f, %f\n  ", \
02363        FuncName, P1[0], P1[1], P1[2]);
02364       fprintf(SUMA_STDOUT,"%s: P2 = %f, %f, %f\n    %f, %f, %f\n", \
02365        FuncName, P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]);
02366       fprintf(SUMA_STDOUT,"%s: U = %f, %f, %f\n  ", \
02367        FuncName, U[0], U[1], U[2]);
02368    }
02369 
02370    
02371    /* make sure 1st point is along the same direction */
02372    Uorig[0] = P2[0][0] - P1[0]; /* use Uorig, not needed anymore */
02373    Uorig[1] = P2[0][1] - P1[1];
02374    Uorig[2] = P2[0][2] - P1[2];
02375 
02376    SUMA_DOTP_VEC(Uorig, U, bf, 3, double, double)
02377    if (LocalHead) fprintf(SUMA_STDOUT,"%s: Dot product = %f\n", FuncName, bf);
02378    if (bf < 0) {
02379       if (LocalHead) fprintf(SUMA_STDOUT,"%s: Flipping at end...\n", FuncName);
02380       for (i=0; i< 3; ++i) {
02381          bf = P2[0][i];
02382          P2[0][i] = P2[1][i]; P2[1][i] = bf;
02383       }
02384    }
02385    
02386    if (LocalHead) {
02387       fprintf(SUMA_STDOUT,"%s: P2 = %f, %f, %f\n    %f, %f, %f\n", \
02388        FuncName, P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]);
02389    }
02390 SUMA_RETURN (P2);
02391    
02392 }/*SUMA_dPoint_At_Distance*/

02980 {/*SUMA_dqsortrow*/
02981    static char FuncName[]={"SUMA_dqsortrow"}; 
02982    int k,  *I;
02983    SUMA_QSORTROW_INT *Z_Q_dStrct;
02984    
02985    SUMA_ENTRY;
02986    
02987    /* allocate for the structure */
02988    Z_Q_dStrct = (SUMA_QSORTROW_INT *) SUMA_calloc(nr, sizeof (SUMA_QSORTROW_INT));
02989    I = (int *) SUMA_calloc (nr,sizeof(int));
02990    
02991    if (!Z_Q_dStrct || !I)
02992       {
02993       fprintf(SUMA_STDERR,"Error %s: Failed to allocate for Z_Q_dStrct || I\n", FuncName);
02994       SUMA_RETURN (NULL);
02995       }
02996 
02997    for (k=0; k < nr; ++k) /* copy the data into a structure */
02998       {
02999          Z_Q_dStrct[k].x = X[k];
03000          Z_Q_dStrct[k].ncol = nc;
03001          Z_Q_dStrct[k].Index = k;
03002       }
03003 
03004    /* sort the structure by comparing the rows in X */
03005    qsort(Z_Q_dStrct, nr, sizeof(SUMA_QSORTROW_INT), (int(*) (const void *, const void *)) compare_SUMA_QSORTROW_INT);
03006 
03007    /* recover the index table */
03008    for (k=0; k < nr; ++k) 
03009       {
03010          X[k] = Z_Q_dStrct[k].x;
03011          I[k] = Z_Q_dStrct[k].Index;
03012       }
03013    
03014    /* free the structure */
03015    SUMA_free(Z_Q_dStrct);
03016 
03017    /* return */
03018    SUMA_RETURN (I);
03019    
03020    
03021 }/*SUMA_dqsortrow*/

01082  {
01083     static char FuncName[]={"SUMA_error_message"};
01084    
01085    SUMA_ENTRY;
01086 
01087    printf ("\n\n\aError: %s\n",s2);
01088    printf ("Error origin: %s\n\n",s1);
01089    if (ext == 1)
01090       {
01091         printf ("Exiting Program ..\n\n");
01092          exit (0);
01093       }
01094       else SUMA_RETURNe;
01095    
01096   }

01616 {/*SUMA_etime*/
01617    static char FuncName[]={"SUMA_etime"}; 
01618    struct  timeval  tn;
01619    float Time_Fact = 1000000.0;
01620    float delta_t;
01621 
01622    SUMA_ENTRY;
01623 
01624    /* get time */
01625    gettimeofday(&tn,0);
01626    
01627    if (Report)
01628       {
01629          delta_t = (((tn.tv_sec - t->tv_sec)*Time_Fact) + (tn.tv_usec - t->tv_usec)) /Time_Fact;
01630       }
01631    else
01632       {
01633          t->tv_sec = tn.tv_sec;
01634          t->tv_usec = tn.tv_usec;
01635          delta_t = 0.0;
01636       }
01637       
01638    SUMA_RETURN (delta_t);
01639    
01640 }/*SUMA_etime*/

04867 {   
04868    static char FuncName[]={"SUMA_FaceSet_Edge_Neighb"};
04869    int i, i1, F0, F1, in0, in1;
04870    SUMA_FACESET_FIRST_EDGE_NEIGHB *SFFN;
04871    
04872    SUMA_ENTRY;
04873 
04874    
04875    SFFN = SUMA_allocate_FaceSet_Edge_Neighb(N_EL/3);
04876    if (SFFN == NULL) {
04877       fprintf (SUMA_STDERR, "Error %s: Failed in SUMA_allocate_FaceSet_Edge_Neighb.\n", FuncName);
04878       SUMA_RETURN (NULL);
04879    }
04880    
04881    i = 0;
04882    while (i < N_EL-1) {
04883       i1 = i + 1;
04884       if (EL[i][0] != EL[i1][0] || EL[i][1] != EL[i1][1]) {
04885          /* edge is part of one triangle only, a cut edge */
04886          i += 1; /* move to next edge */
04887       } else {
04888          F0 = ELps[i][1]; F1 = ELps[i1][1];
04889          in0 = SFFN->N_Neighb[F0]; in1 = SFFN->N_Neighb[F1];
04890          if (in0 > SUMA_MAX_FACESET_EDGE_NEIGHB -1 || in1 > SUMA_MAX_FACESET_EDGE_NEIGHB -1) {
04891             fprintf (SUMA_STDERR, "Error %s: A faceset has more than three neighbors. Bad surface or non triangular mesh\n", FuncName);
04892             SUMA_RETURN (NULL);
04893          }
04894          SFFN->FirstNeighb[F0][in0] = F1; 
04895          SFFN->FirstNeighb[F1][in1] = F0;
04896          SFFN->N_Neighb[F0] ++;
04897          SFFN->N_Neighb[F1] ++;
04898          if (SFFN->N_Neighb[F0] > SFFN->N_Neighb_max) {
04899             SFFN->N_Neighb_max = SFFN->N_Neighb[F0];
04900          }
04901          if (SFFN->N_Neighb[F1] > SFFN->N_Neighb_max) {
04902             SFFN->N_Neighb_max = SFFN->N_Neighb[F1];
04903          }
04904          if (SFFN->N_Neighb[F0] < SFFN->N_Neighb_min) {
04905             SFFN->N_Neighb_min = SFFN->N_Neighb[F0];
04906          }
04907          if (SFFN->N_Neighb[F1] < SFFN->N_Neighb_min) {
04908             SFFN->N_Neighb_min = SFFN->N_Neighb[F1];
04909          }
04910          
04911          i += 2;
04912       }
04913    
04914    }
04915    
04916    fprintf (SUMA_STDERR, "%s: Done with FaceSet neighbors.\nN_Neighb_max = %d, N_Neighb_min = %d.\n", FuncName, SFFN->N_Neighb_max, SFFN->N_Neighb_min);
04917    #if 0
04918    for (i=0; i< N_FL; ++i) {
04919       fprintf (SUMA_STDERR, "%s: Tri %d, %d neighbs = ", FuncName, i, SFFN->N_Neighb[i]);
04920       for (i1=0; i1<SFFN->N_Neighb[i]; ++i1) {
04921          fprintf (SUMA_STDERR, "%d, ", SFFN->FirstNeighb[i][i1]); 
04922       }
04923       fprintf (SUMA_STDERR, "\n");
04924    }
04925    #endif
04926    
04927    SUMA_RETURN (SFFN);
04928 }   

06815 {
06816    static char FuncName[]={"SUMA_file_suck"};
06817    int  fd , ii = 0;
06818    unsigned long len; 
06819    char * buf = NULL;
06820 
06821    SUMA_ENTRY;
06822    
06823    *nread = 0;
06824    if( fname == NULL || fname[0] == '\0') SUMA_RETURN(0) ;
06825 
06826    len = THD_filesize( fname ) ;
06827    if( len <= 0 ) SUMA_RETURN(buf) ;
06828 
06829    buf = (char *) SUMA_malloc( sizeof(char) * (len+4) ) ;
06830    if( buf == NULL ) SUMA_RETURN(buf) ;
06831 
06832    fd = open( fname , O_RDONLY ) ;
06833    if( fd < 0 ) SUMA_RETURN(buf) ;
06834 
06835    ii = read( fd , buf , len ) ;
06836    close( fd ) ;
06837    if( ii <= 0 ){ SUMA_free(buf) ; buf = NULL; SUMA_RETURN(buf); }
06838    *nread = ii ; SUMA_RETURN(buf) ;
06839 }

06540 {/*SUMA_Find_inIntVect*/
06541    int k, *tmp, *ValLocation;
06542    static char FuncName[]={"SUMA_Find_inIntVect"};
06543    SUMA_Boolean LocalHead = NOPE;
06544    
06545    SUMA_ENTRY;
06546 
06547    /* allocate the maximum  space for ValLocation */
06548    tmp = (int *) SUMA_calloc(xsz,sizeof(int));
06549 
06550    *nValLocation = 0;
06551    for (k = 0 ; k < xsz ; ++k)
06552    {
06553       if (x[k] == val)
06554          {
06555             tmp[*nValLocation] = k;
06556             ++*nValLocation;
06557          }
06558 
06559    }
06560 
06561    if (!*nValLocation)
06562       {
06563          SUMA_free (tmp);
06564          SUMA_RETURN (NULL);
06565       }
06566 
06567    /* Now, allocate just enough space for the SUMA_RETURNing vector */
06568       ValLocation = (int *) SUMA_calloc(*nValLocation,sizeof(int));
06569    /*copy the data into ValLocation*/
06570       SUMA_SCALE_VEC(tmp,ValLocation,1,*nValLocation,int,int);
06571    /* get rid of big array */
06572       SUMA_free(tmp);
06573 
06574    SUMA_RETURN (ValLocation);
06575 
06576 }/*SUMA_Find_inIntVect*/

04672 {
04673    static char FuncName[]={"SUMA_FindEdge"};
04674    int eloc;
04675    SUMA_Boolean LocalHead = NOPE;
04676    
04677    SUMA_ENTRY;
04678 
04679    /* make sure n1 is smallest*/
04680    if (n2 < n1) {
04681       eloc = n2; 
04682       n2 = n1; 
04683       n1 = eloc;
04684    }
04685    
04686    /* first location of edge starting with n1 */
04687    if ((eloc = EL->ELloc[n1]) < 0) {
04688       SUMA_S_Err ("Edge location of n1 not found. WEIRD");
04689       SUMA_RETURN (-1);
04690    }
04691    
04692    /* from there on, look for first occurence of n2 */
04693    do {
04694       /* if (LocalHead) fprintf (SUMA_STDERR,"%s: eloc %d, N_EL %d\n", FuncName, eloc, EL->N_EL);*/
04695       if (EL->EL[eloc][1] == n2) SUMA_RETURN (eloc);
04696       ++eloc;
04697    } while (eloc < EL->N_EL && EL->EL[eloc][0] == n1); 
04698    
04699    /* not found */
04700    SUMA_RETURN (-1);
04701 }

04635 {
04636    static char FuncName[]={"SUMA_FindEdgeInTri"};
04637    int eloc;
04638    
04639    SUMA_ENTRY;
04640 
04641    /* make sure n1 is smallest*/
04642    if (n2 < n1) {
04643       eloc = n2; 
04644       n2 = n1; 
04645       n1 = eloc;
04646    }
04647    
04648    /* first location of edge starting with n1 */
04649    eloc = EL->ELloc[n1];
04650    
04651    /* from there on, look for first occurence of n2 and Tri for hosting triangle*/
04652    do {
04653       if (EL->EL[eloc][1] == n2 && EL->ELps[eloc][1] == Tri) SUMA_RETURN (eloc);
04654       ++eloc;
04655    } while (eloc < EL->N_EL && EL->EL[eloc][0] == n1); 
04656    
04657    /* not found */
04658    SUMA_RETURN (-1);
04659 }

00804 { 
00805    int cnt=0,ex;
00806    float buf;
00807    static char FuncName[]={"SUMA_float_file_size"};
00808    FILE*internal_file;
00809    
00810    SUMA_ENTRY;
00811 
00812    internal_file = fopen (f_name,"r");
00813    if (internal_file == NULL) {
00814                           printf ("\aCould not open %s \n",f_name);
00815                           SUMA_RETURN (-1);
00816                           }
00817    ex = fscanf (internal_file,"%f",&buf);                     
00818    while (ex != EOF)
00819    {
00820      ++cnt;
00821      ex = fscanf (internal_file,"%f",&buf);
00822    }
00823 
00824 
00825    fclose (internal_file);
00826    SUMA_RETURN (cnt);                            
00827 }

02891 {/*SUMA_fqsortrow*/
02892    static char FuncName[]={"SUMA_fqsortrow"}; 
02893    int k, *I;
02894    SUMA_QSORTROW_FLOAT *Z_Q_fStrct;
02895    
02896       
02897    SUMA_ENTRY;
02898 
02899    /* allocate for the structure */
02900    Z_Q_fStrct = (SUMA_QSORTROW_FLOAT *) SUMA_calloc(nr, sizeof (SUMA_QSORTROW_FLOAT));
02901    I = (int *) SUMA_calloc (nr,sizeof(int));
02902    
02903    if (!Z_Q_fStrct || !I)
02904       {
02905       fprintf(SUMA_STDERR,"Error %s: Failed to allocate for Z_Q_fStrct || I\n", FuncName);
02906       SUMA_RETURN (NULL);
02907       }
02908 
02909    for (k=0; k < nr; ++k) /* copy the data into a structure */
02910       {
02911          Z_Q_fStrct[k].x = X[k];
02912          Z_Q_fStrct[k].ncol = nc;
02913          Z_Q_fStrct[k].Index = k;
02914       }
02915 
02916    /* sort the structure by comparing the rows in X */
02917    qsort(Z_Q_fStrct, nr, sizeof(SUMA_QSORTROW_FLOAT), (int(*) (const void *, const void *)) compare_SUMA_QSORTROW_FLOAT);
02918 
02919    /* recover the index table */
02920    for (k=0; k < nr; ++k) 
02921       {
02922          X[k] = Z_Q_fStrct[k].x;
02923          I[k] = Z_Q_fStrct[k].Index;
02924       }
02925    
02926    /* free the structure */
02927    SUMA_free(Z_Q_fStrct);
02928 
02929    /* return */
02930    SUMA_RETURN (I);
02931    
02932    
02933 }/*SUMA_fqsortrow*/

01008 {
01009    int i;
01010    static char FuncName[]={"SUMA_free2D"};
01011    
01012    SUMA_ENTRY;
01013 
01014 
01015       #ifdef USE_SUMA_MALLOC
01016          SUMA_SL_Err("NO LONGER SUPPORTED");
01017          SUMA_RETURNe;
01018 
01019       #if SUMA_MEMTRACE_FLAG
01020          if (SUMAg_CF->MemTrace && a) {
01021             SUMA_Boolean Found = NOPE;
01022             for (i=0; i < SUMAg_CF->Mem->N_alloc && !Found; ++i) {
01023                if (SUMAg_CF->Mem->Pointers[i] == a) {
01024                   SUMAg_CF->Mem->Pointers[i] = SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1];
01025                   SUMAg_CF->Mem->Size[i] = SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1];
01026                   SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1] = NULL;
01027                   SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1] = 0;
01028                   --SUMAg_CF->Mem->N_alloc;
01029                   Found = YUP;
01030                }
01031             }
01032             if (!Found) {
01033               fprintf (SUMA_STDERR, "Error %s: Pointer %p not found in Mem struct. \n", FuncName,a); 
01034             }
01035          }
01036       #endif
01037       #else
01038          pause_mcw_malloc();
01039       #endif
01040       
01041    /* free each row of data */
01042    for(i = 0 ; i < rows ; i++) free(a[i]);
01043 
01044    /* free each row pointer */
01045    free((char *)a);
01046    a = NULL;           /* set to null for error */
01047 
01048    #ifdef USE_SUMA_MALLOC
01049       /* don't use ifndef, keep it parallel with stuff above */
01050       SUMA_SL_Err("NO LONGER SUPPORTED");
01051       SUMA_RETURNe;
01052 
01053    #else
01054       resume_mcw_malloc();
01055    #endif
01056    
01057    SUMA_RETURNe;
01058 }

04309 {
04310    static char FuncName[]={"SUMA_free_Edge_List"};
04311    
04312    SUMA_ENTRY;
04313    if (!SEL) SUMA_RETURNe;
04314    if (SEL->N_links) {
04315       SEL = (SUMA_EDGE_LIST*)SUMA_UnlinkFromPointer((void *)SEL);
04316       SUMA_RETURNe;
04317    }
04318    
04319    if (SEL->EL) SUMA_free2D((char **)SEL->EL, SEL->N_EL);
04320    if (SEL->ELloc) SUMA_free(SEL->ELloc);
04321    if (SEL->ELps) SUMA_free2D((char **)SEL->ELps, SEL->N_EL);
04322    if (SEL->Tri_limb) SUMA_free2D((char **)SEL->Tri_limb, SEL->N_EL/3);
04323    if (SEL->idcode_str) SUMA_free(SEL->idcode_str);
04324    if (SEL) SUMA_free(SEL);
04325    SUMA_RETURNe;
04326 }

05476 {
05477    static char FuncName[]={"SUMA_Free_FirstNeighb"};
05478    SUMA_Boolean LocalHead = NOPE;
05479    
05480    SUMA_ENTRY;
05481    SUMA_LH("Entered");
05482    if (!FN) SUMA_RETURN(YUP);
05483 
05484    if (FN->N_links) {
05485       SUMA_LH("Just a link release");
05486       FN = (SUMA_NODE_FIRST_NEIGHB *)SUMA_UnlinkFromPointer((void *)FN);
05487       SUMA_RETURN (YUP);
05488    }
05489    
05490    /* no more links, go for it */
05491    SUMA_LH("No more links, here we go");
05492    if (FN->idcode_str) SUMA_free(FN->idcode_str); 
05493    if (FN->NodeId) SUMA_free(FN->NodeId);
05494    if (FN->N_Neighb) SUMA_free(FN->N_Neighb);
05495    if (FN->FirstNeighb) SUMA_free2D ((char **)FN->FirstNeighb, FN->N_Node);
05496    if (FN) SUMA_free(FN);
05497    SUMA_RETURN (YUP);
05498 }

00213    {
00214       free(ptr);
00215       return (NULL);
00216    }

00667 {
00668    static char FuncName[]={"SUMA_Free_IRGB"};
00669    
00670    SUMA_ENTRY;
00671    
00672    if (irgb) {
00673       if (irgb->i) SUMA_free(irgb->i);
00674       if (irgb->r) SUMA_free(irgb->r);
00675       if (irgb->g) SUMA_free(irgb->g);
00676       if (irgb->b) SUMA_free(irgb->b);
00677       SUMA_free(irgb);
00678    }
00679    
00680    SUMA_RETURN(NULL);
00681 }

01944 {
01945    static char FuncName[]={"SUMA_Free_IsInBox"};
01946    
01947    SUMA_ENTRY;
01948 
01949    if (IB == NULL) {
01950       fprintf (SUMA_STDERR,"Error SUMA_Free_IsInBox: pointer to null cannot be freed\n");
01951       SUMA_RETURN (NOPE);
01952    }
01953    if (IB->IsIn != NULL) SUMA_free(IB->IsIn);
01954    if (IB->d != NULL) SUMA_free(IB->d);
01955    IB->nIsIn = 0;
01956    SUMA_RETURN (YUP);   
01957 }

01771 {
01772    static char FuncName[]={"SUMA_Free_IsInSphere"};
01773    
01774    SUMA_ENTRY;
01775 
01776    if (IB == NULL) {
01777       fprintf (SUMA_STDERR,"Error SUMA_Free_IsInSphere: pointer to null cannot be freed\n");
01778       SUMA_RETURN (NOPE);
01779    }
01780    if (IB->IsIn != NULL) SUMA_free(IB->IsIn);
01781    if (IB->d != NULL) SUMA_free(IB->d);
01782    IB->nIsIn = 0;
01783    SUMA_RETURN (YUP);   
01784 }

00413                                                              {
00414    static char FuncName[]={"SUMA_Free_MemTrace"};
00415          
00416    #ifdef USE_SUMA_MALLOC
00417    SUMA_SL_Err("NO LONGER SUPPORTED");
00418    return(NOPE);
00419    /* DO NOT USE SUMA_free function here ! */
00420    if (Mem->Pointers) free (Mem->Pointers);
00421    if (Mem->Size) free(Mem->Size);
00422    if (Mem) free (Mem);
00423    #endif
00424    return(YUP);
00425 }

03657 {
03658    static char FuncName[]={"SUMA_Free_MT_intersect_triangle"};
03659    
03660    SUMA_ENTRY;
03661 
03662    if (MTI->t) SUMA_free(MTI->t);
03663    if (MTI->u) SUMA_free(MTI->u);
03664    if (MTI->v) SUMA_free(MTI->v);
03665    if (MTI->isHit) SUMA_free(MTI->isHit);
03666    if (MTI) SUMA_free(MTI);
03667    SUMA_RETURN(NULL);
03668 }

06086 {
06087    static char FuncName[]={"SUMA_Free_SURFACE_CURVATURE"};
06088    
06089    SUMA_ENTRY;
06090 
06091    if (SC == NULL) SUMA_RETURNe;
06092    if (SC->Kp1) SUMA_free(SC->Kp1);
06093    if (SC->Kp2) SUMA_free(SC->Kp2);
06094    if (SC->T1) SUMA_free2D ((char **)SC->T1, SC->N_Node);
06095    if (SC->T2) SUMA_free2D ((char **)SC->T2, SC->N_Node);
06096    if (SC) SUMA_free(SC);
06097    SUMA_RETURNe;
06098 }

03688 {/* SUMA_FromToRotation */
03689    static char FuncName[]={"SUMA_FromToRotation"};
03690    float v[3], vn;
03691    float e, h, f;
03692 
03693    SUMA_ENTRY;
03694 
03695    /*normalize both vectors */
03696    vn = sqrt(v0[0]*v0[0] + v0[1]*v0[1] + v0[2]*v0[2]);
03697    if (vn == 0.0) {
03698       fprintf(SUMA_STDERR,"Error %s: v0 is null.\n",FuncName);
03699       SUMA_RETURN (NOPE);
03700    }
03701    v0[0] /= vn;
03702    v0[1] /= vn;
03703    v0[2] /= vn;
03704 
03705    vn = sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
03706    if (vn == 0.0) {
03707       fprintf(SUMA_STDERR,"Error %s: v1 is null.\n",FuncName);
03708       SUMA_RETURN (NOPE);
03709    }
03710    v1[0] /= vn;
03711    v1[1] /= vn;
03712    v1[2] /= vn;
03713 
03714    SUMA_MT_CROSS(v, v0, v1);
03715    e = SUMA_MT_DOT(v0, v1);
03716    f = (e < 0)? -e:e;
03717    if (f > 1.0 - SUMA_EPSILON)     /* "v0" and "v1"-vector almost parallel */
03718    {
03719       float u[3], v[3]; /* temporary storage vectors */
03720       float x[3];       /* vector most nearly orthogonal v1 "v0" */
03721       float c1, c2, c3; /* coefficients for later use */
03722       int i, j;
03723 
03724       x[0] = (v0[0] > 0.0)? v0[0] : -v0[0];
03725       x[1] = (v0[1] > 0.0)? v0[1] : -v0[1];
03726       x[2] = (v0[2] > 0.0)? v0[2] : -v0[2];
03727 
03728       if (x[0] < x[1]) 
03729       {
03730          if (x[0] < x[2]) 
03731          {
03732            x[0] = 1.0; x[1] = x[2] = 0.0;
03733          }
03734          else 
03735          {
03736            x[2] = 1.0; x[0] = x[1] = 0.0;
03737          }
03738       }
03739       else 
03740       {
03741          if (x[1] < x[2]) 
03742          {
03743            x[1] = 1.0; x[0] = x[2] = 0.0;
03744          }
03745          else 
03746          {
03747            x[2] = 1.0; x[0] = x[1] = 0.0;
03748          }
03749       }
03750 
03751       u[0] = x[0] - v0[0]; u[1] = x[1] - v0[1]; u[2] = x[2] - v0[2];
03752       v[0] = x[0] - v1[0];   v[1] = x[1] - v1[1];   v[2] = x[2] - v1[2];
03753 
03754       c1 = 2.0 / SUMA_MT_DOT(u, u);
03755       c2 = 2.0 / SUMA_MT_DOT(v, v);
03756       c3 = c1 * c2  * SUMA_MT_DOT(u, v);
03757 
03758       for (i = 0; i < 3; i++) {
03759          for (j = 0; j < 3; j++) {
03760            mtx[i][j] =  - c1 * u[i] * u[j]
03761                         - c2 * v[i] * v[j]
03762                         + c3 * v[i] * u[j];
03763          }
03764          mtx[i][i] += 1.0;
03765       }
03766    }
03767    else  /* the most common case, unless "v0"="v1", or "v0"=-"v1" */
03768    {
03769       #if 0
03770          /* unoptimized version - a good compiler will optimize this. */
03771          h = (1.0 - e)/SUMA_MT_DOT(v, v);
03772          mtx[0][0] = e + h * v[0] * v[0];  
03773          mtx[0][1] = h * v[0] * v[1] - v[2]; 
03774          mtx[0][2] = h * v[0] * v[2] + v[1];
03775 
03776          mtx[1][0] = h * v[0] * v[1] + v[2]; 
03777          mtx[1][1] = e + h * v[1] * v[1];    
03778          mtx[1][2] = h * v[1] * v[2] - v[0];
03779 
03780          mtx[2][0] = h * v[0] * v[2] - v[1]; 
03781          mtx[2][1] = h * v[1] * v[2] + v[0]; 
03782          mtx[2][2] = e + h * v[2] * v[2];
03783       #else
03784          /* ...otherwise use this hand optimized version (9 mults less) */
03785          float hvx, hvz, hvxy, hvxz, hvyz;
03786          h = (1.0 - e)/SUMA_MT_DOT(v, v);
03787          hvx = h * v[0];
03788          hvz = h * v[2];
03789          hvxy = hvx * v[1];
03790          hvxz = hvx * v[2];
03791          hvyz = hvz * v[1];
03792          mtx[0][0] = e + hvx * v[0]; 
03793          mtx[0][1] = hvxy - v[2];     
03794          mtx[0][2] = hvxz + v[1];
03795 
03796          mtx[1][0] = hvxy + v[2];  
03797          mtx[1][1] = e + h * v[1] * v[1]; 
03798          mtx[1][2] = hvyz - v[0];
03799 
03800          mtx[2][0] = hvxz - v[1];  
03801          mtx[2][1] = hvyz + v[0];     
03802          mtx[2][2] = e + hvz * v[2];
03803       #endif
03804    }   
03805    
03806    mtx[0][3] = 0.0;
03807    mtx[1][3] = 0.0;
03808    mtx[2][3] = 0.0;
03809    mtx[3][0] = 0.0;
03810    mtx[3][1] = 0.0;
03811    mtx[3][2] = 0.0;
03812    mtx[3][3] = 1.0;
03813    SUMA_RETURN (YUP);
03814 }

04764 {
04765    static char FuncName[] = {"SUMA_Get_Incident"};
04766    int nt, in1, iseek, m_N_EL;
04767    
04768    if (IOtrace) SUMA_ENTRY;
04769 
04770    /*fprintf(SUMA_STDERR,"Entering %s: n1,n2 =%d,%d ...", FuncName,n1,n2);*/
04771    if (n1 > n2) {
04772       /*make the first node be the smallest */
04773       nt = n1;
04774       n1 = n2;
04775       n2 = nt;
04776    }
04777    
04778    /* find the location of the first edge with n1 */
04779    in1 = SEL->ELloc[n1];
04780    iseek = in1;
04781    m_N_EL = SEL->N_EL -1;
04782    *N_Incident = 0;
04783    while (SEL->EL[iseek][0] == n1) {
04784       if (SEL->EL[iseek][1] == n2) {
04785          Incident[*N_Incident] = SEL->ELps[iseek][1]; /* store the faceset index containing the edge */
04786          *N_Incident = *N_Incident + 1;
04787       }
04788       ++iseek;
04789       if (iseek > m_N_EL) {
04790          if (!*N_Incident) fprintf(SUMA_STDERR,"Warning %s: No Incident FaceSets found!\n", FuncName);
04791          if (IOtrace) { SUMA_RETURN (YUP); }
04792          else return(YUP);
04793       }
04794       
04795    }
04796    if (!*N_Incident) fprintf(SUMA_STDERR,"Warning %s: No Incident FaceSets found!\n", FuncName);
04797    /*fprintf(SUMA_STDERR,"Leaving %s.\n", FuncName);*/
04798    if (IOtrace) { SUMA_RETURN(YUP); }
04799    else return(YUP);   
04800 }

04718 {
04719    static char FuncName[] = {"SUMA_Get_NodeIncident"};
04720    int i, n3, N_Neighb;
04721    
04722    SUMA_ENTRY;
04723    
04724    *N_Incident = 0;
04725    
04726    N_Neighb = SO->FN->N_Neighb[n1];
04727    if (N_Neighb < 3) {
04728       fprintf (SUMA_STDERR, "Warning %s: Node %d has less than 3 neighbors.\n", FuncName, n1);
04729       /* nothing found */
04730       SUMA_RETURN(YUP);
04731    }
04732    
04733    i = 0;
04734    while ((i < N_Neighb )) { 
04735       if ( i+1 == N_Neighb) n3 = SO->FN->FirstNeighb[n1][0];
04736       else n3 = SO->FN->FirstNeighb[n1][i+1];
04737       if ((Incident[*N_Incident] = SUMA_whichTri (SO->EL, n1, SO->FN->FirstNeighb[n1][i], n3, 1)) < 0) {
04738          fprintf (SUMA_STDERR, "Error %s: Triangle formed by nodes %d %d %d not found.\n", 
04739             FuncName, n1, SO->FN->FirstNeighb[n1][i], n3);
04740          SUMA_RETURN(NOPE);
04741       }
04742       ++*N_Incident;
04743       ++i;
04744    }
04745 
04746    SUMA_RETURN(YUP);   
04747 }

06117 {
06118    static char FuncName[] = {"SUMA_Householder"};
06119    float mNi, e[3], b[3], mb;
06120    int ii;
06121    #ifdef TAUBIN_Householder
06122    float d[3], s[3], nd, ns;
06123    #endif
06124 
06125    SUMA_ENTRY;
06126    
06127    e[0] = 1.0; e[1] = 0.0; e[2] = 0.0;
06128    
06129    #ifndef TAUBIN_Householder
06130       /* generic algorithm */
06131       mNi = sqrt(Ni[0] * Ni[0] + Ni[1] * Ni[1] + Ni[2] * Ni[2]);
06132       for (ii=0; ii < 3; ++ii) 
06133          b[ii] = Ni[ii] + mNi * e[ii];
06134       mb = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
06135 
06136       if (mb == 0) {
06137          fprintf (SUMA_STDERR,"Error %s: mb = 0\n",FuncName);
06138          SUMA_RETURN (NOPE);
06139       }
06140 
06141       b[0] /= mb; b[1] /= mb; b[2] /= mb;
06142    #else
06143       /* Taubin's algorithm Estimating the tensor of curvature of a surface from a polyhedral approximation */ 
06144       /* calculate difference and sum vectors with their norms (save sqrt for later) */
06145       
06146       d[0] = e[0] - Ni[0]; d[1] = e[1] - Ni[1]; d[2] = e[2] - Ni[2]; 
06147       nd = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
06148       
06149       s[0] = e[0] + Ni[0]; s[1] = e[1] + Ni[1]; s[2] = e[2] + Ni[2];
06150       ns = s[0]*s[0] + s[1]*s[1] + s[2]*s[2];
06151       
06152       if (!nd || !ns) {
06153          fprintf (SUMA_STDERR,"Error %s: nd || ns = 0\n",FuncName);
06154          SUMA_RETURN (NOPE);
06155       }
06156       
06157       if (nd > ns) {
06158          nd = sqrt(nd);
06159          b[0] = d[0] / nd;
06160          b[1] = d[1] / nd;
06161          b[2] = d[2] / nd;
06162          /*Q(:,1) will be equal to -Ni*/
06163          
06164       } else {
06165          ns = sqrt(ns);
06166          b[0] = s[0] / ns;
06167          b[1] = s[1] / ns;
06168          b[2] = s[2] / ns;
06169          /*Q(:,1) will be equal to  Ni */
06170       }
06171       
06172    #endif 
06173    
06174    /* calc Q = I - 2 b * b' */
06175    Q[0][0] = 1 - 2 * b[0] * b[0];
06176    Q[1][0] = - 2 * b[1] * b[0];
06177    Q[2][0] = - 2 * b[2] * b[0];
06178 
06179    Q[0][1] = - 2 * b[0] * b[1];
06180    Q[1][1] = 1 - 2 * b[1] * b[1];
06181    Q[2][1] = - 2 * b[2] * b[1];
06182 
06183    Q[0][2] = - 2 * b[0] * b[2];
06184    Q[1][2] = - 2 * b[1] * b[2];
06185    Q[2][2] = 1 - 2 * b[2] * b[2];
06186 
06187    SUMA_RETURN (YUP);   
06188 }

03970 {
03971    static char FuncName[]={"SUMA_isConsistent"};
03972    static int ic, in, LOC[2], loc[2], d, D;
03973    
03974    SUMA_ENTRY;
03975 
03976    ic = 0;   /* common node index*/
03977    in = 0; /* number of node searched in T */
03978    while (ic < 2 && in < 3) {
03979       if (t[0] == T[in]) {
03980          LOC[ic] = in; /* location of icth node in 1st triangle */
03981          loc[ic] = 0; /* location of icth node in 2nt triangle */
03982          ++ic;
03983       }else {
03984          if (t[1] == T[in]) {
03985             LOC[ic] = in; /* location of icth node in 1st triangle */
03986             loc[ic] = 1; /* location of icth node in 2nt triangle */
03987             ++ic;
03988          }else {
03989             if (t[2] == T[in]) {
03990                LOC[ic] = in; /* location of icth node in 1st triangle */
03991                loc[ic] = 2; /* location of icth node in 2nt triangle */
03992                ++ic;
03993             }
03994          }
03995       }
03996       ++in; /* look for next node */
03997    }
03998    if (ic != 2) {
03999       fprintf(SUMA_STDERR,"Error %s: Triangles do not share 2 nodes.\n", FuncName);
04000       SUMA_RETURN (0);
04001    }
04002    
04003    D = (LOC[1]-LOC[0]);
04004    d = (loc[1]-loc[0]);
04005    /*fprintf(SUMA_STDERR,"%s: T[%d, %d, %d], t[%d, %d, %d]\n", FuncName, T[0], T[1], T[2], t[0], t[1], t[2]);
04006    fprintf(SUMA_STDERR,"%s: LOC[0,1]=[%d, %d], loc[0,1] = [%d, %d]\n", FuncName, LOC[0], LOC[1], loc[0], loc[1]);
04007    fprintf(SUMA_STDERR,"%s: D = %d, d = %d\n", FuncName, D, d);*/   
04008    if (d > 1 || d < -1) d = - d /2 ;
04009    if (D > 1 || D < -1) D = - D /2 ;
04010    /*fprintf(SUMA_STDERR,"%s: D = %d, d = %d\n", FuncName, D, d);*/   
04011 
04012    if (d != D) {
04013       /*fprintf(SUMA_STDERR,"%s: Triangles consistent.\n", FuncName);*/
04014       SUMA_RETURN (1);
04015    }
04016    
04017       
04018    /*fprintf(SUMA_STDERR,"%s: Triangles NOT consistent.\n", FuncName);*/
04019    in = t[0];
04020    t[0] = t[2];
04021    t[2] = in;
04022    SUMA_RETURN (-1);
04023 }

01829 {/*SUMA_isinbox*/
01830    
01831    static char FuncName[]={"SUMA_isinbox"}; 
01832    float t0, t1, t2, hdim0, hdim1, hdim2, *d;
01833    int k , *IsIn, id, ND;
01834    SUMA_ISINBOX IsIn_strct;
01835 
01836    SUMA_ENTRY;
01837    
01838    ND = 3;
01839    /*
01840    fprintf(SUMA_STDOUT,"%f %f %f, %f %f %f, %d, %f, %f, %f\n",\
01841       S_cent[0], S_cent[1], S_cent[2], S_dim[0], S_dim[1], S_dim[2], nr, XYZ[0], XYZ[1], XYZ[2]);
01842    */
01843       
01844    IsIn_strct.nIsIn = 0;   
01845 
01846    hdim0 = S_dim[0]/2;
01847    hdim1 = S_dim[1]/2;
01848    hdim2 = S_dim[2]/2;
01849    
01850    IsIn = (int *) SUMA_calloc (nr, sizeof(int));
01851    d = (float *)SUMA_calloc(nr, sizeof(float));
01852    
01853    if (!IsIn || !d)
01854       {
01855          SUMA_alloc_problem (FuncName);
01856          SUMA_RETURN (IsIn_strct);
01857       }
01858 
01859    if (BoundIn) /* split into two to avoid checking for this condition all the time */
01860       {
01861          /*fprintf(SUMA_STDERR,"%s: inbound\n", FuncName);*/
01862          for (k=0; k < nr; ++k)
01863             {
01864             /*fprintf(SUMA_STDERR,"%s: inbound %d\n", FuncName, k);*/
01865             /* relative distance to center */
01866                id = ND * k;
01867                t0 = hdim0 - fabs(XYZ[id] - S_cent[0]);   
01868                
01869                if (t0 >= 0) {
01870                   t1 = hdim1 - fabs(XYZ[id+1] - S_cent[1]);   
01871                   if (t1 >= 0) {
01872                      t2 = hdim2 - fabs(XYZ[id+2] - S_cent[2]);   
01873                      if (t2 >= 0)
01874                         {
01875                            IsIn[IsIn_strct.nIsIn] = k;
01876                            d[IsIn_strct.nIsIn] = sqrt(t0*t0+t1*t1+t2*t2);
01877                            ++(IsIn_strct.nIsIn);
01878                         }
01879                   }
01880                }
01881             }         
01882             /*fprintf(SUMA_STDERR,"%s: outbound\n", FuncName);*/
01883 
01884       }
01885    else
01886       {
01887          for (k=0; k < nr; ++k)
01888             {
01889                /* relative distance to center */
01890                id = ND * k;
01891                t0 = hdim0 - fabs(XYZ[id] - S_cent[0]);   
01892                
01893                if (t0 > 0) {
01894                   t1 = hdim1 - fabs(XYZ[id+1] - S_cent[1]);   
01895                   if (t1 > 0) {
01896                      t2 = hdim2 - fabs(XYZ[id+2] - S_cent[2]);   
01897                      if (t2 > 0)
01898                         {
01899                            IsIn[IsIn_strct.nIsIn] = k;
01900                            d[IsIn_strct.nIsIn] = sqrt(t0*t0+t1*t1+t2*t2);
01901                            ++(IsIn_strct.nIsIn);
01902                         }
01903                   }
01904                }
01905             }
01906       }
01907    
01908    if (IsIn_strct.nIsIn) {
01909       /*fprintf(SUMA_STDERR,"%s: SUMA_realloc\n", FuncName);*/
01910 
01911       /* get ridd of extra allocation space*/
01912       IsIn_strct.IsIn = (int *) SUMA_calloc (IsIn_strct.nIsIn, sizeof(int));
01913       IsIn_strct.d = (float *)SUMA_calloc(IsIn_strct.nIsIn, sizeof(float));
01914 
01915       if (!IsIn_strct.IsIn || !IsIn_strct.d)
01916          {
01917             IsIn_strct.nIsIn = 0;
01918             SUMA_alloc_problem(FuncName);
01919             SUMA_RETURN (IsIn_strct);
01920          }
01921 
01922       SUMA_COPY_VEC (IsIn, IsIn_strct.IsIn , IsIn_strct.nIsIn, int , int);
01923       SUMA_COPY_VEC (d, IsIn_strct.d, IsIn_strct.nIsIn, float, float);
01924    } else {
01925       /*fprintf(SUMA_STDERR,"%s: NADA\n", FuncName);*/
01926       IsIn_strct.IsIn = NULL;
01927       IsIn_strct.d = NULL;
01928    }
01929    
01930    /*fprintf(SUMA_STDERR,"%s: freeing\n", FuncName);*/
01931    SUMA_free(IsIn);
01932    SUMA_free(d);
01933    /*fprintf(SUMA_STDERR,"%s: freed\n", FuncName);*/
01934 
01935    SUMA_RETURN (IsIn_strct) ;
01936 
01937 }/*SUMA_isinbox*/

01982 {
01983    static char FuncName[]={"SUMA_isinpoly"};
01984    byte *isin=NULL;
01985    int iv, i, ip, counter, ni;
01986    double xinters;
01987    float p1[2], p2[2], p[2], poly[300];
01988    SUMA_Boolean LocalHead = NOPE;
01989    
01990    SUMA_ENTRY;
01991    
01992    *N_in = 0;
01993    if (!usethis) {
01994       isin = (byte *)SUMA_malloc(sizeof(byte)*N_FaceSet);
01995       if (!isin) {
01996          SUMA_SL_Crit("Failed to allocate!");
01997          SUMA_RETURN(NOPE);
01998       }
01999    } else isin = usethis;
02000    if (FaceSetDim > 99) {
02001       SUMA_SL_Err("max FaceSetDim = 99");
02002       SUMA_RETURN(NULL);
02003    }
02004    if (dims[0] < 0 || dims[0] > 2 || dims[1] < 0 || dims[1] > 2) {
02005       SUMA_SL_Err("dims is a 2x1 vector with allowed values of 0 1 or 2 only.");
02006       SUMA_RETURN(NULL);
02007    }
02008 
02009    p[0] = P[dims[0]]; p[1] = P[dims[1]]; /* the point of interest */
02010    for (iv = 0; iv < N_FaceSet; ++iv) {
02011       counter = 0;
02012       for (i=0; i<FaceSetDim; ++i) { /* form the polygon coordinate vector */
02013          ni = FaceSetList[FaceSetDim*iv+i];
02014          poly[3*i] = NodeList[3*ni];  poly[3*i+1] = NodeList[3*ni+1]; poly[3*i+2] = NodeList[3*ni+2];
02015       }
02016       if (culled) if (culled[iv]) continue;
02017       
02018       p1[0] = poly[dims[0]]; p1[1] = poly[dims[1]]; /* the very first point */
02019       for (i=1; i <=FaceSetDim; ++i) {
02020          ip = i % FaceSetDim;
02021          p2[0] = poly[3*ip+dims[0]]; p2[1] = poly[3*ip+dims[1]];
02022          if (p[1] > SUMA_MIN_PAIR(p1[1], p2[1])) {
02023             if (p[1] <= SUMA_MAX_PAIR(p1[1], p2[1])) {
02024                if (p[0] <= SUMA_MAX_PAIR(p1[0], p2[0])) {
02025                   if (p1[1] != p2[1]) {
02026                      xinters = (p[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0];
02027                      if (p1[0] == p2[0] || p[0] <= xinters) {
02028                         counter++; 
02029                      }
02030                   }
02031                }
02032             }
02033          }
02034          p1[0] = p2[0]; p1[1] = p2[1];
02035       }
02036    
02037       if (counter % 2 == 0) { 
02038          isin[iv] = 0;
02039       } else {
02040          isin[iv] = 1; ++(*N_in); /* p is inside polygon iv */
02041          #if 0
02042          if (LocalHead) 
02043          {
02044             int kk;
02045             fprintf(SUMA_STDERR,"\n%%hit!\nPoly = [");
02046             for (kk=0; kk < FaceSetDim; ++kk) {
02047                fprintf(SUMA_STDERR,"%.2f %.2f; ", poly[3*kk+dims[0]] , poly[3*kk+dims[1]]);
02048             } fprintf(SUMA_STDERR,"%.2f %.2f] \np = [%.3f %.3f];", poly[dims[0]], poly[dims[1]], p[0], p[1]);   
02049          }
02050          #endif
02051       }
02052    }
02053 
02054    SUMA_RETURN(isin);
02055 }

01684 {/*SUMA_isinsphere*/
01685    static char FuncName[]={"SUMA_isinsphere"}; 
01686    float *t, t0, t1, t2, ta;
01687    int k, *IsIn, id, ND;
01688    SUMA_ISINSPHERE IsIn_strct;
01689    
01690    SUMA_ENTRY;
01691 
01692    ND = 3;
01693    IsIn_strct.nIsIn = 0;
01694       
01695    t = (float *) SUMA_calloc (nr, sizeof(float));
01696    IsIn = (int *) SUMA_calloc (nr, sizeof(int));
01697    
01698    if (!t || !IsIn)
01699       {
01700          SUMA_alloc_problem (FuncName);
01701          SUMA_RETURN (IsIn_strct);
01702       }
01703    
01704    
01705    if (BoundIn) /* split into two to avoid checking for this condition all the time */
01706       {
01707          for (k=0; k < nr; ++k)
01708             {
01709                id = ND * k;
01710                /* Net distance to center */
01711                t0 = NodeList[id] - S_cent[0];   
01712                t1 = NodeList[id+1] - S_cent[1];   
01713                t2 = NodeList[id+2] - S_cent[2];   
01714 
01715                ta = sqrt (t0 * t0 + t1 * t1 + t2 * t2);
01716                
01717                if (ta <= S_rad)
01718                   {
01719                      IsIn[IsIn_strct.nIsIn] = k;
01720                      t[IsIn_strct.nIsIn] = ta;
01721                      ++(IsIn_strct.nIsIn);
01722                   }
01723             }
01724       }
01725    else
01726       {
01727          for (k=0; k < nr; ++k)
01728             {
01729                id = ND * k;
01730                /* Net distance to center */
01731                t0 = NodeList[id] - S_cent[0];   
01732                t1 = NodeList[id+1] - S_cent[1];   
01733                t2 = NodeList[id+2] - S_cent[2];   
01734 
01735                ta = sqrt (t0 * t0 + t1 * t1 + t2 * t2);
01736                
01737                if (ta < S_rad)
01738                   {
01739                      IsIn[IsIn_strct.nIsIn] = k;
01740                      t[IsIn_strct.nIsIn] = ta;
01741                      ++(IsIn_strct.nIsIn);
01742                   }
01743             }
01744       }
01745          
01746    /* get ridd of extra allocation space*/
01747    IsIn_strct.d = (float *) SUMA_calloc (IsIn_strct.nIsIn, sizeof(float));
01748    IsIn_strct.IsIn = (int *) SUMA_calloc (IsIn_strct.nIsIn, sizeof(int));
01749    
01750    if (!IsIn_strct.d || !IsIn_strct.IsIn )
01751       {
01752          IsIn_strct.nIsIn = 0;
01753          SUMA_alloc_problem(FuncName);
01754          SUMA_RETURN (IsIn_strct);
01755       }
01756    
01757    SUMA_COPY_VEC (t, IsIn_strct.d, IsIn_strct.nIsIn, float , float);
01758    SUMA_COPY_VEC (IsIn, IsIn_strct.IsIn , IsIn_strct.nIsIn, int , int);
01759    
01760    SUMA_free(t);
01761    SUMA_free(IsIn);
01762    
01763    SUMA_RETURN (IsIn_strct);
01764    
01765 }/*SUMA_isinsphere*/

03930 {
03931    static char FuncName[]={"SUMA_isTriLinked"};
03932    int ic, in;
03933    
03934    SUMA_ENTRY;
03935 
03936    ic = 0;   /* common node index*/
03937    in = 0; /* number of node searched in T */
03938    while (ic < 2 && in < 3) {
03939       if (t[0] == T[in]) {
03940          cn[ic] = t[0]; 
03941          ++ic;
03942       }else {
03943          if (t[1] == T[in]) {
03944             cn[ic] = t[1]; 
03945             ++ic;
03946          }else {
03947             if (t[2] == T[in]) {
03948                cn[ic] = t[2]; 
03949                ++ic;
03950             }
03951          }
03952       }
03953       ++in; /* look for next node */
03954    }
03955    
03956    SUMA_RETURN (ic);
03957 }

03167 {
03168    static char FuncName[]={"SUMA_isVoxelIntersect_Triangle"};
03169    int i = 0;
03170    float P0[3], P1[3], iP[3];
03171    
03172    SUMA_ENTRY;
03173    
03174    /* loop accross all 12 edges and find out which pierces the triangle */
03175    for (i=0; i<12; ++i) {
03176       SUMA_EDGE_OF_VOXEL(center, dxyz, i, P0, P1);
03177       if (SUMA_MT_isIntersect_Triangle (P0, P1, vert0, vert1, vert2, iP, NULL, NULL)) {
03178          /* intersects, make sure intersection is between P0 and P1 */
03179          if (SUMA_IS_POINT_IN_SEGMENT(iP, P0, P1)) {
03180             if (0) fprintf(SUMA_STDERR, "%s:\n"
03181                                                 "Triangle %.3f, %.3f, %.3f\n"
03182                                                 "         %.3f, %.3f, %.3f\n"
03183                                                 "         %.3f, %.3f, %.3f\n"
03184                                                 "Intersects voxel at:\n"
03185                                                 "         %.3f, %.3f, %.3f\n", 
03186                                                 FuncName,
03187                                                 vert0[0], vert0[1], vert0[2],
03188                                                 vert1[0], vert1[1], vert1[2],
03189                                                 vert2[0], vert2[1], vert2[2],
03190                                                 center[0], center[1], center[2]);
03191             SUMA_RETURN(YUP);
03192          }
03193       }
03194    }  
03195    SUMA_RETURN(NOPE);
03196 }

01155 {/*SUMA_iswordin*/
01156    int i=0,j=0;
01157    static char FuncName[]={"SUMA_iswordin"};
01158 
01159    SUMA_ENTRY;
01160 
01161    if (sbig == NULL && ssub == NULL) SUMA_RETURN (-2);
01162    if (sbig == NULL || ssub == NULL) SUMA_RETURN (-1);
01163 
01164    if (strlen(sbig) < strlen(ssub))
01165        SUMA_RETURN (0);
01166 
01167    j=0;
01168    while (sbig[i] != '\0' && ssub[j] != '\0')
01169    {
01170        if (sbig[i] == ssub[j])
01171           {
01172              ++j;
01173              /*printf ("j=%d ",j);*/
01174           }
01175        else j=0;
01176    ++i;
01177    }
01178 
01179    if (j == strlen (ssub)) {
01180       SUMA_RETURN (1);
01181    }
01182    else {
01183       SUMA_RETURN (0);
01184    }
01185 
01186 }/*SUMA_iswordin*/

01102 {
01103    static char FuncName[]={"SUMA_iswordin_ci"};
01104    char *sbigc, *ssubc;
01105    int ans;
01106    
01107    SUMA_ENTRY;
01108    sbigc = SUMA_copy_string((char *)sbig);
01109    ssubc = SUMA_copy_string((char *)ssub);
01110    
01111    ans = SUMA_iswordin (sbigc, ssubc);
01112    if (sbigc) SUMA_free(sbigc); sbigc = NULL;
01113    if (ssubc) SUMA_free(ssubc); ssubc = NULL;
01114    
01115    SUMA_RETURN(ans);
01116    
01117 } 

04332 {
04333    static char FuncName[]={"SUMA_Make_Edge_List"};
04334    
04335    SUMA_ENTRY;
04336 
04337    SUMA_RETURN(SUMA_Make_Edge_List_eng(FL, N_FL, N_Node, NodeList, 1, ownerid));
04338 }

04371 {
04372    static char FuncName[]={"SUMA_Make_Edge_List_eng"};
04373    int i, ie, ip, *isort_EL, **ELp, lu, ht, *iTri_limb, icur, in1, in2;
04374    float dx, dy, dz;
04375    SUMA_EDGE_LIST *SEL;
04376    SUMA_Boolean LocalHead = NOPE;
04377    
04378    SUMA_ENTRY;
04379 
04380    if (!FL) {
04381       SUMA_SL_Err("Null FL");
04382       SUMA_RETURN(NULL);
04383    }
04384    if (LocalHead) {
04385       fprintf(SUMA_STDERR,"%s: %d Facesets to work with.\n", FuncName, N_FL);
04386    }
04387    if (!N_FL) {
04388       SUMA_SL_Err("zero N_FL");
04389       SUMA_RETURN(NULL);
04390    }
04391    /* allocate and form the List of edges */
04392    SEL = (SUMA_EDGE_LIST *) SUMA_malloc(sizeof(SUMA_EDGE_LIST));
04393    SEL->idcode_str = NULL;
04394    SUMA_NEW_ID(SEL->idcode_str, NULL); 
04395    SEL->N_links = 0;
04396    if (ownerid) sprintf(SEL->owner_id, "%s", ownerid);
04397    else SEL->owner_id[0] = '\0';
04398    SEL->LinkedPtrType = SUMA_LINKED_OVERLAY_TYPE;
04399 
04400    SEL->N_EL = 3 * N_FL;
04401    SEL->EL = (int **) SUMA_allocate2D (SEL->N_EL, 2, sizeof(int)); /* edge list */
04402    SEL->ELloc = (int *)SUMA_calloc(N_Node, sizeof(int));
04403    SEL->Le = (float *) SUMA_calloc (SEL->N_EL, sizeof(float)); /* length of each edge */
04404    ELp = (int **) SUMA_allocate2D (SEL->N_EL, 2, sizeof(int)); /* edge property list */
04405                                                          /* 1st column, 1 = is flipped from orientation in triangle, -1 as present in triangle 
04406                                                               2nd column, index of triangle (FaceSet) that edge is a part of */    
04407    SEL->ELps = (int **) SUMA_allocate2D (SEL->N_EL, 3, sizeof(int)); /*sorted edge property list */
04408    
04409    /*fprintf(SUMA_STDERR, "%s: SEL->NEL %d\n", FuncName, SEL->N_EL/3);*/
04410    
04411    SEL->Tri_limb = (int **) SUMA_allocate2D (SEL->N_EL/3, 3, sizeof(int)); 
04412    iTri_limb = (int *)SUMA_calloc (SEL->N_EL/3,sizeof(int)); 
04413    
04414    if (SEL == NULL || SEL->EL == NULL || ELp == NULL || SEL->ELps == NULL || SEL->Tri_limb == NULL || iTri_limb== NULL || SEL->ELloc == NULL) {
04415       fprintf(SUMA_STDERR, "Error %s: Failed to allocate for EL, ELp.\n", FuncName);
04416       SUMA_RETURN (NULL);
04417    }
04418 
04419    /* form the edge list */
04420    SUMA_LH("Forming Edge List...\n");
04421    for (i=0; i< N_FL; ++i) {/* begin, form edge list */
04422       /* first edge, 0->1*/
04423       ie = 3*i;
04424       ip = 3*i;
04425       if (FL[ip] > FL[ip+1]) {
04426          /* flip it, to make sorting easier */
04427          SEL->EL[ie][0] = FL[ip+1];
04428          SEL->EL[ie][1] = FL[ip];
04429          /* store parameters */
04430          ELp[ie][0] = 1; /* flip happened */
04431       } else { 
04432          /* no flip necessary */
04433          SEL->EL[ie][0] = FL[ip];
04434          SEL->EL[ie][1] = FL[ip+1];
04435          ELp[ie][0] = -1; /* NO flip happened */
04436       }
04437       ELp[ie][1] = i; /* FaceSetMember */
04438 
04439       /* second edge, 1->2*/
04440       ie += 1;
04441       if (FL[ip+1] > FL[ip+2]) {
04442          /* flip it, to make sorting easier */
04443          SEL->EL[ie][0] = FL[ip+2];
04444          SEL->EL[ie][1] = FL[ip+1];
04445          /* store parameters */
04446          ELp[ie][0] = 1; /* flip happened */
04447       } else { 
04448          /* no flip necessary */
04449          SEL->EL[ie][0] = FL[ip+1];
04450          SEL->EL[ie][1] = FL[ip+2];
04451          ELp[ie][0] = -1; /* NO flip happened */
04452       }
04453       ELp[ie][1] = i; /* FaceSetMember */
04454       
04455       /* third edge, 2->0*/
04456       ie += 1;
04457       if (FL[ip+2] > FL[ip]) {
04458          /* flip it, to make sorting easier */
04459          SEL->EL[ie][0] = FL[ip];
04460          SEL->EL[ie][1] = FL[ip+2];
04461          /* store parameters */
04462          ELp[ie][0] = 1; /* flip happened */
04463       } else { 
04464          /* no flip necessary */
04465          SEL->EL[ie][0] = FL[ip+2];
04466          SEL->EL[ie][1] = FL[ip];
04467          ELp[ie][0] = -1; /* NO flip happened */
04468       }
04469       ELp[ie][1] = i; /* FaceSetMember */
04470       
04471    }/* end, form edge list */
04472    SUMA_LH("Edge list done.");
04473    
04474    if (LocalHead) SUMA_Show_Edge_List(SEL,NULL);
04475    
04476    #if 0
04477       fprintf(SUMA_STDERR,"%s: Node1 Node2 | FlipVal Triangle\n", FuncName); 
04478       for (i=0; i < SEL->N_EL; ++i) {
04479          fprintf (SUMA_STDERR, "%d %d | %d %d\n", SEL->EL[i][0], SEL->EL[i][1], ELp[i][0], ELp[i][1]);
04480       }
04481    #endif
04482 
04483    /* now sort the Edge list */
04484    SUMA_LH("Sorting edge list...");
04485    isort_EL = SUMA_dqsortrow (SEL->EL, SEL->N_EL, 2);
04486    
04487    /* reorder ELp to match sorted EL */
04488    for (i=0; i< SEL->N_EL; ++i) {
04489       SEL->ELps[i][0] = ELp[isort_EL[i]][0];
04490       SEL->ELps[i][1] = ELp[isort_EL[i]][1];
04491    }
04492    
04493    SUMA_LH("Sorting edge list done.");
04494    
04495    if (isort_EL) SUMA_free(isort_EL);
04496    isort_EL = NULL;
04497    
04498    
04499    #if 0
04500       fprintf(SUMA_STDERR,"%s: Node1 Node2 | FlipVal Triangle\n", FuncName); 
04501       for (i=0; i < SEL->N_EL; ++i) {
04502          fprintf (SUMA_STDERR, "%d %d | %d %d\n", SEL->EL[i][0], SEL->EL[i][1], SEL->ELps[i][0], SEL->ELps[i][1]);
04503       }
04504    #endif
04505    
04506    /* calculate the length of each edge */
04507    for (ie=0; ie < SEL->N_EL; ++ie) {
04508       in1 = 3 * SEL->EL[ie][0]; in2 = 3 * SEL->EL[ie][1];
04509       dx = (NodeList[in2] - NodeList[in1]);
04510       dy = (NodeList[in2+1] - NodeList[in1+1]);
04511       dz = (NodeList[in2+2] - NodeList[in1+2]);
04512       SEL->Le[ie] = (float) sqrt (  dx * dx + dy * dy + dz * dz );
04513    }
04514    
04515    /* free unsorted ELp */
04516    if (ELp) SUMA_free2D((char **)ELp, SEL->N_EL);
04517    ELp = NULL;
04518    
04519    SEL->max_N_Hosts = -1;
04520    SEL->min_N_Hosts = 1000;
04521    
04522    /* do a search for some funky stuff */
04523    SUMA_LH("Searching SEL for funky stuff");
04524    SEL->N_Distinct_Edges = 0;
04525    i=0;
04526    while (i < SEL->N_EL) {
04527       /* store the location of this edge for the triangle hosting it */
04528       ht = SEL->ELps[i][1]; /* host triangle index */
04529       SEL->Tri_limb[ht][iTri_limb[ht]] = i;
04530       iTri_limb[ht] += 1;
04531       SEL->N_Distinct_Edges += 1; /* a new edge */ 
04532       SEL->ELps[i][2] = 1; /* number of triangles hosting edge */
04533       lu = 1; 
04534       while (i+lu < SEL->N_EL) {
04535          if (SEL->EL[i+lu][0] == SEL->EL[i][0] && SEL->EL[i+lu][1] == SEL->EL[i][1]) {/* found matching edge */
04536             SEL->ELps[i][2] += 1; /* number of triangles hosting edge */
04537             SEL->ELps[i+lu][2] = -1; /* flag to mean that this edge is a duplicte in the list */
04538 
04539             /* store the location of this edge for the triangle hosting it */
04540             ht = SEL->ELps[i+lu][1]; /* host triangle index */
04541             SEL->Tri_limb[ht][iTri_limb[ht]] = i+lu;
04542             iTri_limb[ht] += 1;
04543 
04544             ++lu;
04545          }else break;
04546          
04547       }
04548       if (SEL->max_N_Hosts < SEL->ELps[i][2]) SEL->max_N_Hosts = SEL->ELps[i][2];
04549       if (SEL->min_N_Hosts > SEL->ELps[i][2]) SEL->min_N_Hosts = SEL->ELps[i][2]; 
04550       i += lu;
04551    }
04552    SUMA_LH("Do adjust your radio.\nNo funk here");
04553    
04554    if (SEL->max_N_Hosts == -1 || SEL->min_N_Hosts == 1000) {
04555       SUMA_SL_Crit("Bad bad news.\nCould not calculate max_N_Hosts &/| min_N_Hosts");
04556       SUMA_free_Edge_List(SEL); SEL = NULL;
04557       SUMA_RETURN(SEL);
04558    }
04559    
04560    {
04561       int winedonce = 0;
04562       if (debug && (SEL->min_N_Hosts == 1 || SEL->max_N_Hosts == 1)) {
04563          winedonce = 1;
04564          fprintf(SUMA_STDERR,"Warning %s:\n Min/Max number of edge hosting triangles: [%d/%d] \n", FuncName, SEL->min_N_Hosts, SEL->max_N_Hosts);
04565          fprintf(SUMA_STDERR," You have edges that form a border in the surface.\n");
04566       }
04567       if (SEL->min_N_Hosts > 2 || SEL->max_N_Hosts > 2) {
04568          winedonce = 1;
04569          fprintf(SUMA_STDERR, "Warning %s:\n"
04570                               "Min/Max number of edge hosting triangles: [%d/%d] \n", FuncName, SEL->min_N_Hosts, SEL->max_N_Hosts);
04571          fprintf(SUMA_STDERR, "Warning %s:\n"
04572                               " You have edges that belong to more than two triangles.\n"
04573                               " Bad for analysis assuming surface is a 2-manifold.\n", FuncName);
04574          if (debug) {
04575             int iii=0;
04576             fprintf(SUMA_STDERR, " These edges are formed by the following nodes:\n");
04577             for (iii = 0; iii < SEL->N_EL; ++iii) { 
04578                if (SEL->ELps[iii][2] > 2) fprintf (SUMA_STDERR," %d: Edge [%d %d] shared by %d triangles.\n", 
04579                                                 iii+1, SEL->EL[iii][0], SEL->EL[iii][1] , SEL->ELps[iii][2] );
04580             }
04581          }
04582       }
04583       if (debug && !winedonce) 
04584          fprintf(SUMA_STDERR,"%s: Min/Max number of edge hosting triangles: [%d/%d] \n", FuncName, SEL->min_N_Hosts, SEL->max_N_Hosts);
04585    }
04586    #if 0
04587       fprintf(SUMA_STDERR,"%s:(ELindex) Node1 Node2 | FlipVal Triangle N_hosts\n", FuncName); 
04588       for (i=0; i < SEL->N_EL; ++i) {
04589          fprintf (SUMA_STDERR, "(%d) %d %d | %d %d %d\n", i, SEL->EL[i][0], SEL->EL[i][1], SEL->ELps[i][0], SEL->ELps[i][1], SEL->ELps[i][2]);
04590       }
04591       fprintf(SUMA_STDERR,"%s:Tri_limb\n", FuncName); 
04592       for (i=0; i < SEL->N_EL/3; ++i) {
04593          fprintf (SUMA_STDERR, "%d %d %d\n", SEL->Tri_limb[i][0], SEL->Tri_limb[i][1], SEL->Tri_limb[i][2]);
04594       }
04595    #endif
04596    
04597    /* store where each node's listing begins
04598    ELloc is used to quickly find a certain edge in EL
04599    to find the edge formed by nodes na-nb
04600    find the minimum of na and nb (say it's nb)
04601    the first reference of an edge containing nb starts at EL(ELloc(nb),:)
04602    NOTE: ELloc contains an entry for each node in FaceSetList, except the largest node index since that's never in the 
04603    first column of EL */
04604    
04605    SUMA_LH("storing locations ...");
04606    for (i=0; i < N_Node; ++i) SEL->ELloc[i] = -1;
04607    i = 0;
04608    icur = SEL->EL[0][0];
04609    SEL->ELloc[icur] = i; 
04610    while (i < SEL->N_EL) {
04611       if (SEL->EL[i][0] != icur) {
04612          icur = SEL->EL[i][0];
04613          SEL->ELloc[icur] = i;
04614       }
04615       ++i;
04616    }
04617    
04618    if (iTri_limb) SUMA_free(iTri_limb); /* Thanks B. Argall */
04619    SUMA_LH("Done with storage, returning...\n");
04620    
04621    SUMA_RETURN (SEL);
04622 }

04950 {
04951    static char FuncName[]={"SUMA_MakeConsistent"};
04952    /* see for more documentation labbook NIH-2 test mesh  p61 */
04953    int i, it, NP, ip, N_flip=0, *isflip, *ischecked, ht0, ht1, NotConsistent, miss, miss_cur, N_iter, EdgeSeed, TriSeed, N_checked;
04954    SUMA_FACESET_FIRST_EDGE_NEIGHB *SFFN;
04955    SUMA_Boolean LocalHead = NOPE;
04956    
04957    SUMA_ENTRY;
04958 
04959    if (detail > 1) LocalHead = YUP;
04960    
04961    if (!SEL || !FL) {
04962       SUMA_SL_Err("NULL input!");
04963       SUMA_RETURN (NOPE);
04964    }
04965    
04966    NP = 3;
04967    isflip = (int *)SUMA_calloc(SEL->N_EL/3, sizeof(int));
04968    ischecked = (int *)SUMA_calloc(SEL->N_EL/3, sizeof(int));
04969    
04970    if (isflip == NULL || ischecked == NULL ) {
04971       fprintf(SUMA_STDERR, "Error %s: Failed to allocate for isflip\n", FuncName);
04972       SUMA_RETURN (NOPE);
04973    }
04974    
04975    
04976    /* Now, go through the sorted edge list and flip what needs flipping*/
04977    N_iter = 0;
04978    miss = 0;
04979    miss_cur = SEL->N_EL;
04980    N_checked = 1;
04981    while (miss_cur != miss) {
04982       miss_cur = miss;
04983       miss = 0;
04984       
04985       /* both methods work just fine here */
04986       #if 0
04987          /*start with the first edge as an edge seed */
04988          EdgeSeed = 0;
04989          i=EdgeSeed;
04990          ht0 = SEL->ELps[i][1];
04991       #else      
04992          /* start with the first edge of the seed triangle as an edge seed */
04993          TriSeed = 0;
04994          i = SEL->Tri_limb[TriSeed][0];
04995          ht0 = TriSeed;
04996       #endif
04997       
04998       ischecked[ht0] = 1;
04999       while (i < SEL->N_EL) {
05000          ht0 = SEL->ELps[i][1];
05001          /* make sure edge is not part of three triangles, if it is, skip it */
05002          if (SEL->ELps[i][2] > 3) {
05003             ++i;
05004             fprintf(SUMA_STDERR, "%s: Bad edge (#%d: %d--%d), part of more than 2 triangles, skip it\n", FuncName, i, SEL->EL[i][0], SEL->EL[i][1]); 
05005             continue;
05006          }
05007          if (SEL->ELps[i][2] == 2) {
05008             /* that's a good edge, see if the next edge after it is consistent */
05009             NotConsistent = SEL->ELps[i][0] * SEL->ELps[i+1][0]; /* if 1 then edges were either both flipped or not flipped in the list */
05010             ht1 = SEL->ELps[i+1][1];
05011             if (ischecked[ht0] && !ischecked[ht1]) {
05012                if (NotConsistent == 0) {
05013                   fprintf(SUMA_STDERR, "Error %s: NotConsistent = 0 here. This should not be.\n", FuncName);
05014                   SUMA_RETURN (NOPE);
05015                }
05016                if (NotConsistent < 0) {
05017                   /* triangles hosting these edges are consistent */
05018                   /* next triangle needs no flipping */
05019                   ischecked[ht1] = 1;
05020                   ++N_checked;
05021                } else {
05022                   /* triangles hosting these edges are NOT consistent */
05023                   /* flip the next triangle */
05024                   ip = NP * ht1;
05025                   it = FL[ip];
05026                   FL[ip] = FL[ip+2];
05027                   FL[ip+2] = it;
05028                   /* Now make sure the flip is reflected in ELps */
05029                   it = SEL->Tri_limb[ht1][0]; SEL->ELps[it][0] *= -1;
05030                   it = SEL->Tri_limb[ht1][1]; SEL->ELps[it][0] *= -1;
05031                   it = SEL->Tri_limb[ht1][2]; SEL->ELps[it][0] *= -1;
05032                   N_flip += 1;
05033                   isflip[ht1] = 1;
05034                   ischecked[ht1] = 1;
05035                   ++N_checked;
05036                }
05037                ++i; 
05038                continue;
05039             } 
05040             
05041             /*try if next edge's host is in good shape */
05042             if (ischecked [ht1] && !ischecked[ht0]) {
05043                if (NotConsistent == 0) {
05044                   fprintf(SUMA_STDERR, "Error %s: NotConsistent = 0 here. This should not be.\n", FuncName);
05045                   SUMA_RETURN (NOPE);
05046                }
05047                if (NotConsistent < 0) {
05048                   /* triangles hosting these edges are consistent */
05049                   /* 1st triangle needs no flipping */
05050                   ischecked[ht0] = 1;
05051                   ++N_checked;
05052                } else {
05053                   /* triangles hosting these edges are NOT consistent */
05054                   /* flip the 1st triangle */
05055                   ip = NP * ht0;
05056                   it = FL[ip];
05057                   FL[ip] = FL[ip+2];
05058                   FL[ip+2] = it;
05059                   /* Now make sure the flip is reflected in ELps */
05060                   it = SEL->Tri_limb[ht0][0]; SEL->ELps[it][0] *= -1;
05061                   it = SEL->Tri_limb[ht0][1]; SEL->ELps[it][0] *= -1;
05062                   it = SEL->Tri_limb[ht0][2]; SEL->ELps[it][0] *= -1;
05063                   N_flip += 1;
05064                   isflip[ht0] = 1;
05065                   ischecked[ht0] = 1;
05066                   ++N_checked;
05067                }
05068                ++i; 
05069                continue;
05070             } 
05071             if (!ischecked[ht0] && !ischecked [ht1]) { /* a good lead that was missed on this pass */
05072                if (LocalHead) fprintf(SUMA_STDERR,"%s: Miss = %d, MissCur = %d\n", FuncName, miss, miss_cur); 
05073                ++miss;
05074             }
05075          }
05076          ++i;   
05077       }
05078       if (LocalHead) fprintf(SUMA_STDERR,"%s: Miss = %d, MissCur = %d\n", FuncName, miss, miss_cur);
05079       ++N_iter;
05080    }
05081 
05082    *trouble = 0;
05083    if (LocalHead) fprintf(SUMA_STDERR,"%s: %d iterations required to check the surface.\n", FuncName, N_iter);
05084    if (detail) fprintf(SUMA_STDERR,"%s: %d/%d (%f%%) triangles checked.\n", FuncName, N_checked, SEL->N_EL/3, (float)N_checked/(SEL->N_EL/3)*100.0);
05085    if (N_checked != SEL->N_EL/3) {
05086       *trouble = 1;
05087    }
05088    if (N_flip) {
05089       *trouble = 1;
05090       if (detail) fprintf(SUMA_STDERR,"%s: %d triangles were flipped to make them consistent with the triangle containing the first edge in the list.\n", FuncName, N_flip);
05091    } else if (detail) fprintf(SUMA_STDERR,"%s: All checked triangles were consistent with the triangle containing the first edge in the list.\n", FuncName);
05092    if (miss) {
05093       if (detail) fprintf(SUMA_STDERR,"%s: %d segments with two neighbors were skipped. Not good in general.\n", FuncName, miss);
05094       *trouble = 1;
05095    }
05096    
05097    #if 0
05098       /* now show the fixed mesh list */
05099       fprintf (SUMA_STDERR,"%s: %d triangles were flipped \n", FuncName, N_flip);
05100       for (i=0; i < SEL->N_EL/3; ++i) {
05101          ip = NP * i;
05102          if (isflip[i]) fprintf (SUMA_STDERR,"\t%d %d %d\t(%d)\t*\n", FL[ip], FL[ip+1], FL[ip+2], ischecked[i]);
05103             else fprintf (SUMA_STDERR,"\t%d %d %d\t(%d)\n", FL[ip], FL[ip+1], FL[ip+2], ischecked[i]);
05104       }
05105    #endif
05106       
05107    /* freedom */
05108    if (LocalHead) fprintf(SUMA_STDERR,"%s: Free time \n", FuncName);
05109    if (isflip) SUMA_free(isflip);
05110    if (ischecked) SUMA_free(ischecked);
05111    if (LocalHead) fprintf(SUMA_STDERR,"%s: returning.\n", FuncName);
05112    
05113    SUMA_RETURN (YUP);
05114 }

00209    {
00210       return (malloc(size));
00211    }

03830 {
03831    double tr, s;
03832    int i,j,k, nxt[3] = {1, 2, 0};
03833    static char FuncName[]={"SUMA_mattoquat"};
03834    
03835    SUMA_ENTRY;
03836 
03837    /* calculate the trace */
03838    tr = mat[0][0] + mat[1][1] + mat[2][2];
03839    if (tr > 0.0) {
03840       s = sqrt(tr + 1.0);
03841       q[3] = s * 0.5;
03842       s = 0.5/s;
03843       
03844       q[0] = (mat[1][2] - mat[2][1])*s;
03845       q[1] = (mat[2][0] - mat[0][2])*s;
03846       q[2] = (mat[0][1] - mat[1][0])*s;
03847    
03848    } /* tr > 0.0 */ else {
03849       i = 0;
03850       if (mat[1][1] > mat[0][0]) i = 1;
03851       if (mat[2][2] > mat[i][i]) i = 2;
03852       j = nxt[i]; k = nxt[j];
03853    
03854       s = sqrt( (mat[i][i] - (mat[j][j]+mat[k][k])) + 1.0);
03855       q[i] = s * 0.5;
03856       s = 0.5/s;
03857       q[3] = (mat[j][k] - mat[k][j])*s;
03858       q[j] = (mat[i][j] + mat[j][i])*s;
03859       q[k] = (mat[i][k] + mat[k][i])*s;
03860    } /* tr < 0.0 */
03861    SUMA_RETURN (YUP);
03862 }

03355 {
03356    static char FuncName[]={"SUMA_MT_intersect_triangle"};
03357    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
03358    double det,inv_det;
03359    int iface, ND, id, NP, ip;
03360    double vert0[3],vert1[3], vert2[3], dir[3], dirn, orig[3];
03361    float tmin, tmax, dii, disttest;
03362    static SUMA_MT_INTERSECT_TRIANGLE *MTI = NULL;
03363    static int N_FaceSet_Previous = 0, entry = 0;
03364    SUMA_Boolean LocalHead = NOPE;
03365    
03366    SUMA_ENTRY;
03367 
03368    tmin = 10000000.0;
03369    tmax = 0.0;
03370    
03371    if (!PrevMTI) { /* nothing preallocated */
03372       entry = 0;
03373       if (LocalHead) fprintf(SUMA_STDERR,"%s: First entry or nothing pre-allocated.\n", FuncName);
03374    } else { /* returning a used MTI, check number of facesets */
03375       if (N_FaceSet_Previous != N_FaceSet) { /* must reallocate */
03376          if (LocalHead) fprintf(SUMA_STDERR,"%s: Reallocating for MTI, a change in number of FaceSets.\n", FuncName);
03377          /* free current MTI */
03378          PrevMTI = SUMA_Free_MT_intersect_triangle (PrevMTI);
03379          entry = 0;
03380       }else if (LocalHead) fprintf(SUMA_STDERR,"%s: Reusing.\n", FuncName);
03381    }
03382    
03383    if (!entry) {
03384       MTI = (SUMA_MT_INTERSECT_TRIANGLE *)SUMA_malloc(sizeof(SUMA_MT_INTERSECT_TRIANGLE));
03385       if (MTI == NULL) {
03386          fprintf(SUMA_STDERR,"Error %s: Failed to allocate for MTI\n", FuncName);
03387          SUMA_RETURN (NULL);
03388       }
03389       MTI->t = NULL;
03390       MTI->u = NULL;
03391       MTI->v = NULL;
03392       MTI->isHit = NULL;
03393    } else {
03394       MTI = PrevMTI;
03395    }
03396 
03397    /* direction from two points */
03398    orig[0] = (double)P0[0];
03399    orig[1] = (double)P0[1];
03400    orig[2] = (double)P0[2];
03401    
03402    dir[0] = (double)P1[0] - orig[0];
03403    dir[1] = (double)P1[1] - orig[1];
03404    dir[2] = (double)P1[2] - orig[2];
03405    dirn = sqrt(dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]);
03406    dir[0] /= dirn;
03407    dir[1] /= dirn;
03408    dir[2] /= dirn;
03409    
03410    if (!entry) {
03411       MTI->isHit = (SUMA_Boolean *)SUMA_malloc(N_FaceSet*sizeof(SUMA_Boolean));
03412       MTI->t = (float *)SUMA_calloc(N_FaceSet, sizeof(float));
03413       MTI->u = (float *)SUMA_calloc(N_FaceSet, sizeof(float));
03414       MTI->v = (float *)SUMA_calloc(N_FaceSet, sizeof(float));
03415    
03416       if (MTI->isHit == NULL || MTI->t == NULL || MTI->u == NULL || MTI->v == NULL) {
03417          fprintf(SUMA_STDERR,"Error : Failed to allocate for MTI->isHit | MTI->t | MTI->u | MTI->v\n");
03418          SUMA_RETURN (NULL);
03419       }
03420    }
03421    
03422    MTI->N_hits = 0; MTI->N_poshits = 0; 
03423    ND = 3;
03424    NP = 3;
03425    for (iface= 0; iface < N_FaceSet; ++iface) {/* iface */
03426       /* set up the coordinates in a humane nomenclature */
03427       ip = NP * iface;
03428       id = ND * FaceSetList[ip];
03429       vert0[0] = (double)NodeList[id];
03430        vert0[1] = (double)NodeList[id+1];
03431       vert0[2] = (double)NodeList[id+2];
03432       
03433       id = ND * FaceSetList[ip+1];
03434       vert1[0] = (double)NodeList[id];
03435        vert1[1] = (double)NodeList[id+1];
03436       vert1[2] = (double)NodeList[id+2];
03437       
03438       id = ND * FaceSetList[ip+2];
03439       vert2[0] = (double)NodeList[id];
03440        vert2[1] = (double)NodeList[id+1];
03441       vert2[2] = (double)NodeList[id+2];
03442       
03443       
03444       /* find vectors for two edges sharing vert0 */
03445       SUMA_MT_SUB(edge1, vert1, vert0);
03446       SUMA_MT_SUB(edge2, vert2, vert0);
03447 
03448       /* begin calculating determinant - also used to calculate U parameter */
03449       SUMA_MT_CROSS(pvec, dir, edge2);
03450 
03451       /* if determinant is near zero, ray lies in plane of triangle */
03452       det = SUMA_MT_DOT(edge1, pvec);
03453 
03454    #ifdef SUMA_MT_TEST_CULL           /* define TEST_CULL if culling is desired */
03455       if (det > -SUMA_EPSILON && det < SUMA_EPSILON)
03456          MTI->isHit[iface] = NOPE;
03457       else {
03458          /* calculate distance from vert0 to ray origin */
03459          SUMA_MT_SUB(tvec, orig, vert0);
03460 
03461          /* calculate U parameter and test bounds */
03462          MTI->u[iface] = (float)SUMA_MT_DOT(tvec, pvec);
03463          if (MTI->u[iface] < 0.0 || MTI->u[iface] > det)
03464             MTI->isHit[iface] = NOPE;
03465          else {
03466             /* prepare to test V parameter */
03467             SUMA_MT_CROSS(qvec, tvec, edge1);
03468 
03469              /* calculate V parameter and test bounds */
03470             MTI->v[iface] = (float)SUMA_MT_DOT(dir, qvec);
03471             if (MTI->v[iface] < 0.0 || MTI->u[iface] + MTI->v[iface] > det)
03472                MTI->isHit[iface] = NOPE;
03473             else {
03474                /* calculate t, scale parameters, ray intersects triangle */
03475                MTI->t[iface] = (float)SUMA_MT_DOT(edge2, qvec);
03476                inv_det = 1.0 / det;
03477                MTI->t[iface] *= (float)inv_det;
03478                MTI->u[iface] *= (float)inv_det;
03479                MTI->v[iface] *= (float)inv_det;         
03480                MTI->isHit[iface] = YUP;
03481                ++MTI->N_hits; 
03482                /* store shortest distance triangle info */
03483                if (MTI->t[iface] < 0) disttest = -MTI->t[iface];
03484                   else { disttest = MTI->t[iface]; ++MTI->N_poshits;}
03485                    
03486                if (disttest < tmin) {
03487                   tmin = disttest;
03488                   MTI->ifacemin = iface;
03489                   /* calculate the location of the intersection in XYZ coords */
03490                   MTI->P[0] = vert0[0] + MTI->u[iface] * (vert1[0] - vert0[0] ) + MTI->v[iface] * (vert2[0] - vert0[0] );
03491                   MTI->P[1] = vert0[1] + MTI->u[iface] * (vert1[1] - vert0[1] ) + MTI->v[iface] * (vert2[1] - vert0[1] );
03492                   MTI->P[2] = vert0[2] + MTI->u[iface] * (vert1[2] - vert0[2] ) + MTI->v[iface] * (vert2[2] - vert0[2] );
03493                   /* find out which node is closest to P */
03494                   MTI->inodeminlocal = 0;
03495                   MTI->d = (vert0[0] - MTI->P[0])*(vert0[0] - MTI->P[0]) + (vert0[1] - MTI->P[1])*(vert0[1] - MTI->P[1]) + (vert0[2] - MTI->P[2])*(vert0[2] - MTI->P[2]);
03496                   dii = (vert1[0] - MTI->P[0])*(vert1[0] - MTI->P[0]) + (vert1[1] - MTI->P[1])*(vert1[1] - MTI->P[1]) + (vert1[2] - MTI->P[2])*(vert1[2] - MTI->P[2]);
03497                   if (dii < MTI->d) {
03498                      MTI->d = dii;
03499                      MTI->inodeminlocal = 1;
03500                   }
03501                   dii = (vert2[0] - MTI->P[0])*(vert2[0] - MTI->P[0]) + (vert2[1] - MTI->P[1])*(vert2[1] - MTI->P[1]) + (vert2[2] - MTI->P[2])*(vert2[2] - MTI->P[2]);
03502                   if (dii < MTI->d) {
03503                      MTI->d = dii;
03504                      MTI->inodeminlocal = 2;
03505                   }
03506                   MTI->d = (float)sqrt((double)MTI->d);
03507                }
03508                if (disttest > tmax) {
03509                   tmax = disttest;
03510                   MTI->ifacemax = iface;
03511                }
03512             }
03513          }
03514       }
03515    #else                    /* the non-culling branch */
03516       if (det > -SUMA_EPSILON && det < SUMA_EPSILON)
03517          MTI->isHit[iface] = NOPE;
03518       else {
03519          inv_det = 1.0 / det;
03520 
03521          /* calculate distance from vert0 to ray origin */
03522          SUMA_MT_SUB(tvec, orig, vert0);
03523 
03524          /* calculate U parameter and test bounds */
03525          MTI->u[iface] = (float)SUMA_MT_DOT(tvec, pvec) * inv_det;
03526          if (MTI->u[iface] < 0.0 || MTI->u[iface] > 1.0)
03527             MTI->isHit[iface] = NOPE;
03528          else {
03529             /* prepare to test V parameter */
03530             SUMA_MT_CROSS(qvec, tvec, edge1);
03531 
03532             /* calculate V parameter and test bounds */
03533             MTI->v[iface] = (float)SUMA_MT_DOT(dir, qvec) * inv_det;
03534             if (MTI->v[iface] < 0.0 || MTI->u[iface] + MTI->v[iface] > 1.0)
03535                MTI->isHit[iface] = NOPE;
03536             else {
03537                /* calculate t, ray intersects triangle */
03538                MTI->t[iface] = (float)SUMA_MT_DOT(edge2, qvec) * inv_det;         
03539                MTI->isHit[iface] = YUP;
03540                ++MTI->N_hits;
03541                /* store shortest distance triangle info */
03542                if (MTI->t[iface] < 0) disttest = -MTI->t[iface];
03543                   else  { disttest = MTI->t[iface]; ++MTI->N_poshits;}
03544                
03545                if (disttest < tmin) {
03546                   tmin = disttest;
03547                   MTI->ifacemin = iface;
03548                   /* calculate the location of the intersection in XYZ coords */
03549                   MTI->P[0] = vert0[0] + MTI->u[iface] * (vert1[0] - vert0[0] ) + MTI->v[iface] * (vert2[0] - vert0[0] );
03550                   MTI->P[1] = vert0[1] + MTI->u[iface] * (vert1[1] - vert0[1] ) + MTI->v[iface] * (vert2[1] - vert0[1] );
03551                   MTI->P[2] = vert0[2] + MTI->u[iface] * (vert1[2] - vert0[2] ) + MTI->v[iface] * (vert2[2] - vert0[2] );
03552                   /* find out which node is closest to P */
03553                   MTI->inodeminlocal = 0;
03554                   MTI->d = (vert0[0] - MTI->P[0])*(vert0[0] - MTI->P[0]) + (vert0[1] - MTI->P[1])*(vert0[1] - MTI->P[1]) + (vert0[2] - MTI->P[2])*(vert0[2] - MTI->P[2]);
03555                   dii = (vert1[0] - MTI->P[0])*(vert1[0] - MTI->P[0]) + (vert1[1] - MTI->P[1])*(vert1[1] - MTI->P[1]) + (vert1[2] - MTI->P[2])*(vert1[2] - MTI->P[2]);
03556                   if (dii < MTI->d) {
03557                      MTI->d = dii;
03558                      MTI->inodeminlocal = 1;
03559                   }
03560                   dii = (vert2[0] - MTI->P[0])*(vert2[0] - MTI->P[0]) + (vert2[1] - MTI->P[1])*(vert2[1] - MTI->P[1]) + (vert2[2] - MTI->P[2])*(vert2[2] - MTI->P[2]);
03561                   if (dii < MTI->d) {
03562                      MTI->d = dii;
03563                      MTI->inodeminlocal = 2;
03564                   }
03565                   MTI->d = (float)sqrt((double)MTI->d);
03566                   ip = NP * iface + MTI->inodeminlocal;
03567                   MTI->inodemin = FaceSetList[ip];
03568                }
03569                if (disttest > tmax) {
03570                   tmax = disttest;
03571                   MTI->ifacemax = iface;
03572                }
03573             }
03574          }
03575       }
03576    #endif
03577    }/*iface */
03578    MTI->N_el = N_FaceSet;
03579    
03580    ++entry;
03581    N_FaceSet_Previous = N_FaceSet;
03582 
03583    SUMA_RETURN (MTI);
03584 }