00001 #include "SUMA_suma.h"
00002 
00003 extern SUMA_CommonFields *SUMAg_CF; 
00004 
00005 #ifdef USE_SUMA_MALLOC
00006    
00007    
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023    void *SUMA_malloc_fn (const char *CF, size_t size)
00024    {
00025       void *ptr;
00026       static char FuncName[]={"SUMA_malloc_fn"};
00027 
00028       #if SUMA_MEMTRACE_FLAG
00029          SUMA_ENTRY;
00030       #endif
00031       
00032       ptr = malloc (size);
00033 
00034       #if SUMA_MEMTRACE_FLAG
00035          if (SUMAg_CF->MemTrace) {
00036             ++SUMAg_CF->Mem->N_alloc;
00037             if (SUMAg_CF->Mem->N_MaxPointers <= SUMAg_CF->Mem->N_alloc) {
00038                
00039                SUMAg_CF->Mem->N_MaxPointers += SUMA_MEMTRACE_BLOCK;
00040                SUMAg_CF->Mem->Pointers = (void **)realloc (SUMAg_CF->Mem->Pointers, sizeof(void*) * SUMAg_CF->Mem->N_MaxPointers);
00041                SUMAg_CF->Mem->Size  = (int *)realloc (SUMAg_CF->Mem->Size, sizeof(int) * SUMAg_CF->Mem->N_MaxPointers);
00042               if (!SUMAg_CF->Mem->Pointers || !SUMAg_CF->Mem->Pointers) {
00043                   fprintf (SUMA_STDERR, "Error %s: Failed to reallocate.\nTurning off memory tracing.\n", \
00044                            FuncName);
00045                   
00046                   if (SUMAg_CF->Mem->Pointers) free(SUMAg_CF->Mem->Pointers); 
00047                   if (SUMAg_CF->Mem->Size) free(SUMAg_CF->Mem->Size); 
00048                   SUMAg_CF->MemTrace = 0;
00049                   SUMAg_CF->Mem->N_alloc = 0;
00050                   SUMAg_CF->Mem->N_MaxPointers = 0;
00051                   SUMA_RETURN(ptr);
00052                }
00053             }
00054             SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1] = ptr;
00055             SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1] = size;
00056          }
00057          SUMA_RETURN(ptr);
00058       #else
00059          return(ptr);
00060       #endif
00061    }
00062 
00063    void *SUMA_realloc_fn (const char *CF, void *ptr, size_t size) 
00064    {
00065       void *ptr2;
00066       int i;
00067       static char FuncName[]={"SUMA_realloc_fn"};
00068 
00069       #if SUMA_MEMTRACE_FLAG
00070          SUMA_ENTRY;
00071       #endif
00072 
00073       
00074       ptr2 = realloc(ptr, size);
00075 
00076       #if SUMA_MEMTRACE_FLAG
00077          if (SUMAg_CF->MemTrace) {
00078             SUMA_Boolean Found = NOPE;
00079             
00080             for (i=0; i < SUMAg_CF->Mem->N_alloc && !Found; ++i) {
00081                if (SUMAg_CF->Mem->Pointers[i] == ptr) {
00082                  
00083                   SUMAg_CF->Mem->Pointers[i] = ptr2;
00084                   SUMAg_CF->Mem->Size[i] = size;
00085                   Found = YUP;
00086                }   
00087             }
00088 
00089             if (!Found) {
00090               fprintf (SUMA_STDERR, "Error %s: Pointer %p not found in Mem struct. \n", FuncName,ptr); 
00091             }
00092          }
00093 
00094          SUMA_RETURN(ptr2);
00095       #else
00096          return(ptr2);
00097       #endif
00098 
00099    }
00100 
00101 
00102 
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00110    void *SUMA_calloc_fn (const char *CF, size_t nmemb, size_t size) {
00111       void *ptr;
00112       static char FuncName[]={"SUMA_calloc_fn"};
00113 
00114       #if SUMA_MEMTRACE_FLAG
00115          SUMA_ENTRY;
00116       #endif
00117 
00118       
00119       ptr = calloc (nmemb, size);
00120 
00121       
00122       #if SUMA_MEMTRACE_FLAG
00123          if (SUMAg_CF->MemTrace) {
00124             ++SUMAg_CF->Mem->N_alloc;
00125             if (SUMAg_CF->Mem->N_MaxPointers <= SUMAg_CF->Mem->N_alloc) {
00126                
00127                
00128                SUMAg_CF->Mem->N_MaxPointers += SUMA_MEMTRACE_BLOCK;
00129 
00130                SUMAg_CF->Mem->Pointers = (void **)realloc (SUMAg_CF->Mem->Pointers, sizeof(void*) * SUMAg_CF->Mem->N_MaxPointers);
00131                SUMAg_CF->Mem->Size  = (int *)realloc ((void *)SUMAg_CF->Mem->Size, sizeof(int) * SUMAg_CF->Mem->N_MaxPointers);
00132                if (!SUMAg_CF->Mem->Pointers || !SUMAg_CF->Mem->Pointers) {
00133                   fprintf (SUMA_STDERR, "Error %s: Failed to reallocate.\nTurning off memory tracing.\n", \
00134                      FuncName);
00135                   
00136                   if (SUMAg_CF->Mem->Pointers) free(SUMAg_CF->Mem->Pointers); SUMAg_CF->Mem->Pointers = NULL;
00137                   if (SUMAg_CF->Mem->Size) free(SUMAg_CF->Mem->Size); SUMAg_CF->Mem->Size = NULL;
00138                   SUMAg_CF->MemTrace = 0;
00139                   SUMAg_CF->Mem->N_alloc = 0;
00140                   SUMAg_CF->Mem->N_MaxPointers =0;
00141                   SUMA_RETURN(ptr);
00142                }
00143             }
00144             SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1] = ptr;
00145             SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1] = nmemb * size;
00146          }
00147          SUMA_RETURN(ptr);
00148       #else
00149          return(ptr);
00150       #endif
00151 
00152    }
00153 
00154 
00155 
00156 
00157 
00158 
00159 
00160 
00161 
00162 
00163 
00164 
00165 
00166 
00167 
00168 
00169    void* SUMA_free_fn(const char *CF, void *ptr)
00170    {
00171       static char FuncName[]={"SUMA_free_fn"};
00172       int i;
00173 
00174       #if SUMA_MEMTRACE_FLAG
00175          SUMA_ENTRY;
00176       #endif
00177 
00178 
00179       
00180       #if SUMA_MEMTRACE_FLAG
00181 
00182          if (SUMAg_CF->MemTrace && ptr) {
00183             SUMA_Boolean Found = NOPE;
00184             for (i=0; i < SUMAg_CF->Mem->N_alloc && !Found; ++i) {
00185                if (SUMAg_CF->Mem->Pointers[i] == ptr) {
00186                   SUMAg_CF->Mem->Pointers[i] = SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1];
00187                   SUMAg_CF->Mem->Size[i] = SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1];
00188                   SUMAg_CF->Mem->Pointers[SUMAg_CF->Mem->N_alloc-1] = NULL;
00189                   SUMAg_CF->Mem->Size[SUMAg_CF->Mem->N_alloc-1] = 0;
00190                   --SUMAg_CF->Mem->N_alloc;
00191                   Found = YUP;
00192                }
00193             }
00194             if (!Found) {
00195               fprintf (SUMA_STDERR, "Error %s: Pointer %p not found in Mem struct. \n\tOffending Calling Function is %s\n", FuncName,ptr, CF); 
00196             }
00197          }
00198 
00199          if (ptr) free (ptr); 
00200          SUMA_RETURN(NULL);
00201       #else
00202          if (ptr) free (ptr); 
00203          return (NULL);
00204       #endif
00205    }
00206 #else
00207    
00208    void *SUMA_malloc_fn (const char *CF, size_t size)
00209    {
00210       return (malloc(size));
00211    }
00212    void* SUMA_free_fn(const char *CF, void *ptr)
00213    {
00214       free(ptr);
00215       return (NULL);
00216    }
00217    void *SUMA_calloc_fn (const char *CF, size_t nmemb, size_t size) 
00218    {
00219       return (calloc(nmemb, size));
00220    }
00221    void *SUMA_realloc_fn (const char *CF, void *ptr, size_t size)
00222    {
00223       return (realloc(ptr, size));
00224    }
00225 #endif
00226 
00227 SUMA_MEMTRACE_STRUCT * SUMA_Create_MemTrace (void) {
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    
00235    
00236    
00237    Mem = malloc (sizeof(SUMA_MEMTRACE_STRUCT));
00238    
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 }
00254  
00255 
00256 
00257 
00258 
00259 
00260 
00261 
00262 
00263 
00264 
00265 
00266 
00267 float * SUMA_Sph2Cart (float *sph, int Nval, float *center ) 
00268 {
00269    static char FuncName[]={"SUMA_Sph2Cart"};
00270    float v[3], *f;
00271    int i, i3;
00272    float *coord=NULL;
00273    
00274    SUMA_ENTRY;
00275    
00276    if (Nval <= 0) {
00277       SUMA_RETURN(NULL);
00278    }
00279    
00280    coord = (float *)SUMA_malloc(Nval*sizeof(float)*3);
00281    if (!coord) {
00282       SUMA_SL_Crit("Failed to allocate");
00283       SUMA_RETURN(NULL);
00284    }
00285    
00286    for (i=0; i<Nval; ++i) {
00287       i3 = 3*i;
00288       f = &(sph[i3]);
00289       SUMA_SPH_2_CART(f, v);
00290       
00291       if (center) {
00292          coord[i3+0] = v[0] + center[0]; 
00293          coord[i3+1] = v[1] + center[1]; 
00294          coord[i3+2] = v[2] + center[2]; 
00295       } else {
00296          coord[i3+0] = v[0]; 
00297          coord[i3+1] = v[1]; 
00298          coord[i3+2] = v[2]; 
00299       }
00300    
00301    }
00302    
00303    SUMA_RETURN(coord);
00304 }
00305 
00306 
00307 
00308 
00309 
00310 
00311 
00312 
00313 
00314 
00315 
00316 
00317 
00318 
00319 float * SUMA_Cart2Sph (float *coord, int Nval, float *center ) 
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 }
00355 
00356 void SUMA_ShowMemTrace (SUMA_MEMTRACE_STRUCT *Mem, FILE *Out) 
00357 {
00358    static char FuncName[]={"SUMA_ShowMemTrace"};
00359    int i, *isort = NULL, *mem_sz_sort = NULL, Tot;
00360    
00361    SUMA_ENTRY;
00362    
00363    #ifdef USE_SUMA_MALLOC
00364    SUMA_SL_Err("NO LONGER SUPPORTED");
00365    SUMA_RETURNe;
00366    if (!Out) Out = SUMA_STDERR;
00367    if (!Mem) {
00368       fprintf (Out,"\nNull struct. Nothing to show.\n");
00369       SUMA_RETURNe;
00370    }
00371    
00372    fprintf (Out,"\nShowing SUMA_MEMTRACE_STRUCT: %p\n", Mem);    
00373    fprintf (Out,"->N_alloc: %d allocated elements.\n", Mem->N_alloc);
00374    fprintf (Out,"->N_MaxPointers: %d\n", Mem->N_MaxPointers);
00375    
00376    
00377    
00378 
00379    mem_sz_sort = (int *)calloc(Mem->N_alloc, sizeof(int));
00380    if (!mem_sz_sort) {
00381       fprintf (SUMA_STDERR, "Error %s: Could not allocate for mem_sz_sort.\n", FuncName);
00382       SUMA_RETURNe;
00383    }
00384    
00385    #if 1
00386    for (i=0; i < Mem->N_alloc; ++i) mem_sz_sort[i] = Mem->Size[i];
00387    isort = SUMA_z_dqsort_nsc (mem_sz_sort, Mem->N_alloc); 
00388    
00389    Tot = 0;
00390    for (i=0; i < Mem->N_alloc; ++i) {
00391       fprintf (Out,"->[%d]\tPointer %p\t %d bytes.\n", i, Mem->Pointers[isort[i]], Mem->Size[isort[i]]);
00392       Tot += Mem->Size[isort[i]];
00393    }
00394    #else
00395      
00396    Tot = 0;
00397    for (i=0; i < Mem->N_alloc; ++i) {
00398       fprintf (Out,"->[%d]\tPointer %p\t %d bytes.\n", i, Mem->Pointers[i], Mem->Size[i]);
00399       Tot += Mem->Size[i];
00400    }
00401    #endif
00402    
00403    fprintf (Out,"Total Memory Allocated %f Mbytes.\n", (float)Tot/1000000.0);
00404    if (mem_sz_sort) free(mem_sz_sort); 
00405    if (isort) free(isort); 
00406    
00407    #endif
00408    
00409    SUMA_RETURNe;
00410    
00411 }
00412 
00413 SUMA_Boolean SUMA_Free_MemTrace (SUMA_MEMTRACE_STRUCT * Mem) {
00414    static char FuncName[]={"SUMA_Free_MemTrace"};
00415          
00416    #ifdef USE_SUMA_MALLOC
00417    SUMA_SL_Err("NO LONGER SUPPORTED");
00418    return(NOPE);
00419    
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 }
00426 
00427 
00428 
00429 
00430 
00431 
00432 
00433 
00434 
00435 
00436 
00437 
00438 
00439 
00440 
00441 
00442 
00443 
00444 
00445 
00446 
00447 
00448 
00449 
00450 
00451 
00452 
00453    
00454 int SUMA_Read_dfile (int *x,char *f_name,int n_points)
00455    
00456 { 
00457    int cnt=0,ex,dec;
00458    static char FuncName[]={"SUMA_Read_dfile"};
00459    FILE*internal_file;
00460 
00461    SUMA_ENTRY;
00462 
00463    internal_file = fopen (f_name,"r");
00464    if (internal_file == NULL) {
00465                           fprintf(SUMA_STDERR, "\aCould not open %s \n",f_name);
00466                           fprintf(SUMA_STDERR, "Exiting @ SUMA_Read_file function\n");
00467                           exit (0);
00468                           }
00469    ex = fscanf (internal_file,"%d",&x[cnt]);                     
00470    while (ex != EOF)
00471    {
00472      ++cnt;
00473      
00474 
00475 
00476 
00477 
00478 
00479 
00480 
00481 
00482      ex = fscanf (internal_file,"%d",&x[cnt]);
00483 
00484      if ((n_points != 0) && (cnt == n_points)) ex = EOF;
00485    }
00486 
00487    if (cnt < n_points) 
00488       {
00489        fprintf(SUMA_STDERR, "\a\nAttempt to read %d points failed,\n",n_points);
00490        fprintf(SUMA_STDERR, " file contains %d points only.\n",cnt);
00491        do {
00492 
00493        fprintf(SUMA_STDERR, "End Execution (Yes (1) No (0) ? : ");
00494        ex=scanf ("%d",&dec);
00495        } while (ex != 1 || (dec != 1 && dec !=0));
00496        if (dec)
00497         {
00498           fprintf(SUMA_STDERR, "Exiting @ SUMA_Read_file function\n");
00499       exit (0);
00500       }
00501       else fprintf(SUMA_STDERR, "\nContinuing execution with %d points\n",cnt);
00502 
00503       }
00504 
00505    fclose (internal_file);
00506    SUMA_RETURN (cnt);                            
00507 }
00508 int SUMA_Read_file (float *x,char *f_name,int n_points)
00509    
00510 { 
00511    int cnt=0,ex,dec;
00512    FILE*internal_file;
00513    static char FuncName[]={"SUMA_Read_file"};
00514    
00515    SUMA_ENTRY;
00516 
00517    internal_file = fopen (f_name,"r");
00518    if (internal_file == NULL) {
00519                           fprintf(SUMA_STDERR, "\aCould not open %s \n",f_name);
00520                           fprintf(SUMA_STDERR, "Exiting @ SUMA_Read_file function\n");
00521                           exit (0);
00522                           }
00523    ex = fscanf (internal_file,"%f",&x[cnt]);                     
00524    while (ex != EOF)
00525    {
00526      ++cnt;
00527      
00528      ex = fscanf (internal_file,"%f",&x[cnt]);
00529 
00530      if ((n_points != 0) && (cnt == n_points)) ex = EOF;
00531    }
00532 
00533    if (cnt < n_points) 
00534       {
00535        fprintf(SUMA_STDERR, "\a\nAttempt to read %d points failed,\n",n_points);
00536        fprintf(SUMA_STDERR, " file contains %d points only.\n",cnt);
00537        do {
00538 
00539        fprintf(SUMA_STDERR, "End Execution (Yes (1) No (0) ? : ");
00540        ex=scanf ("%d",&dec);
00541        } while (ex != 1 || (dec != 1 && dec !=0));
00542        if (dec)
00543         {
00544           fprintf(SUMA_STDERR, "Exiting @ SUMA_Read_file function\n");
00545               exit (0);
00546       }
00547       else fprintf(SUMA_STDERR, "\nContinuing execution with %d points\n",cnt);
00548 
00549       }
00550 
00551    fclose (internal_file);
00552    return (cnt);                            
00553 }
00554 
00555 
00556 
00557 
00558 
00559 
00560 
00561 
00562 
00563 
00564 
00565 
00566 
00567 
00568 
00569 
00570 
00571 
00572 
00573 
00574 
00575 
00576 
00577 
00578 
00579 
00580 
00581 
00582 
00583 
00584 
00585 
00586 
00587 
00588 
00589 
00590 
00591 
00592   
00593  
00594 int SUMA_Read_2Dfile (char *f_name, float **x,  int n_cols, int n_rows)
00595 {
00596    int ir=0, ic=0, ex;
00597    FILE*internal_file;
00598    static char FuncName[]={"SUMA_Read_2Dfile"};
00599 
00600    SUMA_ENTRY;
00601 
00602    internal_file = fopen (f_name,"r");
00603    if (internal_file == NULL) {
00604                           fprintf (SUMA_STDERR,"%s: \aCould not open %s \n",FuncName, f_name);
00605                           SUMA_RETURN (-1);
00606                           }
00607    ir = 0;
00608    while (ir < n_rows)
00609    {
00610        ic = 0;
00611       while (ic < n_cols)
00612          {
00613             ex = fscanf (internal_file,"%f",&x[ir][ic]);   
00614             if (ex == EOF)
00615                {
00616                   fprintf(stderr,"Error SUMA_Read_2Dfile: Premature EOF\n");
00617                   fclose (internal_file);
00618                   SUMA_RETURN (n_rows);
00619                }
00620             ++ic;
00621          }
00622       ++ir;
00623    }
00624 
00625    fclose (internal_file);
00626    SUMA_RETURN (ir);      
00627       
00628 }
00629 
00630 
00631 
00632 
00633 
00634 
00635 SUMA_IRGB *SUMA_Create_IRGB(int n_el)
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 }
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 
00666 SUMA_IRGB *SUMA_Free_IRGB(SUMA_IRGB *irgb)
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 }
00682 
00683 
00684 
00685 
00686 
00687 
00688 
00689 
00690 
00691 
00692 SUMA_IRGB *SUMA_Read_IRGB_file (char *f_name)
00693 {
00694    int i=0, ncol = 0, nrow = 0;
00695    MRI_IMAGE *im = NULL;
00696    float *far=NULL;
00697    SUMA_IRGB *irgb=NULL;
00698    static char FuncName[]={"SUMA_Read_IRGB_file"};
00699 
00700    SUMA_ENTRY;
00701 
00702    im = mri_read_1D (f_name);
00703    
00704    if (!im) {
00705       SUMA_SLP_Err("Failed to read 1D file");
00706       SUMA_RETURN(NULL);
00707    }
00708    
00709    far = MRI_FLOAT_PTR(im);
00710    ncol = im->nx;
00711    nrow = im->ny;
00712    
00713    if (!ncol) {
00714       SUMA_SL_Err("Empty file");
00715       SUMA_RETURN(NULL);
00716    }
00717    if (nrow !=  4 ) {
00718       SUMA_SL_Err("File must have\n"
00719                   "4 columns.");
00720       mri_free(im); im = NULL;   
00721       SUMA_RETURN(NULL);
00722    }
00723   
00724    if (!(irgb = SUMA_Create_IRGB(ncol))) {
00725       fprintf (SUMA_STDERR,"%s: Failed to create irgb.\n",FuncName);
00726       SUMA_RETURN (NULL);
00727    }
00728    
00729    for (i=0; i < ncol; ++i) {
00730       irgb->i[i] = (int)far[i];
00731       irgb->r[i] = far[i+ncol];
00732       irgb->g[i] = far[i+2*ncol];
00733       irgb->b[i] = far[i+3*ncol];
00734    }   
00735    
00736    mri_free(im); im = NULL;
00737    
00738    SUMA_RETURN (irgb);      
00739       
00740 }
00741 
00742 
00743 
00744 
00745 
00746 
00747 
00748 
00749 
00750 
00751 
00752 
00753 
00754 
00755 
00756 
00757 
00758 
00759 int SUMA_Read_2Ddfile (char *f_name, int **x, int n_rows, int n_cols)
00760 {
00761    int ir, ic, ex;
00762    FILE*internal_file;
00763    static char FuncName[]={"SUMA_Read_2Ddfile"};
00764 
00765    SUMA_ENTRY;
00766 
00767    internal_file = fopen (f_name,"r");
00768    if (internal_file == NULL) {
00769       fprintf (SUMA_STDERR,"%s: \aCould not open %s \n",FuncName, f_name);
00770       SUMA_RETURN (-1);
00771    }
00772 
00773 
00774 
00775    ir = 0;
00776    while (ir < n_rows)
00777    {
00778        ic = 0;
00779       while (ic < n_cols)
00780          {
00781             ex = fscanf (internal_file,"%d",&x[ir][ic]);   
00782             if (ex == EOF)
00783                {
00784                   fprintf(stderr,"Error SUMA_Read_2Ddfile: Premature EOF\n");
00785                   fclose (internal_file);
00786                   SUMA_RETURN(ir);
00787                }
00788             ++ic;
00789          }
00790       ++ir;
00791    }
00792 
00793    fclose (internal_file);
00794    SUMA_RETURN (ir);      
00795       
00796 }
00797 
00798 
00799 
00800 
00801 
00802  
00803 int SUMA_float_file_size (char *f_name)
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 }
00828 
00829 
00830 
00831 void SUMA_alloc_problem (char *s1)
00832  
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 }
00842 
00843 
00844 
00845 
00846 
00847 
00848 
00849 
00850 
00851 
00852 
00853 
00854 
00855 
00856 
00857 
00858 
00859 
00860 
00861 
00862 
00863 
00864 
00865 
00866 char **SUMA_allocate2D (int rows,int cols,int element_size)
00867 
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    
00883    switch(element_size) {
00884      case sizeof(short): {    
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): {    
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): {   
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          
00955          
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             
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 }
00982 
00983 
00984 
00985 
00986 
00987 
00988 
00989 
00990 
00991 
00992 
00993 
00994 
00995 
00996 
00997 
00998 
00999 
01000 
01001 
01002 
01003 
01004 
01005 
01006 
01007 void SUMA_free2D(char **a,int rows)
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    
01042    for(i = 0 ; i < rows ; i++) free(a[i]);
01043 
01044    
01045    free((char *)a);
01046    a = NULL;           
01047 
01048    #ifdef USE_SUMA_MALLOC
01049       
01050       SUMA_SL_Err("NO LONGER SUPPORTED");
01051       SUMA_RETURNe;
01052 
01053    #else
01054       resume_mcw_malloc();
01055    #endif
01056    
01057    SUMA_RETURNe;
01058 }
01059 
01060 
01061 
01062 
01063 
01064 
01065 
01066 
01067 
01068 
01069 
01070 
01071 
01072 
01073 
01074 
01075 
01076 
01077 
01078    
01079 
01080 void SUMA_error_message (char *s1,char *s2,int ext)
01081  
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   }
01097 
01098 
01099 
01100 
01101 int SUMA_iswordin_ci ( const char *sbig, const char *ssub)
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 } 
01118 
01119 
01120 
01121 
01122 
01123 
01124 
01125 
01126 
01127 
01128 
01129 
01130 
01131 
01132 
01133 
01134 
01135 
01136 
01137 
01138 
01139 
01140 
01141 
01142 
01143 
01144 
01145 
01146 
01147 
01148 
01149 
01150 
01151 
01152 
01153 
01154 int SUMA_iswordin (const char *sbig, const char *ssub)
01155 {
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              
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 }
01187 
01188 
01189 
01190 
01191 
01192 
01193 
01194 
01195 
01196 
01197 
01198 
01199 
01200 
01201 
01202 
01203 
01204 
01205 
01206 
01207 
01208 
01209 
01210 
01211 
01212 
01213 
01214 
01215 
01216 
01217 
01218 
01219 
01220 
01221 
01222 
01223  
01224 void SUMA_disp_dmat (int **v,int nr, int nc , int SpcOpt)
01225 {
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 }
01248 
01249 
01250 
01251 
01252 
01253 
01254 
01255 
01256 
01257 
01258 
01259 
01260 
01261 
01262 
01263 
01264 
01265 
01266 
01267 
01268 
01269 
01270 
01271  
01272 void SUMA_disp_mat (float **v,int nr, int nc , int SpcOpt)
01273 {
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 }
01296 
01297 
01298 
01299 
01300 
01301 
01302 
01303 
01304 
01305 
01306 
01307 
01308 
01309 
01310 
01311 
01312 
01313 
01314 
01315 
01316 
01317 
01318 
01319 
01320 
01321 
01322 
01323 
01324 
01325 
01326  
01327 void SUMA_disp_vecmat (float *v,int nr, int nc , int SpcOpt, 
01328                         SUMA_INDEXING_ORDER d_order, FILE *fout, SUMA_Boolean AddRowInd)
01329 {
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"); 
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 }
01374 
01375 
01376 
01377 
01378 
01379 
01380 
01381 
01382 
01383 
01384 
01385 
01386 
01387 
01388 
01389 
01390 
01391 
01392 
01393 
01394 
01395 
01396 
01397 
01398 
01399 
01400 
01401 
01402 
01403 
01404 
01405 
01406  
01407 void SUMA_disp_vecdmat (int *v,int nr, int nc , int SpcOpt, 
01408                         SUMA_INDEXING_ORDER d_order, FILE *fout, SUMA_Boolean AddRowInd)
01409 {
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"); 
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 }
01450 
01451 
01452 
01453 
01454 void SUMA_disp_vecucmat (unsigned char *v,int nr, int nc , int SpcOpt, 
01455                         SUMA_INDEXING_ORDER d_order, FILE *fout, SUMA_Boolean AddRowInd)
01456 {
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"); 
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 }
01497 
01498 
01499 
01500 
01501 
01502 
01503 
01504 
01505 
01506 
01507 
01508 
01509 
01510 
01511 
01512 
01513 
01514 
01515 
01516 void SUMA_disp_vect (float *v,int l)
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 }
01533 
01534 
01535 
01536 
01537 
01538 
01539 
01540 
01541 
01542 
01543 
01544 
01545 
01546 
01547 
01548 
01549 
01550 
01551 
01552 void SUMA_disp_dvect (int *v,int l)
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 }
01570 
01571 
01572 
01573 
01574 
01575 
01576 
01577 
01578 
01579 
01580 
01581 
01582 
01583 
01584 
01585 
01586 
01587 
01588 
01589 
01590 
01591 
01592 
01593 
01594 
01595 
01596 
01597 
01598 
01599 
01600 
01601 
01602 
01603 
01604 
01605 
01606 
01607 
01608 
01609 
01610 
01611 
01612 
01613 
01614 
01615 float SUMA_etime (struct  timeval  *t, int Report  )
01616 {
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    
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 }
01641 
01642   
01643 
01644 
01645 
01646 
01647 
01648 
01649 
01650 
01651 
01652 
01653 
01654 
01655 
01656 
01657 
01658 
01659 
01660 
01661 
01662 
01663 
01664 
01665 
01666 
01667 
01668 
01669 
01670 
01671 
01672 
01673 
01674 
01675 
01676 
01677 
01678 
01679 
01680 
01681 
01682 
01683 SUMA_ISINSPHERE SUMA_isinsphere (float * NodeList, int nr, float *S_cent , float S_rad , int BoundIn )
01684 {
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) 
01706       {
01707          for (k=0; k < nr; ++k)
01708             {
01709                id = ND * k;
01710                
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                
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    
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 }
01766 
01767 
01768 
01769 
01770 SUMA_Boolean SUMA_Free_IsInSphere (SUMA_ISINSPHERE *IB)
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 }
01785 
01786 
01787 
01788 
01789 
01790 
01791 
01792 
01793 
01794 
01795 
01796 
01797 
01798 
01799 
01800 
01801 
01802 
01803 
01804 
01805 
01806 
01807 
01808 
01809 
01810 
01811 
01812 
01813 
01814 
01815 
01816 
01817 
01818 
01819 
01820 
01821 
01822 
01823 
01824 
01825 
01826 
01827 
01828 SUMA_ISINBOX SUMA_isinbox (float * XYZ, int nr, float *S_cent , float *S_dim , int BoundIn )
01829 {
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 
01841 
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) 
01860       {
01861          
01862          for (k=0; k < nr; ++k)
01863             {
01864             
01865             
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             
01883 
01884       }
01885    else
01886       {
01887          for (k=0; k < nr; ++k)
01888             {
01889                
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       
01910 
01911       
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       
01926       IsIn_strct.IsIn = NULL;
01927       IsIn_strct.d = NULL;
01928    }
01929    
01930    
01931    SUMA_free(IsIn);
01932    SUMA_free(d);
01933    
01934 
01935    SUMA_RETURN (IsIn_strct) ;
01936 
01937 }
01938 
01939 
01940 
01941 
01942 
01943 SUMA_Boolean SUMA_Free_IsInBox (SUMA_ISINBOX *IB)
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 }
01958 
01959 
01960 
01961 
01962 
01963 
01964 
01965 
01966 
01967 
01968 
01969 
01970 
01971 
01972 
01973 
01974 
01975 
01976 
01977 
01978 
01979 
01980 
01981 byte * SUMA_isinpoly(float *P, float *NodeList, int *FaceSetList, int N_FaceSet, int FaceSetDim, int *dims, int *N_in, byte *usethis, byte *culled)
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]]; 
02010    for (iv = 0; iv < N_FaceSet; ++iv) {
02011       counter = 0;
02012       for (i=0; i<FaceSetDim; ++i) { 
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]]; 
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); 
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 }
02056 
02057        
02058 
02059 
02060 
02061 
02062 
02063 
02064 
02065 
02066 
02067 
02068 
02069 
02070 
02071 
02072 
02073 
02074 
02075 
02076 
02077 
02078 
02079 
02080 
02081 
02082 
02083 float **SUMA_Point_At_Distance(float *U, float *P1, float d)
02084 {
02085    static char FuncName[]={"SUMA_Point_At_Distance"}; 
02086    float bf, **P2, P1orig[3], Uorig[3];
02087    float m, n, p, q, D, A, B, C, epsi = 0.0001;
02088    int flip, i;
02089    SUMA_Boolean LocalHead = NOPE;
02090    
02091    SUMA_ENTRY;
02092 
02093    SUMA_SL_Warn ("useless piece of junk, use SUMA_POINT_AT_DISTANCE instead!");
02094    
02095    if (d == 0) {
02096       fprintf(SUMA_STDERR,"Error %s: d is 0. Not good, Not good at all.\n", FuncName);
02097       SUMA_RETURN (NULL);
02098    }
02099    
02100    if (LocalHead) {
02101       fprintf (SUMA_STDOUT,"%s: U %f, %f, %f, P1 %f %f %f, d %f\n", FuncName,\
02102          U[0], U[1], U[2], P1[0], P1[1], P1[2], d);
02103    }
02104          
02105    
02106    P1orig[0] = P1[0];    
02107    P1orig[1] = P1[1]; 
02108    P1orig[2] = P1[2]; 
02109 
02110    Uorig[0] = U[0];
02111    Uorig[1] = U[1];
02112    Uorig[2] = U[2];
02113    
02114    
02115    flip = 0;
02116    if (fabs(U[0]) < epsi) { 
02117       if (fabs(U[1]) > epsi) {
02118          U[0] = U[1]; U[1] = 0;
02119          bf = P1[0]; P1[0] = P1[1]; P1[1] = bf;
02120          flip = 1;
02121       } else {   
02122          if (fabs(U[2]) > epsi) { 
02123             U[0] = U[2]; U[2] = 0;
02124             bf = P1[0]; P1[0] = P1[2]; P1[2] = bf;
02125             flip = 2;
02126          } else { 
02127             fprintf(SUMA_STDERR, "Error %s: 0 direction vector.\n", FuncName);
02128             SUMA_RETURN (NULL);
02129          }
02130       }
02131    }
02132 
02133    if (LocalHead) fprintf (SUMA_STDERR, "%s: flip = %d\n", FuncName, flip);
02134       
02135    if (LocalHead) fprintf (SUMA_STDERR, "%s: U original: %f, %f, %f\n", FuncName, U[0], U[1], U[2]);
02136    U[1] /= U[0];
02137    U[2] /= U[0];
02138    U[0] = 1.0; 
02139    if (LocalHead) fprintf (SUMA_STDERR, "%s: U normalized: %f, %f, %f\n", FuncName, U[0], U[1], U[2]);
02140    
02141       
02142    m = U[1];
02143    n = U[2];
02144 
02145    q = P1[1] - m*P1[0];
02146    p = P1[2] - n*P1[0];
02147 
02148    if (LocalHead) fprintf (SUMA_STDERR, "%s: m=%f n=%f, p=%f, q=%f\n", FuncName, m, n, p, q);
02149 
02150    
02151    A = (1 + n*n + m*m);
02152    B = -2 * P1[0] + 2 * m * (q - P1[1]) + 2 * n * (p - P1[2]);
02153    C = P1[0]*P1[0] + (q - P1[1])*(q - P1[1]) + (p - P1[2])*(p - P1[2]) - d*d;
02154 
02155    D = B*B - 4*A*C;
02156    
02157    if (LocalHead) fprintf (SUMA_STDERR, "%s: A=%f B=%f, C=%f, D=%f\n", FuncName, A, B, C, D);
02158    
02159    if (D < 0) {
02160       fprintf(SUMA_STDERR, "Error %s: Negative Delta: %f.\n"
02161                            "Input values were: \n"
02162                            "U :[%f %f %f]\n"
02163                            "P1:[%f %f %f]\n"
02164                            "d :[%f]\n"
02165                            , FuncName, D, Uorig[0], Uorig[1], Uorig[2], 
02166                            P1orig[0], P1orig[1], P1orig[2], d);
02167       SUMA_RETURN(NULL);
02168    }
02169 
02170    P2 = (float **)SUMA_allocate2D(2,3, sizeof(float));
02171    if (P2 == NULL) {
02172       fprintf(SUMA_STDERR, "Error %s: Could not allocate for 6 floats! What is this? What is the matter with you?!\n", FuncName);
02173       SUMA_RETURN (NULL);
02174    }
02175 
02176    P2[0][0] = (-B + sqrt(D)) / (2 *A);
02177    P2[1][0] = (-B - sqrt(D)) / (2 *A);
02178 
02179    P2[0][1] = m * P2[0][0] + q;
02180    P2[1][1] = m * P2[1][0] + q;
02181 
02182    P2[0][2] = n * P2[0][0] + p;
02183    P2[1][2] = n * P2[1][0] + p;
02184 
02185 
02186    
02187    if (flip == 1) {
02188     for (i=0; i < 2; ++i) {
02189        bf = P2[i][1];
02190        P2[i][1] = P2[i][0];
02191        P2[i][0] = bf;
02192       }
02193    } else if (flip == 2){
02194     for (i=0; i < 2; ++i) {
02195        bf = P2[i][2]; 
02196        P2[i][2] = P2[i][0];
02197        P2[i][0] = bf;
02198       }
02199    }   
02200 
02201    for (i=0; i < 3; ++i) {
02202       P1[i] = P1orig[i];
02203       U[i] = Uorig[i];
02204    }
02205 
02206    if (LocalHead) {
02207       fprintf(SUMA_STDOUT,"%s: P1 = %f, %f, %f\n  ", \
02208        FuncName, P1[0], P1[1], P1[2]);
02209       fprintf(SUMA_STDOUT,"%s: P2 = %f, %f, %f\n    %f, %f, %f\n", \
02210        FuncName, P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]);
02211       fprintf(SUMA_STDOUT,"%s: U = %f, %f, %f\n  ", \
02212        FuncName, U[0], U[1], U[2]);
02213    }
02214 
02215    
02216    
02217    Uorig[0] = P2[0][0] - P1[0]; 
02218    Uorig[1] = P2[0][1] - P1[1];
02219    Uorig[2] = P2[0][2] - P1[2];
02220 
02221    SUMA_DOTP_VEC(Uorig, U, bf, 3, float, float)
02222    if (LocalHead) fprintf(SUMA_STDOUT,"%s: Dot product = %f\n", FuncName, bf);
02223    if (bf < 0) {
02224       if (LocalHead) fprintf(SUMA_STDOUT,"%s: Flipping at end...\n", FuncName);
02225       for (i=0; i< 3; ++i) {
02226          bf = P2[0][i];
02227          P2[0][i] = P2[1][i]; P2[1][i] = bf;
02228       }
02229    }
02230    
02231    if (LocalHead) {
02232       fprintf(SUMA_STDOUT,"%s: P2 = %f, %f, %f\n    %f, %f, %f\n", \
02233        FuncName, P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]);
02234    }
02235 SUMA_RETURN (P2);
02236    
02237 }
02238 double **SUMA_dPoint_At_Distance(double *U, double *P1, double d)
02239 {
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    
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    
02270    flip = 0;
02271    if (fabs(U[0]) < epsi) { 
02272       if (fabs(U[1]) > epsi) {
02273          U[0] = U[1]; U[1] = 0;
02274          bf = P1[0]; P1[0] = P1[1]; P1[1] = bf;
02275          flip = 1;
02276       } else {   
02277          if (fabs(U[2]) > epsi) { 
02278             U[0] = U[2]; U[2] = 0;
02279             bf = P1[0]; P1[0] = P1[2]; P1[2] = bf;
02280             flip = 2;
02281          } else { 
02282             fprintf(SUMA_STDERR, "Error %s: 0 direction vector.\n", FuncName);
02283             SUMA_RETURN (NULL);
02284          }
02285       }
02286    }
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       
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    
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    
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    
02372    Uorig[0] = P2[0][0] - P1[0]; 
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 }
02393 
02394 
02395 
02396 
02397 
02398 
02399 
02400 
02401 
02402 
02403 
02404 
02405 
02406 
02407 
02408 
02409 
02410 
02411 
02412 
02413 
02414 
02415 
02416 
02417 
02418 
02419 
02420 SUMA_Boolean SUMA_Point_To_Line_Distance (float *NodeList, int N_points, float *P1, float *P2, float *d2, float *d2min, int *i2min)
02421 {
02422    static char FuncName[]={"SUMA_Point_To_Line_Distance"};
02423    float U[3], Un, xn, yn, zn, dx, dy, dz;
02424    int i, id, ND;
02425    
02426    SUMA_ENTRY;
02427    
02428    ND = 3;
02429    if (N_points < 1) {
02430       fprintf(SUMA_STDERR,"Error %s: N_points is 0.\n",FuncName);
02431       SUMA_RETURN (NOPE);
02432    }
02433    
02434    SUMA_UNIT_VEC(P1, P2, U, Un);
02435    if (Un == 0) {
02436       fprintf(SUMA_STDERR,"Error %s: P1 and P2 are identical.\n",FuncName);
02437       SUMA_RETURN (NOPE);
02438    }
02439    
02440    
02441    
02442    
02443    
02444    
02445    
02446 
02447    
02448    if (d2 == NULL) {
02449       fprintf(SUMA_STDERR,"Error %s: d2 not allocated for.\n",FuncName);
02450       SUMA_RETURN (NOPE);
02451    }
02452    
02453    
02454    
02455     i = 0;
02456     xn = NodeList[0] - P1[0];
02457     yn = NodeList[1] - P1[1];
02458     zn = NodeList[2] - P1[2];
02459     
02460     dx = (U[1]*zn - yn*U[2]);
02461     dy = (U[0]*zn - xn*U[2]);
02462     dz = (U[0]*yn - xn*U[1]);
02463     
02464     d2[i] = dx*dx+dy*dy +dz*dz; 
02465     *d2min = d2[i];
02466     *i2min = i;
02467     
02468    for (i=1; i < N_points; ++i) {
02469       id = ND * i;
02470       xn = NodeList[id] - P1[0];
02471       yn = NodeList[id+1] - P1[1];
02472       zn = NodeList[id+2] - P1[2];
02473 
02474       dx = (U[1]*zn - yn*U[2]);
02475       dy = (U[0]*zn - xn*U[2]);
02476       dz = (U[0]*yn - xn*U[1]);
02477 
02478       d2[i] = dx*dx+dy*dy +dz*dz; 
02479       if (d2[i] < *d2min) {
02480          *d2min = d2[i];
02481          *i2min = i;
02482       }
02483    }
02484    SUMA_RETURN (YUP);
02485 }
02486 
02487 
02488 
02489 
02490 
02491 
02492 
02493 
02494 
02495 
02496 
02497 
02498 
02499 
02500 
02501 
02502 
02503 
02504 
02505 
02506 
02507 
02508 
02509 
02510 SUMA_Boolean SUMA_Point_To_Point_Distance (float *NodeList, int N_points, float *P1, float *d2, float *d2min, int *i2min)
02511 {
02512    static char FuncName[]={"SUMA_Point_To_Point_Distance"};
02513    float xn, yn, zn;
02514    int i, id, ND;
02515    
02516    SUMA_ENTRY;
02517 
02518    ND = 3;
02519    if (N_points < 1) {
02520       fprintf(SUMA_STDERR,"Error %s: N_points is 0.\n",FuncName);
02521       SUMA_RETURN (NOPE);
02522    }
02523    
02524    
02525    
02526    
02527    if (d2 == NULL) {
02528       fprintf(SUMA_STDERR,"Error %s: d2 not allocated for.\n",FuncName);
02529       SUMA_RETURN (NOPE);
02530    }
02531    
02532    
02533    
02534     i = 0;
02535     xn = NodeList[0] - P1[0];
02536     yn = NodeList[1] - P1[1];
02537     zn = NodeList[2] - P1[2];
02538     
02539     d2[i] = xn*xn + yn*yn + zn*zn; 
02540     *d2min = d2[i];
02541     *i2min = i;
02542     
02543    for (i=1; i < N_points; ++i) {
02544       id = ND * i;
02545       xn = NodeList[id] - P1[0];
02546       yn = NodeList[id+1] - P1[1];
02547       zn = NodeList[id+2] - P1[2];
02548 
02549 
02550        d2[i] = xn*xn + yn*yn + zn*zn; 
02551       if (d2[i] < *d2min) {
02552          *d2min = d2[i];
02553          *i2min = i;
02554       }
02555    }
02556    SUMA_RETURN (YUP);
02557 }
02558 
02559 
02560 
02561 #define SUMA_Z_QSORT_structs
02562 
02563 
02564 
02565    typedef struct {
02566       float x;
02567       int Index;
02568    } SUMA_Z_QSORT_FLOAT;
02569 
02570    typedef struct {
02571       double x;
02572       int Index;
02573    } SUMA_Z_QSORT_DOUBLE;
02574 
02575    typedef struct {
02576       int x;
02577       int Index;
02578    } SUMA_Z_QSORT_INT;
02579 
02580 int compare_SUMA_Z_QSORT_FLOAT (SUMA_Z_QSORT_FLOAT *a, SUMA_Z_QSORT_FLOAT *b )
02581    {
02582       if (a->x < b->x)
02583          return (-1);
02584       else if (a->x == b->x)
02585          return (0);
02586       else if (a->x > b->x)
02587          return (1);
02588       
02589       return (0);
02590    }
02591 
02592 int compare_SUMA_Z_QSORT_DOUBLE (SUMA_Z_QSORT_DOUBLE *a, SUMA_Z_QSORT_DOUBLE *b )
02593    {
02594       if (a->x < b->x)
02595          return (-1);
02596       else if (a->x == b->x)
02597          return (0);
02598       else if (a->x > b->x)
02599          return (1);
02600       
02601       return (0);
02602    }
02603    
02604    
02605 int compare_SUMA_Z_QSORT_INT (SUMA_Z_QSORT_INT *a, SUMA_Z_QSORT_INT *b )
02606    {
02607       if (a->x < b->x)
02608          return (-1);
02609       else if (a->x == b->x)
02610          return (0);
02611       else if (a->x > b->x)
02612          return (1);
02613       
02614       return (0);
02615    }
02616     
02617 
02618 int SUMA_compare_int (int *a, int *b )
02619 {
02620     if (*a < *b)
02621       return (-1);
02622    else if (*a == *b)
02623       return (0);
02624    else
02625       return (1);
02626    
02627 }
02628    
02629 
02630 
02631 
02632 
02633 
02634 
02635 
02636 
02637 
02638 
02639 
02640 
02641 
02642 
02643 
02644 
02645 
02646 
02647 
02648 
02649 
02650 
02651 
02652 
02653 
02654 
02655 
02656 
02657 
02658 
02659 
02660 
02661 
02662 
02663 
02664 
02665 
02666 
02667 int *SUMA_z_qsort (float *x , int nx )
02668 {
02669    static char FuncName[]={"SUMA_z_qsort"}; 
02670    int *I, k;
02671    SUMA_Z_QSORT_FLOAT *Z_Q_fStrct;
02672    
02673    SUMA_ENTRY;
02674 
02675    
02676    Z_Q_fStrct = (SUMA_Z_QSORT_FLOAT *) SUMA_calloc(nx, sizeof (SUMA_Z_QSORT_FLOAT));
02677    I = (int *) SUMA_calloc (nx, sizeof(int));
02678 
02679    if (!Z_Q_fStrct || !I)
02680       {
02681          fprintf(SUMA_STDERR,"Error %s: Allocation problem.\n",FuncName);
02682          SUMA_RETURN (NULL);
02683       }
02684 
02685    for (k=0; k < nx; ++k) 
02686       {
02687          Z_Q_fStrct[k].x = x[k];
02688          Z_Q_fStrct[k].Index = k;
02689       }
02690 
02691    
02692    qsort(Z_Q_fStrct, nx, sizeof(SUMA_Z_QSORT_FLOAT), (int(*) (const void *, const void *)) compare_SUMA_Z_QSORT_FLOAT);
02693 
02694    
02695    for (k=0; k < nx; ++k) 
02696       {
02697          x[k] = Z_Q_fStrct[k].x;
02698          I[k] = Z_Q_fStrct[k].Index;
02699       }
02700 
02701    
02702    SUMA_free(Z_Q_fStrct);
02703 
02704    
02705    SUMA_RETURN (I);
02706 
02707 
02708 }
02709 
02710 int *SUMA_z_doubqsort (double *x , int nx )
02711 {
02712    static char FuncName[]={"SUMA_z_doubqsort"}; 
02713    int *I, k;
02714    SUMA_Z_QSORT_DOUBLE *Z_Q_doubStrct;
02715    
02716    SUMA_ENTRY;
02717 
02718    
02719    Z_Q_doubStrct = (SUMA_Z_QSORT_DOUBLE *) SUMA_calloc(nx, sizeof (SUMA_Z_QSORT_DOUBLE));
02720    I = (int *) SUMA_calloc (nx, sizeof(int));
02721 
02722    if (!Z_Q_doubStrct || !I)
02723       {
02724          fprintf(SUMA_STDERR,"Error %s: Allocation problem.\n",FuncName);
02725          SUMA_RETURN (NULL);
02726       }
02727 
02728    for (k=0; k < nx; ++k) 
02729       {
02730          Z_Q_doubStrct[k].x = x[k];
02731          Z_Q_doubStrct[k].Index = k;
02732       }
02733 
02734    
02735    qsort(Z_Q_doubStrct, nx, sizeof(SUMA_Z_QSORT_DOUBLE), (int(*) (const void *, const void *)) compare_SUMA_Z_QSORT_DOUBLE);
02736 
02737    
02738    for (k=0; k < nx; ++k) 
02739       {
02740          x[k] = Z_Q_doubStrct[k].x;
02741          I[k] = Z_Q_doubStrct[k].Index;
02742       }
02743 
02744    
02745    SUMA_free(Z_Q_doubStrct);
02746 
02747    
02748    SUMA_RETURN (I);
02749 
02750 
02751 }
02752 
02753 int *SUMA_z_dqsort (int *x , int nx )
02754 {
02755    static char FuncName[]={"SUMA_z_dqsort"}; 
02756    int *I, k;
02757    SUMA_Z_QSORT_INT *Z_Q_iStrct;
02758    
02759    SUMA_ENTRY;
02760 
02761    
02762 
02763    Z_Q_iStrct = (SUMA_Z_QSORT_INT *) SUMA_calloc(nx, sizeof (SUMA_Z_QSORT_INT));
02764    I = (int *) SUMA_calloc (nx,sizeof(int));
02765 
02766    if (!Z_Q_iStrct || !I)
02767       {
02768          fprintf(SUMA_STDERR,"Error %s: Allocation problem.\n",FuncName);
02769          SUMA_RETURN (NULL);
02770       }
02771 
02772    for (k=0; k < nx; ++k) 
02773       {
02774          Z_Q_iStrct[k].x = x[k];
02775          Z_Q_iStrct[k].Index = k;
02776       }
02777 
02778    
02779    qsort(Z_Q_iStrct, nx, sizeof(SUMA_Z_QSORT_INT), (int(*) (const void *, const void *)) compare_SUMA_Z_QSORT_INT);
02780 
02781    
02782    for (k=0; k < nx; ++k) 
02783       {
02784          x[k] = Z_Q_iStrct[k].x;
02785          I[k] = Z_Q_iStrct[k].Index;
02786       }
02787 
02788    
02789    SUMA_free(Z_Q_iStrct);
02790 
02791    
02792    SUMA_RETURN (I);
02793 
02794       
02795 }
02796    
02797 
02798 
02799 
02800 int *SUMA_z_dqsort_nsc (int *x , int nx )
02801 {
02802    static char FuncName[]={"SUMA_z_dqsort_nsc"}; 
02803    int *I, k;
02804    SUMA_Z_QSORT_INT *Z_Q_iStrct;
02805    
02806    SUMA_ENTRY;
02807 
02808    
02809 
02810    Z_Q_iStrct = (SUMA_Z_QSORT_INT *) calloc(nx, sizeof (SUMA_Z_QSORT_INT));
02811    I = (int *) calloc (nx,sizeof(int));
02812 
02813    if (!Z_Q_iStrct || !I)
02814       {
02815          fprintf(SUMA_STDERR,"Error %s: Allocation problem.\n",FuncName);
02816          SUMA_RETURN (NULL);
02817       }
02818 
02819    for (k=0; k < nx; ++k) 
02820       {
02821          Z_Q_iStrct[k].x = x[k];
02822          Z_Q_iStrct[k].Index = k;
02823       }
02824 
02825    
02826    qsort(Z_Q_iStrct, nx, sizeof(SUMA_Z_QSORT_INT), (int(*) (const void *, const void *)) compare_SUMA_Z_QSORT_INT);
02827 
02828    
02829    for (k=0; k < nx; ++k) 
02830       {
02831          x[k] = Z_Q_iStrct[k].x;
02832          I[k] = Z_Q_iStrct[k].Index;
02833       }
02834 
02835    
02836    free(Z_Q_iStrct);
02837 
02838    
02839    SUMA_RETURN (I);
02840 
02841       
02842 }
02843    
02844    
02845 
02846    
02847 typedef struct {
02848       float *x;
02849       int ncol;
02850       int Index;
02851    } SUMA_QSORTROW_FLOAT;
02852 
02853 
02854 
02855 int compare_SUMA_QSORTROW_FLOAT (SUMA_QSORTROW_FLOAT *a, SUMA_QSORTROW_FLOAT *b)
02856    {
02857       int k;
02858       
02859       for (k=0; k < a->ncol ; ++k)
02860          {
02861             if (a->x[k] < b->x[k])
02862                return (-1);
02863             else if (a->x[k] > b->x[k])
02864                return (1);
02865          }
02866       return (0); 
02867    }
02868 
02869    
02870 
02871 
02872 
02873 
02874 
02875 
02876 
02877 
02878 
02879 
02880 
02881 
02882 
02883 
02884 
02885 
02886 
02887 
02888 
02889 
02890 int * SUMA_fqsortrow (float **X , int nr, int nc  )
02891 {
02892    static char FuncName[]={"SUMA_fqsortrow"}; 
02893    int k, *I;
02894    SUMA_QSORTROW_FLOAT *Z_Q_fStrct;
02895    
02896       
02897    SUMA_ENTRY;
02898 
02899    
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) 
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    
02917    qsort(Z_Q_fStrct, nr, sizeof(SUMA_QSORTROW_FLOAT), (int(*) (const void *, const void *)) compare_SUMA_QSORTROW_FLOAT);
02918 
02919    
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    
02927    SUMA_free(Z_Q_fStrct);
02928 
02929    
02930    SUMA_RETURN (I);
02931    
02932    
02933 }
02934 
02935    typedef struct {
02936       int *x;
02937       int ncol;
02938       int Index;
02939    } SUMA_QSORTROW_INT;
02940 
02941 
02942 
02943 int compare_SUMA_QSORTROW_INT (SUMA_QSORTROW_INT *a, SUMA_QSORTROW_INT *b)
02944    {
02945       int k;
02946       
02947       for (k=0; k < a->ncol ; ++k)
02948          {
02949             if (a->x[k] < b->x[k])
02950                return (-1);
02951             else if (a->x[k] > b->x[k])
02952                return (1);
02953          }
02954       return (0); 
02955    }
02956 
02957    
02958 
02959 
02960 
02961 
02962 
02963 
02964 
02965 
02966 
02967 
02968 
02969 
02970 
02971 
02972 
02973 
02974 
02975 
02976 
02977 
02978 
02979 int * SUMA_dqsortrow (int **X , int nr, int nc  )
02980 {
02981    static char FuncName[]={"SUMA_dqsortrow"}; 
02982    int k,  *I;
02983    SUMA_QSORTROW_INT *Z_Q_dStrct;
02984    
02985    SUMA_ENTRY;
02986    
02987    
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) 
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    
03005    qsort(Z_Q_dStrct, nr, sizeof(SUMA_QSORTROW_INT), (int(*) (const void *, const void *)) compare_SUMA_QSORTROW_INT);
03006 
03007    
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    
03015    SUMA_free(Z_Q_dStrct);
03016 
03017    
03018    SUMA_RETURN (I);
03019    
03020    
03021 }
03022 
03023 
03024 
03025 
03026 
03027 
03028 
03029 
03030 
03031 
03032 
03033 #define SUMA_CORNER_OF_VOXEL(center, dxyz, cn, P){\
03034    switch(cn) {   \
03035       case 0:  \
03036          P[0] =  dxyz[0]; P[1] = -dxyz[1]; P[2] =   dxyz[2];   \
03037          break;   \
03038       case 1:  \
03039          P[0] =  dxyz[0]; P[1] = -dxyz[1]; P[2] =  -dxyz[2];   \
03040          break;   \
03041       case 2:  \
03042          P[0] =  dxyz[0]; P[1] =  dxyz[1]; P[2] =  -dxyz[2];   \
03043          break;  \
03044       case 3:  \
03045          P[0] =  dxyz[0]; P[1] =  dxyz[1]; P[2] =   dxyz[2];   \
03046          break;   \
03047       case 4:  \
03048          P[0] = -dxyz[0]; P[1] = -dxyz[1]; P[2] =   dxyz[2];   \
03049          break;   \
03050       case 5:  \
03051          P[0] = -dxyz[0]; P[1] = -dxyz[1]; P[2] =  -dxyz[2];   \
03052          break;   \
03053       case 6:  \
03054          P[0] = -dxyz[0]; P[1] =  dxyz[1]; P[2] =  -dxyz[2];   \
03055          break;  \
03056       case 7:  \
03057          P[0] = -dxyz[0]; P[1] =  dxyz[1]; P[2] =   dxyz[2];   \
03058          break;  \
03059       default: \
03060          P[0] = 0; P[1] = 0; P[2] = 0; \
03061          SUMA_SL_Err("Bad corner index. Returning Center of voxel.");\
03062    }  \
03063      \
03064    P[0] = center[0] + 0.5 * P[0];   \
03065    P[1] = center[1] + 0.5 * P[1];   \
03066    P[2] = center[2] + 0.5 * P[2];   \
03067 }
03068 
03069 
03070 
03071 
03072 
03073 
03074 
03075 
03076 #define SUMA_EDGE_OF_VOXEL(center, dxyz, en, P0, P1){ \
03077    switch(en) {   \
03078       case 0:  \
03079          SUMA_CORNER_OF_VOXEL(center, dxyz, 0, P0);   \
03080          SUMA_CORNER_OF_VOXEL(center, dxyz, 1, P1);   \
03081          break;   \
03082       case 1:  \
03083          SUMA_CORNER_OF_VOXEL(center, dxyz, 0, P0);   \
03084          SUMA_CORNER_OF_VOXEL(center, dxyz, 3, P1);   \
03085          break;   \
03086       case 2:  \
03087          SUMA_CORNER_OF_VOXEL(center, dxyz, 0, P0);   \
03088          SUMA_CORNER_OF_VOXEL(center, dxyz, 4, P1);   \
03089          break;   \
03090       case 3:  \
03091          SUMA_CORNER_OF_VOXEL(center, dxyz, 1, P0);   \
03092          SUMA_CORNER_OF_VOXEL(center, dxyz, 5, P1);   \
03093          break;   \
03094       case 4:  \
03095          SUMA_CORNER_OF_VOXEL(center, dxyz, 1, P0);   \
03096          SUMA_CORNER_OF_VOXEL(center, dxyz, 2, P1);   \
03097          break;   \
03098       case 5:  \
03099          SUMA_CORNER_OF_VOXEL(center, dxyz, 2, P0);   \
03100          SUMA_CORNER_OF_VOXEL(center, dxyz, 6, P1);   \
03101          break;   \
03102       case 6:  \
03103          SUMA_CORNER_OF_VOXEL(center, dxyz, 2, P0);   \
03104          SUMA_CORNER_OF_VOXEL(center, dxyz, 3, P1);   \
03105          break;   \
03106       case 7:  \
03107          SUMA_CORNER_OF_VOXEL(center, dxyz, 3, P0);   \
03108          SUMA_CORNER_OF_VOXEL(center, dxyz, 7, P1);   \
03109          break;   \
03110       case 8:  \
03111          SUMA_CORNER_OF_VOXEL(center, dxyz, 4, P0);   \
03112          SUMA_CORNER_OF_VOXEL(center, dxyz, 7, P1);   \
03113          break;   \
03114       case 9:  \
03115          SUMA_CORNER_OF_VOXEL(center, dxyz, 4, P0);   \
03116          SUMA_CORNER_OF_VOXEL(center, dxyz, 5, P1);   \
03117          break;   \
03118       case 10:  \
03119          SUMA_CORNER_OF_VOXEL(center, dxyz, 5, P0);   \
03120          SUMA_CORNER_OF_VOXEL(center, dxyz, 6, P1);   \
03121          break;   \
03122       case 11:  \
03123          SUMA_CORNER_OF_VOXEL(center, dxyz, 6, P0);   \
03124          SUMA_CORNER_OF_VOXEL(center, dxyz, 7, P1);   \
03125          break;   \
03126       default: \
03127          P0[0] = center[0]; P0[1] = center[1]; P0[2] = center[2]; \
03128          P1[0] = center[0]; P1[1] = center[1]; P1[2] = center[2]; \
03129          SUMA_SL_Err("Bad edge index. Returning Centers of voxel.");\
03130    }  \
03131 }
03132 
03133 
03134 
03135 
03136 
03137 
03138 
03139 
03140 
03141 #define SUMA_IS_POINT_IN_SEGMENT(p, p0, p1)  (  (  (  \
03142                                                    (p[0] >  p0[0] && p[0] <  p1[0]) ||   \
03143                                                    (p[0] <  p0[0] && p[0] >  p1[0]) ||   \
03144                                                    (p[0] == p0[0] || p[0] == p1[0]) ) \
03145                                                    && \
03146                                                    (  \
03147                                                    (p[1] >  p0[1] && p[1] <  p1[1]) ||   \
03148                                                    (p[1] <  p0[1] && p[1] >  p1[1]) ||   \
03149                                                    (p[1] == p0[1] || p[1] == p1[1]) ) \
03150                                                    && \
03151                                                    (  \
03152                                                    (p[2] >  p0[2] && p[2] <  p1[2]) ||   \
03153                                                    (p[2] <  p0[2] && p[2] >  p1[2]) ||   \
03154                                                    (p[2] == p0[2] || p[2] == p1[2]) ) \
03155                                                    ) ? 1 : 0 )
03156                                           
03157 
03158 
03159 
03160 
03161 
03162 
03163 
03164 
03165 
03166 SUMA_Boolean SUMA_isVoxelIntersect_Triangle (float *center, float *dxyz, float *vert0, float *vert1, float *vert2)   
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    
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          
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 }
03197 
03198 
03199 
03200 
03201 
03202 
03203 
03204 
03205 
03206 
03207 
03208 
03209 
03210 
03211 
03212 
03213 
03214 
03215 
03216 
03217 
03218 
03219 
03220 
03221 
03222 
03223 SUMA_Boolean SUMA_MT_isIntersect_Triangle (float *P0, float *P1, float *vert0, float *vert1, float *vert2, float *iP, float *d, int *closest_vert)
03224 {  
03225    static char FuncName[]={"SUMA_MT_isIntersect_Triangle"};
03226    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
03227    double det,inv_det, u, v, t;
03228    double dir[3], dirn, orig[3];
03229    SUMA_Boolean hit = NOPE;
03230    
03231    SUMA_ENTRY;
03232    
03233    
03234    orig[0] = (double)P0[0];
03235    orig[1] = (double)P0[1];
03236    orig[2] = (double)P0[2];
03237    
03238    dir[0] = (double)P1[0] - orig[0];
03239    dir[1] = (double)P1[1] - orig[1];
03240    dir[2] = (double)P1[2] - orig[2];
03241    dirn = sqrt(dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]);
03242    dir[0] /= dirn;
03243    dir[1] /= dirn;
03244    dir[2] /= dirn;
03245 
03246    
03247    SUMA_MT_SUB(edge1, vert1, vert0);
03248    SUMA_MT_SUB(edge2, vert2, vert0);
03249 
03250    
03251    SUMA_MT_CROSS(pvec, dir, edge2);
03252 
03253    
03254    det = SUMA_MT_DOT(edge1, pvec);
03255    
03256    hit = NOPE;
03257    
03258       if (det > -SUMA_EPSILON && det < SUMA_EPSILON) {
03259          
03260          hit = NOPE;
03261       } else {
03262          inv_det = 1.0 / det;
03263 
03264          
03265          SUMA_MT_SUB(tvec, orig, vert0);
03266 
03267          
03268          u = SUMA_MT_DOT(tvec, pvec) * inv_det;
03269          if (u < 0.0 || u > 1.0) {
03270             
03271             hit = NOPE;
03272          } else {
03273             
03274             SUMA_MT_CROSS(qvec, tvec, edge1);
03275 
03276             
03277             v = SUMA_MT_DOT(dir, qvec) * inv_det;
03278             if (v < 0.0 || u + v > 1.0) {
03279                 
03280                 hit = NOPE;
03281             } else {
03282                hit = YUP;
03283                
03284                if (iP) {
03285                   
03286                   t = SUMA_MT_DOT(edge2, qvec) * inv_det;         
03287 
03288                   
03289                   iP[0] = vert0[0] + u * (vert1[0] - vert0[0] ) + v * (vert2[0] - vert0[0] );
03290                   iP[1] = vert0[1] + u * (vert1[1] - vert0[1] ) + v * (vert2[1] - vert0[1] );
03291                   iP[2] = vert0[2] + u * (vert1[2] - vert0[2] ) + v * (vert2[2] - vert0[2] );
03292                   
03293                   if (d) {
03294                      
03295                      d[0] = (vert0[0] - iP[0])*(vert0[0] - iP[0]) + (vert0[1] - iP[1])*(vert0[1] - iP[1]) + (vert0[2] - iP[2])*(vert0[2] - iP[2]);
03296                      *closest_vert = 0;
03297                      d[1] = (vert1[0] - iP[0])*(vert1[0] - iP[0]) + (vert1[1] - iP[1])*(vert1[1] - iP[1]) + (vert1[2] - iP[2])*(vert1[2] - iP[2]);
03298                      if (d[1] < d[*closest_vert]) {
03299                         *closest_vert = 1;
03300                      }
03301                      d[2] = (vert2[0] - iP[0])*(vert2[0] - iP[0]) + (vert2[1] - iP[1])*(vert2[1] - iP[1]) + (vert2[2] - iP[2])*(vert2[2] - iP[2]);
03302                      if (d[2] < d[*closest_vert]) {
03303                         *closest_vert = 2;
03304                      }
03305                      d[0] = (float)sqrt((double)d[0]);
03306                      d[1] = (float)sqrt((double)d[1]);
03307                      d[2] = (float)sqrt((double)d[2]);
03308                   }
03309                }
03310 
03311             }
03312          }
03313       }
03314    
03315    SUMA_RETURN (hit);
03316 }
03317 
03318 
03319 
03320 
03321 
03322 
03323 
03324 
03325 
03326 
03327 
03328 
03329 
03330 
03331 
03332 
03333 
03334 
03335 
03336 
03337 
03338 
03339 
03340 
03341 
03342 
03343 
03344 
03345 
03346 
03347 
03348 
03349 
03350 
03351  
03352  
03353 SUMA_MT_INTERSECT_TRIANGLE *
03354 SUMA_MT_intersect_triangle(float *P0, float *P1, float *NodeList, int N_Node, int *FaceSetList, int N_FaceSet, SUMA_MT_INTERSECT_TRIANGLE *PrevMTI)
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) { 
03372       entry = 0;
03373       if (LocalHead) fprintf(SUMA_STDERR,"%s: First entry or nothing pre-allocated.\n", FuncName);
03374    } else { 
03375       if (N_FaceSet_Previous != N_FaceSet) { 
03376          if (LocalHead) fprintf(SUMA_STDERR,"%s: Reallocating for MTI, a change in number of FaceSets.\n", FuncName);
03377          
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    
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) {
03426       
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       
03445       SUMA_MT_SUB(edge1, vert1, vert0);
03446       SUMA_MT_SUB(edge2, vert2, vert0);
03447 
03448       
03449       SUMA_MT_CROSS(pvec, dir, edge2);
03450 
03451       
03452       det = SUMA_MT_DOT(edge1, pvec);
03453 
03454    #ifdef SUMA_MT_TEST_CULL           
03455       if (det > -SUMA_EPSILON && det < SUMA_EPSILON)
03456          MTI->isHit[iface] = NOPE;
03457       else {
03458          
03459          SUMA_MT_SUB(tvec, orig, vert0);
03460 
03461          
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             
03467             SUMA_MT_CROSS(qvec, tvec, edge1);
03468 
03469              
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                
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                
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                   
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                   
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                    
03516       if (det > -SUMA_EPSILON && det < SUMA_EPSILON)
03517          MTI->isHit[iface] = NOPE;
03518       else {
03519          inv_det = 1.0 / det;
03520 
03521          
03522          SUMA_MT_SUB(tvec, orig, vert0);
03523 
03524          
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             
03530             SUMA_MT_CROSS(qvec, tvec, edge1);
03531 
03532             
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                
03538                MTI->t[iface] = (float)SUMA_MT_DOT(edge2, qvec) * inv_det;         
03539                MTI->isHit[iface] = YUP;
03540                ++MTI->N_hits;
03541                
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                   
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                   
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    }
03578    MTI->N_el = N_FaceSet;
03579    
03580    ++entry;
03581    N_FaceSet_Previous = N_FaceSet;
03582 
03583    SUMA_RETURN (MTI);
03584 }
03585 
03586 
03587 
03588 
03589 
03590 SUMA_Boolean SUMA_Show_MT_intersect_triangle(SUMA_MT_INTERSECT_TRIANGLE *MTI, FILE *Out)
03591 {
03592    static char FuncName[]={"SUMA_Show_MT_intersect_triangle"};
03593    int MaxShow = 5, i,j;
03594    
03595    SUMA_ENTRY;
03596 
03597    if (Out == NULL) Out = stdout;
03598       
03599    if (MTI == NULL) {
03600       fprintf (Out, "NULL Surface Object Pointer\n");
03601       SUMA_RETURN(NOPE);
03602    }
03603    
03604    fprintf (Out,"\n---------------------------------\n");
03605    if (!MTI->N_el) {
03606       fprintf (Out,"Zero elements in structure\n");
03607       SUMA_RETURN (YUP);
03608    }
03609    
03610    if (MTI->isHit == NULL) {
03611       fprintf (SUMA_STDERR,"Error SUMA_Show_MT_intersect_triangle: isHit is NULL\n\n");
03612       SUMA_RETURN (NOPE);
03613    }
03614    else {
03615       if (MaxShow > MTI->N_el) MaxShow = MTI->N_el; 
03616       fprintf (Out, "Intersection results (showing first %d out of %d elements):\n", MaxShow, MTI->N_el);
03617       for (i=0; i < MaxShow; ++i)   {
03618          fprintf (Out, "\tisHit: %d t %f u %f v %f", MTI->isHit[i], MTI->t[i], MTI->u[i],MTI->v[i]);
03619       }
03620          fprintf (Out, "\n");
03621       
03622       if (MTI->N_hits) {
03623          fprintf (Out, "\n%d hits, (%d hists with positive distance).\n", MTI->N_hits, MTI->N_poshits);
03624          fprintf (Out, "Minimum Distance: %d t %f u %f v %f\n", \
03625                   MTI->ifacemin, MTI->t[MTI->ifacemin], MTI->u[MTI->ifacemin],MTI->v[MTI->ifacemin]);
03626          fprintf (Out, "Intersection point P at Minimum Distance FaceSet:\n%f, %f, %f\n", \
03627                   MTI->P[0], MTI->P[1], MTI->P[2]);
03628          fprintf (Out, "Closest node is number %d in Minimum Distance Faceset (%d in NodeList) at %f distance.\n",\
03629                   MTI->inodeminlocal, MTI->inodemin, MTI->d);                           
03630          fprintf (Out, "Maximum Distance: %d t %f u %f v %f\n\n", \
03631                   MTI->ifacemax, MTI->t[MTI->ifacemax], MTI->u[MTI->ifacemax],MTI->v[MTI->ifacemax]);
03632          fprintf (Out, "Intersection of ray with surface (showing first %d out of %d elements):\n", MaxShow, MTI->N_el);
03633          i = 0;
03634          j = 0;
03635          while (i< MTI->N_el && j < MTI->N_hits) {
03636             if (MTI->isHit[i]) {
03637                ++j;
03638                fprintf (Out, "\tisHit: %d t %f u %f v %f\n", MTI->isHit[i], MTI->t[i], MTI->u[i],MTI->v[i]);
03639             }
03640             ++i;
03641          }
03642          fprintf (Out, "\n");
03643       } else {
03644          fprintf (Out, "No Intersection of ray with surface\n");
03645       }
03646 
03647    }
03648    SUMA_RETURN (YUP);
03649 }
03650 
03651 
03652 
03653 
03654 
03655 
03656 void * SUMA_Free_MT_intersect_triangle(SUMA_MT_INTERSECT_TRIANGLE *MTI)
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 }
03669 
03670 
03671 
03672 
03673 
03674 
03675 
03676 
03677 
03678 
03679 
03680 
03681 
03682 
03683 
03684 
03685 
03686 
03687 SUMA_Boolean SUMA_FromToRotation (float *v0, float *v1, float **mtx)
03688 {
03689    static char FuncName[]={"SUMA_FromToRotation"};
03690    float v[3], vn;
03691    float e, h, f;
03692 
03693    SUMA_ENTRY;
03694 
03695    
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)     
03718    {
03719       float u[3], v[3]; 
03720       float x[3];       
03721       float c1, c2, c3; 
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  
03768    {
03769       #if 0
03770          
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          
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 }
03815 
03816 
03817 
03818 
03819 
03820 
03821 
03822 
03823 
03824 
03825 
03826 
03827 
03828 
03829 SUMA_Boolean   SUMA_mattoquat (float **mat, float *q)
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    
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    }  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    } 
03861    SUMA_RETURN (YUP);
03862 }
03863 
03864 
03865 typedef enum {SUMA_NO_NEIGHB, SUMA_NO_MORE_TO_VISIT, SUMA_VISITED_ALL, SUMA_BAD_SEED} SUMA_TAKE_A_HIKE;
03866 
03867 
03868 
03869 
03870 
03871 
03872 
03873 
03874 
03875 
03876 
03877 
03878 
03879 
03880 
03881 int SUMA_whichTri (SUMA_EDGE_LIST * EL, int n1, int n2, int n3, int IOtrace)
03882 {
03883    static char FuncName[]={"SUMA_whichTri"};
03884    int IncTri_E1[100], IncTri_E2[100], N_IncTri_E1 = 0, N_IncTri_E2 = 0, i, j, Tri= -1;
03885    SUMA_Boolean Found = NOPE;
03886    
03887    if (IOtrace) SUMA_ENTRY;
03888    
03889    Tri = -1;
03890    
03891    if (!SUMA_Get_Incident(n1, n2, EL, IncTri_E1, &N_IncTri_E1, IOtrace)) {
03892       fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_Get_Incident.\n", FuncName);
03893    } else if (!SUMA_Get_Incident(n1, n3, EL, IncTri_E2, &N_IncTri_E2, IOtrace)) {
03894       
03895       fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_Get_Incident.\n", FuncName);
03896    } else if (N_IncTri_E1 > 99 || N_IncTri_E2 > 99 ) {
03897       
03898       fprintf (SUMA_STDERR,"Error %s: Exceeded preallocated space.\n", FuncName);
03899    } else {
03900       
03901       i=0;
03902       Found = NOPE;
03903       while (i < N_IncTri_E1 && !Found) {
03904          j = 0;
03905          while (j < N_IncTri_E2 && !Found) {
03906             if (IncTri_E2[j] == IncTri_E1[i]) { 
03907                Found = YUP;
03908                Tri = IncTri_E2[j];
03909             }
03910             ++j;
03911          }
03912          ++i;
03913       }
03914    }
03915    if (IOtrace) { SUMA_RETURN (Tri); }
03916    else return(Tri);
03917 }
03918 
03919 
03920 
03921 
03922 
03923 
03924 
03925 
03926 
03927 
03928 
03929 int SUMA_isTriLinked (int*T, int *t, int *cn)
03930 {
03931    static char FuncName[]={"SUMA_isTriLinked"};
03932    int ic, in;
03933    
03934    SUMA_ENTRY;
03935 
03936    ic = 0;   
03937    in = 0; 
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; 
03954    }
03955    
03956    SUMA_RETURN (ic);
03957 }
03958 
03959 
03960 
03961 
03962 
03963 
03964 
03965 
03966 
03967 
03968 
03969 int SUMA_isConsistent (int *T, int *t)
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;   
03977    in = 0; 
03978    while (ic < 2 && in < 3) {
03979       if (t[0] == T[in]) {
03980          LOC[ic] = in; 
03981          loc[ic] = 0; 
03982          ++ic;
03983       }else {
03984          if (t[1] == T[in]) {
03985             LOC[ic] = in; 
03986             loc[ic] = 1; 
03987             ++ic;
03988          }else {
03989             if (t[2] == T[in]) {
03990                LOC[ic] = in; 
03991                loc[ic] = 2; 
03992                ++ic;
03993             }
03994          }
03995       }
03996       ++in; 
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    
04006 
04007    
04008    if (d > 1 || d < -1) d = - d /2 ;
04009    if (D > 1 || D < -1) D = - D /2 ;
04010       
04011 
04012    if (d != D) {
04013       
04014       SUMA_RETURN (1);
04015    }
04016    
04017       
04018    
04019    in = t[0];
04020    t[0] = t[2];
04021    t[2] = in;
04022    SUMA_RETURN (-1);
04023 }
04024 
04025 #ifdef SUMA_isConsistent_STANDALONE
04026 
04027 
04028 
04029 
04030 int main (int argc,char *argv[])
04031 {
04032    int Ta[6][3] ={{1, 2, 3}, 
04033                   {2, 3, 1}, 
04034                   {3, 1, 2}, 
04035                   {1, 3, 2},  
04036                   {3, 2, 1}, 
04037                   {2, 1, 3},  };
04038    int ta[6][3] ={{4, 3, 2}, 
04039                   {3, 2, 4}, 
04040                   {2, 4, 3}, 
04041                   {2, 3, 4},  
04042                   {3, 4, 2}, 
04043                   {4, 2, 3},  };   
04044    
04045    int T[3], t[3], i, j;
04046    int tc[3];
04047    
04048    
04049    
04050    printf ("ALL THESE SHOULD BE CONSISTENT:\n");   
04051    for (i=0; i < 3; ++i) {
04052       for (j=0; j < 3; ++ j) {
04053          printf ("\tT=[%d %d %d]\n", Ta[i][0], Ta[i][1], Ta[i][2]);
04054          printf ("\tt=[%d %d %d]\n", ta[j][0], ta[j][1], ta[j][2]);
04055          tc[0]=ta[j][0]; 
04056          tc[1]=ta[j][1];
04057          tc[2]=ta[j][2];
04058          if (SUMA_isConsistent(Ta[i], tc) > 0) {
04059             printf ("\ttriangles consistent.\n");
04060          }else {
04061             printf ("\ttriangles inconsistent.Fixed: %d %d %d\n", tc[0], tc[1], tc[2]);
04062          }
04063       }
04064    }   
04065    for (i=3; i < 6; ++i) {
04066       for (j=3; j < 6; ++ j) {
04067          printf ("\tT=[%d %d %d]\n", Ta[i][0], Ta[i][1], Ta[i][2]);
04068          printf ("\tt=[%d %d %d]\n", ta[j][0], ta[j][1], ta[j][2]);
04069          tc[0]=ta[j][0]; 
04070          tc[1]=ta[j][1];
04071          tc[2]=ta[j][2];
04072          if (SUMA_isConsistent(Ta[i], tc) > 0) {
04073             printf ("\ttriangles consistent.\n");
04074          }else {
04075             printf ("\ttriangles inconsistent.Fixed: %d %d %d\n", tc[0], tc[1], tc[2]);
04076          }
04077       }
04078    }   
04079    
04080    
04081    
04082    printf ("ALL THESE SHOULD BE NOT CONSISTENT:\n");   
04083    for (i=3; i < 6; ++i) {
04084       for (j=0; j < 3; ++ j) {
04085          printf ("\tT=[%d %d %d]\n", Ta[i][0], Ta[i][1], Ta[i][2]);
04086          printf ("\tt=[%d %d %d]\n", ta[j][0], ta[j][1], ta[j][2]);
04087          tc[0]=ta[j][0]; 
04088          tc[1]=ta[j][1];
04089          tc[2]=ta[j][2];
04090          if (SUMA_isConsistent(Ta[i], tc) > 0) {
04091             printf ("\ttriangles consistent.\n");
04092          }else {
04093             printf ("\ttriangles inconsistent.Fixed: %d %d %d\n", tc[0], tc[1], tc[2]);
04094          }
04095       }
04096    }   
04097    for (i=0; i < 3; ++i) {
04098       for (j=3; j < 6; ++ j) {
04099          printf ("\tT=[%d %d %d]\n", Ta[i][0], Ta[i][1], Ta[i][2]);
04100          printf ("\tt=[%d %d %d]\n", ta[j][0], ta[j][1], ta[j][2]);
04101          tc[0]=ta[j][0]; 
04102          tc[1]=ta[j][1];
04103          tc[2]=ta[j][2];
04104          if (SUMA_isConsistent(Ta[i], tc) > 0) {
04105             printf ("\ttriangles consistent.\n");
04106          }else {
04107             printf ("\ttriangles inconsistent.Fixed: %d %d %d\n", tc[0], tc[1], tc[2]);
04108          }
04109       }
04110    }   
04111 
04112                               
04113    printf ("Enter triangle 1:\n");
04114    scanf ("%d %d %d", &T[0], &T[1], &T[2]);
04115    printf ("Enter triangle 2:\n");
04116    scanf ("%d %d %d", &t[0], &t[1], &t[2]);
04117 
04118    if (SUMA_isConsistent(T, t) > 0) {
04119       printf ("triangles consistent.\n");
04120    }else {
04121       printf ("triangles inconsistent.Fixed: %d %d %d\n", t[0], t[1], t[2]);
04122    }
04123    return (0);
04124 
04125 }
04126 
04127 #endif
04128 
04129 
04130 
04131 
04132 
04133 
04134 
04135 
04136 int SUMA_Next_Best_Seed (SUMA_FACESET_FIRST_EDGE_NEIGHB *SFFN, int * visited, int N_FL)
04137 {
04138    static int entry = 0, seed=-1;
04139    int Found1 = -1, Found2 = -1, i, N_NotVisNeighb, itry;
04140    static char FuncName[]={"SUMA_Next_Best_Seed"};
04141    
04142    SUMA_ENTRY;
04143 
04144    if (!entry) { 
04145       for (i=0; i < N_FL; ++i) {
04146          if (SFFN->N_Neighb[i] == 3) {
04147             seed = i; ++entry; SUMA_RETURN(seed);
04148          }
04149          if (SFFN->N_Neighb[i] == 2) Found2 = i;
04150          if (SFFN->N_Neighb[i] == 1) Found1 = i;
04151       }   
04152             
04153       if (Found2 > 0) {
04154          ++entry;
04155          SUMA_RETURN (Found2);
04156       }
04157          
04158       if (Found1 > 0) {
04159          ++entry;
04160          SUMA_RETURN (Found1);
04161       }
04162       
04163       SUMA_RETURN (-1);       
04164    }
04165    else {
04166       for (i=0; i < N_FL; ++i) {
04167          if (visited[i]) { 
04168             
04169             N_NotVisNeighb = 0;
04170             itry = 0;
04171             while (itry < SFFN->N_Neighb[i]) {
04172                if (!visited[SFFN->FirstNeighb[i][itry]]) ++N_NotVisNeighb;
04173                ++itry;
04174             }
04175             if (N_NotVisNeighb == 2) {
04176                seed = i; ++entry; SUMA_RETURN (seed);
04177             }
04178             if (N_NotVisNeighb == 1) {
04179                Found1 = i;
04180             }
04181          } 
04182       }
04183       if (Found1 > 0) {
04184          ++entry;
04185          SUMA_RETURN (Found1);
04186       }
04187       SUMA_RETURN (-1);    
04188    }
04189 }
04190 
04191 
04192 
04193 
04194 
04195 
04196 
04197 
04198 
04199 
04200 
04201 
04202 
04203 SUMA_TAKE_A_HIKE SUMA_Take_A_Hike (SUMA_FACESET_FIRST_EDGE_NEIGHB *SFFN, int *visited, int *N_visited, int *Consistent, int *FL, int N_FL, int seed)
04204 {
04205    static char FuncName[]={"SUMA_Take_A_Hike"};
04206    int NotFound, itry, curface, nxtface, ipcur, ipnxt, NP;
04207    static int entry=0;
04208 
04209    SUMA_ENTRY;
04210    NP = 3;
04211    curface = seed;
04212    if (!visited[curface]) { 
04213       if (!entry) {
04214          *N_visited += 1;
04215          visited[curface] = 1;
04216          Consistent[curface] = 1;
04217          
04218       }else {
04219          fprintf (SUMA_STDERR, "Error %s: You should not send unvisited seeds, except at the very first call.\n", FuncName);
04220          SUMA_RETURN (SUMA_BAD_SEED);
04221       }
04222    }
04223    if (SFFN->N_Neighb[curface] == 0) {
04224       SUMA_RETURN (SUMA_NO_NEIGHB);
04225    }
04226    ++entry;
04227 
04228    while (*N_visited <= N_FL) {
04229       
04230       itry = 0; 
04231       NotFound = 1; 
04232       do {
04233          nxtface = SFFN->FirstNeighb[curface][itry];
04234          if (visited[nxtface]) {
04235             
04236             ++itry;
04237          }else {
04238             if (SFFN->N_Neighb[nxtface] == 1) {
04239                
04240                
04241                visited[nxtface] = 1;
04242                Consistent[nxtface] = SUMA_isConsistent (&(FL[curface*NP]), &(FL[nxtface*NP]));
04243                
04244                *N_visited = *N_visited+1;
04245                ++itry;
04246             } else {
04247                
04248                Consistent[nxtface] = SUMA_isConsistent (&(FL[curface*NP]), &(FL[nxtface*NP]));
04249                visited[nxtface] = 1;
04250                curface = nxtface;
04251                
04252                *N_visited = *N_visited+1;
04253                NotFound = 0;
04254                itry = 100; 
04255             }
04256          }
04257       } while (itry < SFFN->N_Neighb[curface]) ;
04258 
04259       if (NotFound) { 
04260          
04261          SUMA_RETURN (SUMA_NO_MORE_TO_VISIT);
04262       }
04263 
04264    }
04265    
04266    SUMA_RETURN (SUMA_VISITED_ALL);
04267 }
04268 
04269 
04270 
04271 
04272 
04273 
04274 void SUMA_Show_Edge_List (SUMA_EDGE_LIST *EL, FILE *Out)
04275 {
04276    static char FuncName[]={"SUMA_Show_Edge_List"};
04277    int i;
04278    
04279    SUMA_ENTRY;
04280    
04281    if (Out == NULL) Out = stdout;
04282    
04283    fprintf(Out,"\nEL contents:\n");
04284    if (EL->idcode_str) fprintf(Out,"IDcode: %s\n", EL->idcode_str);
04285    else fprintf(Out,"IDcode: NULL\n");
04286    
04287    fprintf(Out,"i-\t[EL[i][0] EL[i][1]]\t[ELps[i][0] ELps[i][1] ELps[i][2] ELps[i][3]]\n");
04288    for (i=0; i < EL->N_EL; ++i) {
04289       fprintf(Out,"%d-\t[%d %d]\t[%d %d %d %d]\n", 
04290                i, EL->EL[i][0], EL->EL[i][1], EL->ELps[i][0], EL->ELps[i][1], EL->ELps[i][2], EL->ELps[i][3]);
04291    
04292    }
04293    fprintf(Out,"\nTriLimb contents:\n");
04294    fprintf(Out,"ti-\t[Edge1 Edge2 Edge3]\n");
04295    for (i=0; i < EL->N_EL/3; ++i) { 
04296       fprintf(Out,"t%d-\t[%d %d %d]\n",
04297          i, EL->Tri_limb[i][0], EL->Tri_limb[i][1],EL->Tri_limb[i][2]);
04298    }
04299    
04300    SUMA_RETURNe;
04301 }
04302 
04303 
04304 
04305 
04306 
04307 
04308 void SUMA_free_Edge_List (SUMA_EDGE_LIST *SEL)
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 }
04327 
04328 
04329 
04330 
04331 SUMA_EDGE_LIST * SUMA_Make_Edge_List (int *FL, int N_FL, int N_Node, float *NodeList, char *ownerid)
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 }
04339 
04340 
04341 
04342 
04343 
04344 
04345 
04346 
04347 
04348 
04349 
04350 
04351 
04352 
04353 
04354 
04355 
04356 
04357 
04358 
04359 
04360 
04361 
04362 
04363 
04364 
04365 
04366 
04367 
04368 
04369 
04370 SUMA_EDGE_LIST * SUMA_Make_Edge_List_eng (int *FL, int N_FL, int N_Node, float *NodeList, int debug, char *ownerid)
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    
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)); 
04402    SEL->ELloc = (int *)SUMA_calloc(N_Node, sizeof(int));
04403    SEL->Le = (float *) SUMA_calloc (SEL->N_EL, sizeof(float)); 
04404    ELp = (int **) SUMA_allocate2D (SEL->N_EL, 2, sizeof(int)); 
04405                                                          
04406     
04407    SEL->ELps = (int **) SUMA_allocate2D (SEL->N_EL, 3, sizeof(int)); 
04408    
04409    
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    
04420    SUMA_LH("Forming Edge List...\n");
04421    for (i=0; i< N_FL; ++i) {
04422       
04423       ie = 3*i;
04424       ip = 3*i;
04425       if (FL[ip] > FL[ip+1]) {
04426          
04427          SEL->EL[ie][0] = FL[ip+1];
04428          SEL->EL[ie][1] = FL[ip];
04429          
04430          ELp[ie][0] = 1; 
04431       } else { 
04432          
04433          SEL->EL[ie][0] = FL[ip];
04434          SEL->EL[ie][1] = FL[ip+1];
04435          ELp[ie][0] = -1; 
04436       }
04437       ELp[ie][1] = i; 
04438 
04439       
04440       ie += 1;
04441       if (FL[ip+1] > FL[ip+2]) {
04442          
04443          SEL->EL[ie][0] = FL[ip+2];
04444          SEL->EL[ie][1] = FL[ip+1];
04445          
04446          ELp[ie][0] = 1; 
04447       } else { 
04448          
04449          SEL->EL[ie][0] = FL[ip+1];
04450          SEL->EL[ie][1] = FL[ip+2];
04451          ELp[ie][0] = -1; 
04452       }
04453       ELp[ie][1] = i; 
04454       
04455       
04456       ie += 1;
04457       if (FL[ip+2] > FL[ip]) {
04458          
04459          SEL->EL[ie][0] = FL[ip];
04460          SEL->EL[ie][1] = FL[ip+2];
04461          
04462          ELp[ie][0] = 1; 
04463       } else { 
04464          
04465          SEL->EL[ie][0] = FL[ip+2];
04466          SEL->EL[ie][1] = FL[ip];
04467          ELp[ie][0] = -1; 
04468       }
04469       ELp[ie][1] = i; 
04470       
04471    }
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    
04484    SUMA_LH("Sorting edge list...");
04485    isort_EL = SUMA_dqsortrow (SEL->EL, SEL->N_EL, 2);
04486    
04487    
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    
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    
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    
04523    SUMA_LH("Searching SEL for funky stuff");
04524    SEL->N_Distinct_Edges = 0;
04525    i=0;
04526    while (i < SEL->N_EL) {
04527       
04528       ht = SEL->ELps[i][1]; 
04529       SEL->Tri_limb[ht][iTri_limb[ht]] = i;
04530       iTri_limb[ht] += 1;
04531       SEL->N_Distinct_Edges += 1;  
04532       SEL->ELps[i][2] = 1; 
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]) {
04536             SEL->ELps[i][2] += 1; 
04537             SEL->ELps[i+lu][2] = -1; 
04538 
04539             
04540             ht = SEL->ELps[i+lu][1]; 
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    
04598 
04599 
04600 
04601 
04602 
04603 
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); 
04619    SUMA_LH("Done with storage, returning...\n");
04620    
04621    SUMA_RETURN (SEL);
04622 }
04623 
04624 
04625 
04626 
04627 
04628 
04629 
04630 
04631 
04632 
04633 
04634 int SUMA_FindEdgeInTri (SUMA_EDGE_LIST *EL, int n1, int n2, int Tri) 
04635 {
04636    static char FuncName[]={"SUMA_FindEdgeInTri"};
04637    int eloc;
04638    
04639    SUMA_ENTRY;
04640 
04641    
04642    if (n2 < n1) {
04643       eloc = n2; 
04644       n2 = n1; 
04645       n1 = eloc;
04646    }
04647    
04648    
04649    eloc = EL->ELloc[n1];
04650    
04651    
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    
04658    SUMA_RETURN (-1);
04659 }
04660 
04661 
04662 
04663 
04664 
04665 
04666 
04667 
04668 
04669 
04670 
04671 int SUMA_FindEdge (SUMA_EDGE_LIST *EL, int n1, int n2) 
04672 {
04673    static char FuncName[]={"SUMA_FindEdge"};
04674    int eloc;
04675    SUMA_Boolean LocalHead = NOPE;
04676    
04677    SUMA_ENTRY;
04678 
04679    
04680    if (n2 < n1) {
04681       eloc = n2; 
04682       n2 = n1; 
04683       n1 = eloc;
04684    }
04685    
04686    
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    
04693    do {
04694       
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    
04700    SUMA_RETURN (-1);
04701 }
04702 
04703 
04704 
04705 
04706 
04707 
04708 
04709 
04710 
04711 
04712 
04713 
04714 
04715 
04716 
04717 SUMA_Boolean SUMA_Get_NodeIncident(int n1, SUMA_SurfaceObject *SO, int *Incident, int *N_Incident)
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       
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 }
04748 
04749 
04750 
04751 
04752 
04753 
04754 
04755 
04756 
04757 
04758 
04759 
04760 
04761 
04762 
04763 SUMA_Boolean SUMA_Get_Incident(int n1, int n2, SUMA_EDGE_LIST *SEL, int *Incident, int *N_Incident, int IOtrace)
04764 {
04765    static char FuncName[] = {"SUMA_Get_Incident"};
04766    int nt, in1, iseek, m_N_EL;
04767    
04768    if (IOtrace) SUMA_ENTRY;
04769 
04770    
04771    if (n1 > n2) {
04772       
04773       nt = n1;
04774       n1 = n2;
04775       n2 = nt;
04776    }
04777    
04778    
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]; 
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    
04798    if (IOtrace) { SUMA_RETURN(YUP); }
04799    else return(YUP);   
04800 }
04801 
04802 
04803 
04804 
04805 
04806 
04807 
04808  
04809 void SUMA_free_FaceSet_Edge_Neighb (SUMA_FACESET_FIRST_EDGE_NEIGHB * S)
04810 {
04811    static char FuncName[]={"SUMA_free_FaceSet_Edge_Neighb"};
04812    
04813    SUMA_ENTRY;
04814 
04815    if (S->FirstNeighb) SUMA_free2D((char **)S->FirstNeighb, S->N_FaceSet);
04816    if (S->N_Neighb) SUMA_free(S->N_Neighb);
04817    if (S) SUMA_free(S);
04818    SUMA_RETURNe;
04819 }
04820 
04821 
04822 
04823 
04824 
04825 
04826 
04827 
04828 SUMA_FACESET_FIRST_EDGE_NEIGHB *SUMA_allocate_FaceSet_Edge_Neighb (int N_FaceSet)
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; 
04849    SFFN->N_FaceSet = N_FaceSet;
04850    SFFN->N_Neighb_min = 100; 
04851    SUMA_RETURN (SFFN);
04852 }
04853 
04854 
04855 
04856 
04857 
04858 
04859 
04860 
04861 
04862 
04863 
04864 
04865 
04866 SUMA_FACESET_FIRST_EDGE_NEIGHB *SUMA_FaceSet_Edge_Neighb (int **EL, int **ELps, int N_EL)
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          
04886          i += 1; 
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 }   
04929 
04930 
04931 
04932 
04933 
04934 
04935 
04936 
04937 
04938 
04939 
04940 
04941 
04942 
04943 
04944 
04945 
04946 
04947 
04948 
04949 SUMA_Boolean SUMA_MakeConsistent (int *FL, int N_FL, SUMA_EDGE_LIST *SEL, int detail, int *trouble) 
04950 {
04951    static char FuncName[]={"SUMA_MakeConsistent"};
04952    
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    
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       
04986       #if 0
04987          
04988          EdgeSeed = 0;
04989          i=EdgeSeed;
04990          ht0 = SEL->ELps[i][1];
04991       #else      
04992          
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          
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             
05009             NotConsistent = SEL->ELps[i][0] * SEL->ELps[i+1][0]; 
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                   
05018                   
05019                   ischecked[ht1] = 1;
05020                   ++N_checked;
05021                } else {
05022                   
05023                   
05024                   ip = NP * ht1;
05025                   it = FL[ip];
05026                   FL[ip] = FL[ip+2];
05027                   FL[ip+2] = it;
05028                   
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             
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                   
05049                   
05050                   ischecked[ht0] = 1;
05051                   ++N_checked;
05052                } else {
05053                   
05054                   
05055                   ip = NP * ht0;
05056                   it = FL[ip];
05057                   FL[ip] = FL[ip+2];
05058                   FL[ip+2] = it;
05059                   
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]) { 
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       
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    
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 }
05115 
05116 #ifdef SUMA_MakeConsistent_STANDALONE
05117 void usage ()
05118    
05119   {
05120           printf ("\nUsage:  SUMA_MakeConsistent <FaceSetList file> <NodeList file>\n");
05121           printf ("To compile: \ngcc -DSUMA_MakeConsistent_STANDALONE -Wall -o SUMA_MakeConsistent SUMA_MiscFunc.c ");
05122           printf ("SUMA_lib.a libmri.a -I/usr/X11R6/include -I./ -L/usr/lib -L/usr/X11R6/lib \n");
05123           printf ("-lm -lGL -lGLU -lGLw -lXmu -lXm -lXt -lXext -lX11 -lMesaGLw -lMesaGLwM \n");
05124           printf ("\t\t\t Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \t Wed Mar 20 14:23:42 EST 2002\n");
05125           exit (0);
05126   }
05127    
05128 int main (int argc,char *argv[])
05129 {
05130    char FuncName[100]; 
05131    float *NodeList;
05132    int *FL, N_FL, i, ip, N_Node, trouble;
05133    SUMA_EDGE_LIST *SEL;
05134    
05135    
05136    sprintf (FuncName,"SUMA_MakeConsistent-Main-");
05137    
05138    
05139    if (argc < 2)
05140        {
05141           usage ();
05142           exit (1);
05143        }
05144    
05145    N_FL = SUMA_float_file_size(argv[1]);
05146    N_Node = SUMA_float_file_size(argv[2]);
05147    
05148    N_FL = N_FL / 3;
05149    FL = (int *)SUMA_calloc(N_FL * 3, sizeof(int));
05150    
05151    N_Node = N_Node / 3;
05152    NodeList = (float *)SUMA_calloc(N_Node *3, sizeof(float));
05153    
05154    SUMA_Read_dfile (argv[1], FL,  N_FL * 3);
05155    SUMA_Read_file (argv[2], NodeList, N_Node *3);
05156    
05157    
05158    SEL = SUMA_Make_Edge_List (FL, N_FL, N_Node, NodeList, 1);
05159    if (SEL == NULL) {
05160       fprintf(SUMA_STDERR, "Error %s: Failed in SUMA_Make_Edge_List.\n", FuncName);
05161       return (NOPE);
05162    }
05163 
05164    if (!SUMA_MakeConsistent (FL, N_FL, SEL, 1, &trouble)) {
05165       fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_MakeConsistent.\n", FuncName);
05166       return (1);
05167    }else {
05168       fprintf(SUMA_STDERR,"%s: Eeeexcellent.\n", FuncName);
05169    }
05170    
05171       fprintf(SUMA_STDERR,"%s REARRANGED TRIANGLES:\n", FuncName);
05172    for (i=0; i<N_FL; ++i) {
05173       ip = 3 * i;
05174       fprintf(SUMA_STDERR," %d %d %d\n", FL[ip], FL[ip+1], FL[ip+2]); 
05175    }
05176    SUMA_free_Edge_List(SEL);
05177 
05178    return (0);
05179 }
05180 
05181 #endif
05182 
05183 
05184 
05185 
05186 
05187 
05188 
05189 
05190 
05191 
05192 
05193 
05194 
05195 
05196 
05197 
05198 
05199 
05200 
05201 
05202 
05203 
05204 
05205 float * SUMA_SmoothAttr_Neighb (float *attr, int N_attr, float *attr_sm, SUMA_NODE_FIRST_NEIGHB *fn, int nr)
05206 {
05207    static char FuncName[]={"SUMA_SmoothAttr_Neighb"};
05208    int ni, im, offs, j;
05209     
05210    SUMA_ENTRY;
05211 
05212    if (attr_sm && attr_sm == attr) {
05213       fprintf (SUMA_STDERR, "Error %s: attr and attr_sm point to the same location. BAD!\n",FuncName);
05214       SUMA_RETURN (NULL); 
05215    }
05216    if (fn == NULL) {
05217       fprintf (SUMA_STDERR, "Error %s: fn is null, nothing to do.\n",FuncName);
05218       SUMA_RETURN (NULL); 
05219    }
05220    if (nr*fn->N_Node != N_attr) {
05221       fprintf (SUMA_STDERR, "Error %s: N_attr (%d) must be equal to nr * fn->N_Node (%d * %d = %d).\n",FuncName, N_attr, nr, fn->N_Node, nr * fn->N_Node);
05222       SUMA_RETURN (NULL); 
05223    }
05224    
05225    attr_sm = (float *)attr_sm;
05226    if (attr_sm == NULL) {
05227       attr_sm = (float *)SUMA_calloc (N_attr, sizeof(float));
05228    }
05229    
05230    if (attr_sm == NULL)
05231    {
05232       fprintf (SUMA_STDERR, "Error %s: Failed to allocate for returning variable.\n", FuncName);
05233       SUMA_RETURN (NULL);
05234    } 
05235    
05236    
05237    for (ni=0; ni < fn->N_Node; ++ni) { 
05238       
05239       if (fn->NodeId[ni] != ni) {
05240          
05241          
05242 
05243          
05244 
05245 
05246          continue;
05247       }
05248       offs = nr * ni;
05249       for (im=0; im<nr; ++im) {
05250          attr_sm[offs+im] = attr[offs+im];
05251          for (j=0; j < fn->N_Neighb[ni]; ++j)
05252          {
05253             attr_sm[offs+im] += attr[nr*fn->FirstNeighb[ni][j]+im]; 
05254          }   
05255          attr_sm[offs+im] /= (fn->N_Neighb[ni]+1);
05256       }
05257    }
05258    
05259    SUMA_RETURN (attr_sm);   
05260 }
05261 
05262 
05263 
05264 
05265 
05266 
05267 
05268 
05269 
05270 
05271 
05272 float * SUMA_SmoothAttr_Neighb_Rec (float *attr, int N_attr, float *attr_sm_orig, 
05273                                     SUMA_NODE_FIRST_NEIGHB *fn, int nr, int N_rep)
05274 {
05275    static char FuncName[]={"SUMA_SmoothAttr_Neighb"};
05276    int i;
05277    float *curr_attr=NULL, *attr_sm=NULL;
05278    SUMA_Boolean LocalHead = NOPE;
05279     
05280    SUMA_ENTRY;
05281 
05282    if (N_rep < 1) {
05283       SUMA_SL_Err("N_rep < 1");
05284       SUMA_RETURN(NULL);
05285    }
05286    
05287    if (N_rep == 1 && attr == attr_sm_orig) {
05288       SUMA_SL_Err("attr = attr_sm_orig && N_rep == 1. BAD.\n");
05289       SUMA_RETURN(NULL);
05290    }
05291    
05292    i = 1;
05293    curr_attr = attr; 
05294    while (i < N_rep) {
05295       
05296       attr_sm = SUMA_SmoothAttr_Neighb (curr_attr, N_attr, NULL, fn, nr);
05297       if (i > 1)  { 
05298          
05299          if (curr_attr) SUMA_free(curr_attr);
05300       }
05301       curr_attr = attr_sm; 
05302       ++i;
05303    }      
05304    
05305    
05306    attr_sm = SUMA_SmoothAttr_Neighb (curr_attr, N_attr, attr_sm_orig, fn, nr);
05307    
05308    
05309    if (i > 1) {
05310       if (curr_attr) SUMA_free(curr_attr);
05311    }
05312       
05313    SUMA_RETURN (attr_sm); 
05314 }
05315  
05316 
05317 
05318 
05319 
05320 
05321 
05322 
05323 
05324 
05325 
05326 
05327 SUMA_NODE_FIRST_NEIGHB * SUMA_Build_FirstNeighb (SUMA_EDGE_LIST *el, int N_Node, char *ownerid)
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    
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    
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          
05384          FN->NodeId[n1] = n1; 
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          
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          
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    }
05411 
05412    
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    
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 
05427         
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                
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    
05469    SUMA_RETURN (FN);
05470 }
05471 
05472 
05473 
05474  
05475 SUMA_Boolean SUMA_Free_FirstNeighb (SUMA_NODE_FIRST_NEIGHB *FN)
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    
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 }
05499 
05500 
05501 
05502 
05503 
05504 
05505 
05506 
05507 
05508 
05509 
05510 SUMA_Boolean SUMA_TriNorm (float *n0, float *n1, float *n2, float *norm)
05511 {
05512    static char FuncName[]={"SUMA_TriNorm"};
05513    int i;
05514    float d1[3], d2[3], d;
05515    
05516    for (i=0; i<3; ++i) {
05517          d1[i] = n0[i] - n1[i];
05518          d2[i] = n1[i] - n2[i];
05519    }
05520    norm[0] = d1[1]*d2[2] - d1[2]*d2[1];
05521    norm[1] = d1[2]*d2[0] - d1[0]*d2[2];
05522    norm[2] = d1[0]*d2[1] - d1[1]*d2[0];  
05523    
05524    d = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
05525    
05526    if (d==0.0) {
05527       norm[0] = norm[1] = norm[2] = 1.0;
05528       SUMA_RETURN (NOPE);
05529    }else {
05530       for (i=0; i<3; ++i) norm[i] /= d;
05531       SUMA_RETURN (YUP);
05532    }
05533 }
05534 
05535 
05536 
05537 
05538 
05539 
05540 
05541 
05542 
05543 
05544 
05545 
05546 
05547 float SUMA_TriSurf3 (float *n0, float *n1, float *n2)
05548 {
05549    static char FuncName[]={"SUMA_TriSurf3"};
05550    float dv[3], dw[3], cross[3], A; 
05551    int i, ii, coord, kk, jj;
05552    
05553    SUMA_ENTRY;
05554    
05555    SUMA_MT_SUB (dv, n1, n0);
05556    SUMA_MT_SUB (dw, n2, n0);
05557    SUMA_MT_CROSS(cross,dv,dw);
05558    SUMA_NORM(A, cross);
05559    A *= 0.5;
05560    
05561    SUMA_RETURN (A); 
05562 }
05563 
05564 
05565 
05566 
05567 
05568 
05569 
05570 
05571 
05572 
05573 
05574 
05575 
05576 float * SUMA_TriSurf3v (float *NodeList, int *FaceSets, int N_FaceSet)
05577 {
05578    static char FuncName[]={"SUMA_TriSurf3v"};
05579    float *A = NULL, *n0, *n1, *n2, a;
05580    int i, i3;
05581 
05582    SUMA_ENTRY;
05583    
05584    A = (float *) SUMA_calloc (N_FaceSet, sizeof(float));
05585    if (A == NULL ) {
05586       fprintf(SUMA_STDERR,"Error %s; Failed to allocate for A \n", FuncName);
05587       SUMA_RETURN (NULL);
05588    }  
05589    
05590    for (i=0;  i<N_FaceSet; ++i) {
05591       i3 = 3*i;
05592       n0 = &(NodeList[3*FaceSets[i3]]);
05593       n1 = &(NodeList[3*FaceSets[i3+1]]);
05594       n2 = &(NodeList[3*FaceSets[i3+2]]);
05595       SUMA_TRI_AREA( n0, n1, n2, A[i]); 
05596       
05597    }
05598    
05599    SUMA_RETURN (A);
05600 }
05601 
05602 
05603 
05604 
05605 
05606 
05607 
05608 
05609 
05610 
05611 
05612 
05613 
05614 
05615 
05616 
05617 
05618 
05619 
05620 float * SUMA_PolySurf3 (float *NodeList, int N_Node, int *FaceSets, int N_FaceSet, int PolyDim, float *FaceNormList, SUMA_Boolean SignedArea)
05621 {
05622    static char FuncName[]={"SUMA_PolySurf3"};
05623    float **V, *A, ax, ay, az, an;
05624    int i, ii, coord, kk, jj, id, ND, ip, NP;
05625    
05626    SUMA_ENTRY;
05627 
05628    ND = 3;
05629    NP = PolyDim;
05630    A = (float *) SUMA_calloc (N_FaceSet, sizeof(float));
05631    V = (float **) SUMA_allocate2D(PolyDim+2, 3, sizeof(float));
05632    
05633    if (A == NULL || V == NULL) {
05634       fprintf(SUMA_STDERR,"Error %s; Failed to allocate for A or V\n", FuncName);
05635       SUMA_RETURN (NULL);
05636    }
05637 
05638    for (i=0; i < N_FaceSet; ++i) {
05639       ip = NP * i;
05640       if (FaceNormList[ip] > 0) ax = FaceNormList[ip];
05641          else ax = -FaceNormList[ip];
05642       
05643       if (FaceNormList[ip+1] > 0) ay = FaceNormList[ip+1];
05644          else ay = -FaceNormList[ip+1];
05645       
05646       if (FaceNormList[ip+2] > 0) az = FaceNormList[ip+2];
05647          else az = -FaceNormList[ip+2];
05648    
05649    
05650       coord = 3;
05651       if (ax > ay) {
05652          if (ax > az) coord = 1;
05653       } else {
05654          if (ay > az) coord = 2;
05655       }
05656    
05657       for (ii=0; ii< PolyDim; ++ii) {
05658          ip = NP * i;
05659          id = ND * FaceSets[ip+ii];
05660          V[ii][0] = NodeList[id];
05661          V[ii][1] = NodeList[id+1];
05662          V[ii][2] = NodeList[id+2];
05663       }
05664       ii = PolyDim;
05665       V[ii][0] = V[0][0]; V[ii][1] = V[0][1]; V[ii][2] = V[0][2];
05666       ii = PolyDim + 1;
05667       V[ii][0] = V[1][0]; V[ii][1] = V[1][1]; V[ii][2] = V[1][2];
05668       
05669       
05670       jj = 2;
05671       kk = 0;
05672       for (ii=1; ii < PolyDim+1; ++ii) {
05673          switch (coord) {
05674             case 1:
05675                A[i] = A[i] + ( V[ii][1] * (V[jj][2] - V[kk][2]) );
05676                break;
05677             case 2:
05678                A[i] = A[i] + ( V[ii][0] * (V[jj][2] - V[kk][2]) );
05679                break;
05680             case 3:
05681                A[i] = A[i] + ( V[ii][0] * (V[jj][1] - V[kk][1]) );
05682                break;
05683          }
05684          
05685          ++jj;
05686          ++kk;
05687          
05688       }
05689       
05690       
05691       an = (float) sqrt(ax * ax + ay * ay + az * az);
05692       switch (coord) {
05693          case 1:
05694             A[i] = (A[i] * (an / (2*ax)));
05695             break;
05696          case 2:
05697             A[i] = (A[i] * (an / (2*ay)));
05698             break;
05699          case 3:
05700             A[i] = (A[i] * (an / (2*az)));
05701             break;
05702       }
05703       
05704       if (!SignedArea) {
05705          if (A[i] < 0) A[i] = -A[i];
05706       }
05707    } 
05708    
05709    SUMA_free2D((char **)V, PolyDim+2);
05710    SUMA_RETURN (A);
05711 }
05712 
05713 
05714 #define MAX_INCIDENT_TRI 200
05715 #define DBG_1 
05716 
05717 #ifdef DBG_3
05718    #define DBG_2
05719 #endif
05720 
05721 #ifdef DBG_2
05722    #define DBG_1
05723 #endif
05724 
05725 
05726 
05727 
05728 
05729 
05730 
05731 
05732 
05733 
05734 
05735 
05736 
05737 
05738 
05739 
05740 
05741 
05742 
05743 
05744 
05745 
05746 SUMA_SURFACE_CURVATURE * SUMA_Surface_Curvature (float *NodeList, int N_Node, float *NodeNormList, float *A, int N_FaceSet, SUMA_NODE_FIRST_NEIGHB *FN, SUMA_EDGE_LIST *SEL)
05747 { 
05748    static char FuncName[] = {"SUMA_Surface_Curvature"};
05749    int i, N_Neighb, j, ji, Incident[MAX_INCIDENT_TRI], N_Incident, kk, ii, id, ND; 
05750    float  Ntmp[3],  vi[3], vj[3], *Num, NumNorm, num, denum, sWij, T1e[3], T2e[3], mg, c, s;
05751    float **fa33, **fb33, **fc33, **Ni, **Nit, *Wij, *Kij, **Tij, **I, **Mi, **Q, **Qt, **fa22, **mMi, **mMir;
05752    SUMA_Boolean *SkipNode;
05753    SUMA_SURFACE_CURVATURE *SC;
05754    
05755    SUMA_ENTRY;
05756    
05757    if (!A || !NodeList || !NodeNormList || !FN || !SEL) {
05758       fprintf (SUMA_STDERR, "Error %s: One of your inputs is NULL.\n", FuncName);
05759       SUMA_RETURN(NULL);
05760    }
05761    
05762    SC = (SUMA_SURFACE_CURVATURE *)SUMA_malloc (sizeof(SUMA_SURFACE_CURVATURE));
05763    if (!SC) {
05764       fprintf (SUMA_STDERR, "Error %s: Failed to allocate for SC.\n", FuncName);
05765       SUMA_RETURN(NULL);
05766    }
05767    
05768    Wij = (float *)SUMA_calloc (FN->N_Neighb_max, sizeof(float));
05769    Kij = (float *)SUMA_calloc (FN->N_Neighb_max, sizeof(float));
05770    Num = (float *)SUMA_calloc (3, sizeof(float));
05771    SkipNode = (SUMA_Boolean *) SUMA_calloc (N_Node, sizeof(SUMA_Boolean));
05772    mMi = (float **) SUMA_allocate2D (2,2, sizeof(float));
05773    mMir =(float **) SUMA_allocate2D (2,2, sizeof(float));
05774    fa22 =(float **) SUMA_allocate2D (2,2, sizeof(float));
05775    Tij = (float **) SUMA_allocate2D (FN->N_Neighb_max, 3, sizeof(float));
05776    Ni =  (float **) SUMA_allocate2D (3, 1, sizeof(float));
05777    Nit = (float **) SUMA_allocate2D (1, 3, sizeof(float));
05778    fa33 =(float **) SUMA_allocate2D (3, 3, sizeof(float));
05779    fb33 =(float **) SUMA_allocate2D (3, 3, sizeof(float));
05780    fc33 =(float **) SUMA_allocate2D (3, 3, sizeof(float));
05781    I =   (float **) SUMA_allocate2D (3, 3, sizeof(float));
05782    Q =   (float **) SUMA_allocate2D (3, 3, sizeof(float));
05783    Qt =  (float **) SUMA_allocate2D (3, 3, sizeof(float));
05784    Mi =  (float **) SUMA_allocate2D (3, 3, sizeof(float));
05785    SC->T1 = (float **) SUMA_allocate2D (N_Node, 3, sizeof(float));
05786    SC->T2 = (float **) SUMA_allocate2D (N_Node, 3, sizeof(float));
05787    SC->Kp1 =(float *)SUMA_calloc (N_Node, sizeof(float));
05788    SC->Kp2 =(float *)SUMA_calloc (N_Node, sizeof(float));
05789 
05790    if (!fa22 || !mMir || !mMi || !Wij || !Kij || !Tij || !Ni || !Nit || !fa33 || !fb33 || !fc33 || !I || !Num || !SkipNode || !Mi || !Q || !Qt || !SC->T1 || !SC->T2 || !SC->Kp1 || !SC->Kp2) {
05791       fprintf (SUMA_STDERR, "Error %s: Failed to allocate for Wij, Kij, Tij.\n", FuncName);
05792       if (Wij) SUMA_free(Wij);
05793       if (Kij) SUMA_free(Kij);
05794       if (Num) SUMA_free(Num);
05795       if (SkipNode) SUMA_free(SkipNode);
05796       if (mMi) SUMA_free2D((char **)mMi, 2);
05797       if (mMir) SUMA_free2D((char **)mMir, 2);
05798       if (fa22) SUMA_free2D((char **)fa22, 2);
05799       if (Tij) SUMA_free2D((char **)Tij, FN->N_Neighb_max);
05800       if (Ni) SUMA_free2D((char **)Ni, 3);
05801       if (Nit) SUMA_free2D((char **)Nit, 1);
05802       if (fa33) SUMA_free2D((char **)fa33, 3);
05803       if (fb33) SUMA_free2D((char **)fb33, 3);
05804       if (I) SUMA_free2D((char **)I, 3);
05805       if (Q) SUMA_free2D((char **)Q, 3);
05806       if (Qt) SUMA_free2D((char **)Qt, 3);
05807       if (Mi) SUMA_free2D((char **)Mi, 3);
05808       if (SC) SUMA_Free_SURFACE_CURVATURE (SC);      
05809       SUMA_RETURN(NULL);
05810    }
05811 
05812    
05813    I[0][0] = I[1][1] = I[2][2] = 1.0; I[0][1] = I[0][2] = I[1][0] = I[1][2] = I[2][0] = I[2][1] = 0.0;
05814    
05815    
05816    SC->N_SkipNode = 0;
05817    SC->N_Node = N_Node;
05818    
05819    fprintf (SUMA_STDERR, "%s: Beginning curvature computations:\n", FuncName);
05820    
05821    ND = 3;
05822    SC->N_SkipNode = 0;
05823    for (i=0; i < N_Node; ++i) { 
05824       #ifdef DBG_1
05825          if (!(i%10000)) {
05826             fprintf (SUMA_STDERR, "%s: [%d]/[%d] %.2f/100%% completed\n", FuncName, i, N_Node, (float)i / N_Node * 100);
05827          }
05828       #endif
05829       SkipNode[i] = NOPE;
05830       
05831       N_Neighb = FN->N_Neighb[i];
05832       id = ND * i;
05833       Ni[0][0] = NodeNormList[id]; Ni[1][0] = NodeNormList[id+1]; Ni[2][0] = NodeNormList[id+2]; 
05834       Nit[0][0] = NodeNormList[id]; Nit[0][1] = NodeNormList[id+1]; Nit[0][2] = NodeNormList[id+2];  
05835       vi[0] = NodeList[id]; vi[1] = NodeList[id+1]; vi[2] = NodeList[id+2];   
05836       #ifdef DBG_2
05837          fprintf (SUMA_STDERR, "%s: Looping over neighbors, i = %d\n", FuncName, i);
05838       #endif
05839       j=0;
05840       sWij = 0.0;
05841       while (j < N_Neighb) {
05842          ji = FN->FirstNeighb[i][j]; 
05843          id = ND * ji;
05844          vj[0] = NodeList[id]; vj[1] = NodeList[id+1]; vj[2] = NodeList[id+2];  
05845          
05846          
05847          #ifdef DBG_2
05848             fprintf (SUMA_STDERR, "%s: Mat Op j=%d\n", FuncName, j);
05849          #endif
05850          
05851          
05852          SUMA_MULT_MAT(Ni,Nit,fa33,3,1,3,float,float,float); 
05853          
05854          
05855          SUMA_SUB_MAT(I, fa33, fb33, 3, 3, float, float, float); 
05856 
05857          
05858          fa33[0][0] = vi[0] - vj[0];  fa33[1][0] = vi[1] - vj[1]; fa33[2][0] = vi[2] - vj[2];
05859          
05860          
05861          SUMA_MULT_MAT(fb33, fa33, fc33, 3, 3, 1, float, float, float);
05862          Num[0] = fc33[0][0]; Num[1] = fc33[1][0]; Num[2] = fc33[2][0];
05863 
05864          
05865          NumNorm = (float)sqrt(Num[0]*Num[0] + Num[1]*Num[1] + Num[2]*Num[2]);
05866          if (NumNorm == 0) {
05867             fprintf (SUMA_STDERR, "Warning %s: NumNorm = 0 for node %d.\n", FuncName, i); 
05868             SkipNode[i] = YUP;
05869             SC->N_SkipNode++;
05870             break;
05871          }
05872          
05873          Tij[j][0] = Num[0] / NumNorm; 
05874          Tij[j][1] = Num[1] / NumNorm; 
05875          Tij[j][2] = Num[2] / NumNorm;
05876          
05877          #ifdef DBG_2
05878             fprintf(SUMA_STDOUT,"%s: i,j, ji =%d,%d, %d Ni = %f %f %f\n Tij(%d,:) = %f %f %f.\n", \
05879                FuncName, i, j, ji,Ni[0][0], Ni[1][0], Ni[2][0], j, Tij[j][0], Tij[j][1], Tij[j][2]);
05880          #endif
05881           
05882          
05883          
05884          fa33[0][0] = (vj[0] - vi[0]);  fa33[1][0] = (vj[1] - vi[1]); fa33[2][0] = (vj[2] - vi[2]);
05885          
05886          SUMA_MULT_MAT(Nit, fa33, fb33, 1, 3, 1, float, float, float);
05887          num = fb33[0][0]; 
05888          
05889          denum = fa33[0][0] * fa33[0][0] + fa33[1][0] * fa33[1][0]+ fa33[2][0] * fa33[2][0]; 
05890          
05891          Kij[j] = 2 * num / denum;
05892          #ifdef DBG_2
05893             fprintf(SUMA_STDOUT,"%s: Kij[%d] = %f\n", FuncName, j, Kij[j]);
05894          #endif
05895          
05896          
05897             
05898             if (!SUMA_Get_Incident(i, ji, SEL, Incident, &N_Incident, 1))
05899             {
05900                fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_Get_Incident.\n", FuncName);
05901                if (Wij) SUMA_free(Wij);
05902                if (Kij) SUMA_free(Kij);
05903                if (Num) SUMA_free(Num);
05904                if (SkipNode) SUMA_free(SkipNode);
05905                if (mMi) SUMA_free2D((char **)mMi, 2);
05906                if (mMir) SUMA_free2D((char **)mMir, 2);
05907                if (fa22) SUMA_free2D((char **)fa22, 2);
05908                if (Tij) SUMA_free2D((char **)Tij, FN->N_Neighb_max);
05909                if (Ni) SUMA_free2D((char **)Ni, 3);
05910                if (Nit) SUMA_free2D((char **)Nit, 1);
05911                if (fa33) SUMA_free2D((char **)fa33, 3);
05912                if (fb33) SUMA_free2D((char **)fb33, 3);
05913                if (I) SUMA_free2D((char **)I, 3);
05914                if (Q) SUMA_free2D((char **)Q, 3);
05915                if (Qt) SUMA_free2D((char **)Qt, 3);
05916                if (Mi) SUMA_free2D((char **)Mi, 3);
05917                if (SC) SUMA_Free_SURFACE_CURVATURE (SC);      
05918                SUMA_RETURN(NULL);
05919             }
05920 
05921             #ifdef DBG_2
05922                fprintf (SUMA_STDERR,"%s: Incidents ...\n", FuncName);
05923                for (kk=0; kk < N_Incident; ++kk) {
05924                   fprintf (SUMA_STDERR,"\t %d", Incident[kk]);
05925                }
05926                fprintf (SUMA_STDERR,"\n");
05927             #endif
05928 
05929             if (N_Incident != 2 && N_Incident != 1)
05930             {
05931                fprintf (SUMA_STDERR,"Warning %s: Unexpected N_Incident = %d at i,j = %d,%d\n", FuncName, N_Incident, i, j);
05932                SkipNode[i] = YUP;
05933                ++SC->N_SkipNode;
05934                break;
05935             }
05936             Wij[j] = 0.0; 
05937             for (ii=0; ii < N_Incident; ++ii) {
05938                Wij[j] = Wij[j] + fabs(A[Incident[ii]]);
05939             }
05940             sWij += Wij[j];
05941             if (Wij[j] == 0.0) {
05942                fprintf (SUMA_STDERR,"Warning %s: Null Wij[%d] at i,j=%d,%d\n", FuncName, j, i, j);
05943                SkipNode[i] = YUP;
05944                ++SC->N_SkipNode;
05945                break; 
05946             }
05947          
05948          ++j;
05949          
05950       }
05951       if (!SkipNode[i]) {
05952             
05953             #ifdef DBG_2   
05954                fprintf (SUMA_STDERR,"%s: Wij:\n", FuncName);
05955             #endif
05956             for (ii=0; ii < N_Neighb; ++ii) {
05957                Wij[ii] /= sWij; 
05958                
05959             }
05960             #ifdef DBG_2   
05961                fprintf (SUMA_STDERR,"\n");
05962             #endif
05963             
05964             Mi[0][0] = Mi[1][0] = Mi[2][0] = Mi[0][1] = Mi[1][1] = Mi[2][1] = Mi[0][2] = Mi[1][2] = Mi[2][2] = 0.0;
05965             for (j=0; j < N_Neighb; ++j) {
05966                
05967                 fa33[0][0] = Tij[j][0]; fa33[1][0] = Tij[j][1]; fa33[2][0] = Tij[j][2];
05968                fb33[0][0] = Tij[j][0]; fb33[0][1] = Tij[j][1]; fb33[0][2] = Tij[j][2];
05969                SUMA_MULT_MAT (fa33, fb33, fc33, 3, 1, 3, float, float, float);
05970                
05971                for (ii=0; ii < 3; ++ii) {
05972                   for (kk=0; kk < 3; ++kk) {
05973                      Mi[ii][kk] = Mi[ii][kk] + Wij[j] * Kij[j] * fc33[ii][kk];
05974                   }
05975                }
05976             }
05977             #ifdef DBG_2   
05978                SUMA_disp_mat (Mi, 3, 3, 1);
05979             #endif
05980             
05981             Ntmp[0] = Ni[0][0]; Ntmp[1] = Ni[1][0]; Ntmp[2] = Ni[2][0];
05982             if (!SUMA_Householder(Ntmp, Q)) {
05983                fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_Householder for node %d.\n ", FuncName, i);
05984                mg = 0.0;
05985                SkipNode[i] = YUP;
05986                SC->N_SkipNode++;
05987             } else {
05988                T1e[0] = Q[0][1]; T1e[1] = Q[1][1]; T1e[2] = Q[2][1];
05989                T2e[0] = Q[0][2]; T2e[1] = Q[1][2]; T2e[2] = Q[2][2];  
05990                SUMA_TRANSP_MAT (Q, Qt, 3, 3, float, float);
05991                #ifdef DBG_2   
05992                   SUMA_disp_mat (Q, 3, 3, 1);
05993                   SUMA_disp_mat (Qt, 3, 3, 1);
05994                #endif
05995                
05996                
05997                SUMA_MULT_MAT (Qt, Mi, fa33, 3, 3, 3, float, float, float);
05998                SUMA_MULT_MAT (fa33, Q, Mi, 3, 3, 3, float, float, float);
05999                #ifdef DBG_2   
06000                   SUMA_disp_mat (Mi, 3, 3, 1);
06001                #endif
06002                mMi[0][0] = Mi[1][1]; mMi[0][1] = Mi[1][2];
06003                mMi[1][0] = Mi[2][1]; mMi[1][1] = Mi[2][2]; 
06004                
06005                
06006                mg = sqrt(mMi[0][0]*mMi[0][0] + mMi[1][0]*mMi[1][0]);
06007                c = mMi[0][0] / mg;
06008                s = mMi[1][0] / mg;
06009                
06010                fa22[0][0] =  c; fa22[0][1] = s; 
06011                fa22[1][0] = -s; fa22[1][1] = c; 
06012                SUMA_MULT_MAT(fa22, mMi, mMir, 2, 2, 2, float, float, float);
06013                #ifdef DBG_2   
06014                   fprintf (SUMA_STDERR,"%s: mg = %f, c = %f, s = %f\n", FuncName,  mg, c, s);
06015                #endif   
06016                
06017                SC->T1[i][0] = c * T1e[0] - s * T2e[0];               
06018                SC->T1[i][1] = c * T1e[1] - s * T2e[1];               
06019                SC->T1[i][2] = c * T1e[2] - s * T2e[2];   
06020                            
06021                SC->T2[i][0] = s * T1e[0] + c * T2e[0];               
06022                SC->T2[i][1] = s * T1e[1] + c * T2e[1];               
06023                SC->T2[i][2] = s * T1e[2] + c * T2e[2];   
06024                
06025                
06026                SC->Kp1[i] = 3 * mMir[0][0] - mMir[1][1];
06027                SC->Kp2[i] = 3 * mMir[1][1] - mMir[0][0];
06028                #ifdef DBG_2   
06029                   fprintf (SUMA_STDERR,"%s: SC->Kp1[i] = %f, SC->Kp2[i] = %f, mKp[i] = %f\n", FuncName,  SC->Kp1[i], SC->Kp2[i], (SC->Kp1[i]+SC->Kp2[i])/2);
06030                #endif   
06031             }
06032             
06033          #ifdef DBG_3
06034             SUMA_PAUSE_PROMPT("Done with node, waiting to move to next");
06035          #endif
06036          
06037       } 
06038       if (SkipNode[i]) {
06039          SC->T1[i][0] = SC->T1[i][1] = SC->T1[i][2] = 0.0;
06040          SC->T2[i][0] = SC->T2[i][1] = SC->T2[i][2] = 0.0;
06041          SC->Kp1[i] = SC->Kp2[i] = 0.0;
06042       }
06043    }
06044 
06045    
06046    {
06047       FILE *fid;
06048       fprintf(SUMA_STDOUT,"%s: Writing Kp1 & Kp2 to Curvs_c.txt ...", FuncName);
06049       fid = fopen("Curvs_c.txt","w");
06050       for (ii=0; ii < SC->N_Node; ++ii) {
06051          
06052          fprintf(fid,"%f %f\n", SC->Kp1[ii], SC->Kp2[ii]);
06053       }
06054       fclose (fid);
06055       
06056       fprintf(SUMA_STDOUT,"%s: Done.\n", FuncName);
06057    }
06058    
06059    
06060    if (Wij) SUMA_free(Wij);
06061    if (Kij) SUMA_free(Kij);
06062    if (Num) SUMA_free(Num);
06063    if (SkipNode) SUMA_free(SkipNode);
06064    if (mMi) SUMA_free2D((char **)mMi, 2);
06065    if (mMir) SUMA_free2D((char **)mMir, 2);
06066    if (fa22) SUMA_free2D((char **)fa22, 2);
06067    if (Tij) SUMA_free2D((char **)Tij, FN->N_Neighb_max);
06068    if (Ni) SUMA_free2D((char **)Ni, 3);
06069    if (Nit) SUMA_free2D((char **)Nit, 1);
06070    if (fa33) SUMA_free2D((char **)fa33, 3);
06071    if (fb33) SUMA_free2D((char **)fb33, 3);
06072    if (I) SUMA_free2D((char **)I, 3);
06073    if (Q) SUMA_free2D((char **)Q, 3);
06074    if (Qt) SUMA_free2D((char **)Qt, 3);
06075    if (Mi) SUMA_free2D((char **)Mi, 3);
06076    
06077    fprintf (SUMA_STDERR, "%s: Done with curvature computations.\n", FuncName);
06078 
06079    SUMA_RETURN (SC);
06080 }
06081 
06082 
06083 
06084 
06085 void SUMA_Free_SURFACE_CURVATURE (SUMA_SURFACE_CURVATURE *SC)
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 }
06099 
06100 
06101 
06102 
06103 
06104 
06105 
06106 
06107 
06108 
06109 
06110 
06111 
06112 
06113 
06114 
06115 #define TAUBIN_Householder
06116 SUMA_Boolean SUMA_Householder (float *Ni, float **Q)
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       
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        
06144       
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          
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          
06170       }
06171       
06172    #endif 
06173    
06174    
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 }
06189 
06190 
06191 
06192 
06193 
06194 
06195 
06196 
06197 
06198 
06199 
06200 
06201 
06202 
06203 
06204 
06205 
06206 
06207 
06208 
06209 
06210 
06211 
06212 
06213 
06214 
06215 
06216 
06217 float * SUMA_Convexity (float *NL, int N_N, float *NNL, SUMA_NODE_FIRST_NEIGHB *FN) 
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 }
06229 
06230 
06231 
06232 
06233 
06234 
06235 
06236 
06237 
06238 
06239 
06240 
06241 
06242 
06243 
06244 
06245 
06246 
06247 
06248 
06249 
06250 
06251              
06252 float * SUMA_Convexity_Engine (float *NL, int N_N, float *NNL, SUMA_NODE_FIRST_NEIGHB *FN, char *DetailFile)
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    
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       
06281 
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          
06289 
06290 
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          
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          
06300          
06301          
06302          
06303 
06304           
06305          C[i] -= d/dij; 
06306          
06307          if (DetailFile) fprintf(fid,"%f\t%f\t%f\t", d, dij, d/dij);
06308          
06309       }
06310       
06311       if (DetailFile) {
06312          
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);  
06320    
06321    #if 0
06322    {
06323       
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 } 
06337 
06338 
06339 
06340 
06341 
06342 
06343 
06344 
06345 
06346 
06347 
06348 
06349 
06350 
06351 
06352 
06353 char * SUMA_pad_str ( char *str, char pad_val , int pad_ln , int opt)
06354 {
06355     static char FuncName[]={"SUMA_pad_str"};
06356    int lo,i;
06357     char *strp , *buf1;
06358 
06359    SUMA_ENTRY;
06360 
06361     assert (str);
06362 
06363     lo = (int)strlen(str);
06364 
06365    buf1 = (char *)SUMA_calloc (pad_ln-lo+2,sizeof (char));
06366    strp = (char *)SUMA_calloc (pad_ln+lo+2,sizeof (char));
06367 
06368    for (i=0;i<pad_ln-lo;++i)
06369        {
06370           if (i == 0) sprintf (buf1,"%c",pad_val);
06371              else sprintf (buf1,"%s%c",buf1,pad_val);
06372 
06373        }
06374     if (opt == 0)
06375        sprintf (strp,"%s%s",buf1,str);
06376     else if (opt == 1)
06377        {
06378           sprintf (strp,"%s%s",str,buf1);
06379 
06380       }         
06381        else 
06382           {
06383              fprintf (SUMA_STDERR, "Error %s: Wrong opt paramter, only (0,1) allowed\n", FuncName);
06384              SUMA_free(strp);
06385             SUMA_free(buf1);
06386             SUMA_RETURN (NULL);
06387           }
06388 
06389     SUMA_free(buf1);
06390 
06391     SUMA_RETURN (strp);
06392 
06393 }
06394 
06395 
06396 char SUMA_ReadCharStdin (char def, int case_sensitive, char *allowed)
06397 {
06398    static char FuncName[]={"SUMA_ReadCharStdin"};
06399    char str[10], *strback;
06400    char cbuf;
06401    int Done, i, nc;
06402    
06403    SUMA_ENTRY;
06404    
06405    do {
06406       Done = 1;
06407        
06408       str[0] = def;
06409       strback = fgets(str, 2*sizeof(char), stdin);
06410       cbuf = str[0];
06411       if (SUMA_IS_BLANK(str[0])) {
06412          cbuf = def;
06413       }
06414 
06415       if (!case_sensitive) {
06416          if (cbuf >= 'A' && cbuf <= 'Z') cbuf = cbuf + 'a' - 'A';  
06417       }
06418       
06419       if (allowed && cbuf) {
06420          
06421          nc = strlen(allowed);
06422          for (i=0; i<nc;++i) {
06423             if (cbuf == allowed[i]) SUMA_RETURN(cbuf); 
06424          }
06425          Done = 0;
06426          
06427          fprintf(stdout,"\abad input, try again: "); fflush(stdout);
06428       }
06429    } while (!Done);
06430    SUMA_RETURN(cbuf);
06431 }
06432 
06433 
06434 
06435 
06436 
06437 
06438 
06439 
06440 
06441 
06442 
06443 
06444 
06445 
06446 int SUMA_ReadNumStdin (float *fv, int nv)
06447 {   
06448    int i=0, nvr = 0;
06449    char *endp, *strtp, s[SUMA_MAX_STRING_LENGTH], cbuf;
06450    static char FuncName[]={"SUMA_ReadNumStdin"};
06451    SUMA_Boolean eos, LocalHead = NOPE;
06452    
06453    SUMA_ENTRY;
06454 
06455    fflush (stdin);
06456    
06457    while ((cbuf = getc(stdin)) != '\n' && i < SUMA_MAX_STRING_LENGTH-1) {
06458       if (cbuf == ',' || cbuf == '\t') {
06459          cbuf = ' ';
06460       }
06461          s[i] = cbuf;
06462          ++ i;
06463    }
06464    
06465    if (i == SUMA_MAX_STRING_LENGTH-1) {
06466       fprintf(SUMA_STDERR,"Error %s: No more than %d characters are allowed on stdin.\n", FuncName, SUMA_MAX_STRING_LENGTH-1);
06467       fflush(stdin);
06468       SUMA_RETURN(-1);
06469    }
06470    
06471    s[i] = '\0';
06472    
06473    if (!i) SUMA_RETURN(0);
06474    
06475    
06476    strtp = s;
06477    endp = NULL;
06478    nvr = 0;
06479    eos = NOPE;
06480    while (nvr < nv && !eos) {
06481       fv[nvr] = strtod(strtp, &endp);
06482       if (LocalHead) fprintf (SUMA_STDERR, "Local Debug %s: ERANGE: %d, EDOM %d, errno %d\n", FuncName, ERANGE, EDOM, errno); 
06483       
06484       if (endp == strtp) { 
06485          eos = YUP;
06486       } else {
06487          ++nvr;
06488          strtp = endp;
06489       }
06490    }
06491    
06492    if (eos && nvr < nv) {
06493       fprintf (SUMA_STDERR, "Warning %s: Expected to read %d elements, read only %d.\n", FuncName, nv, nvr);
06494    }
06495    
06496    SUMA_RETURN(nvr);
06497 }
06498 
06499         
06500 
06501 
06502 
06503 
06504 
06505 
06506 
06507 
06508 
06509 
06510 
06511 
06512 
06513 
06514 
06515 
06516 
06517 
06518 
06519 
06520 
06521 
06522 
06523 
06524 
06525 
06526 
06527 
06528 
06529 
06530 
06531 
06532 
06533 
06534 
06535 
06536 
06537 
06538 
06539 int * SUMA_Find_inIntVect (int *x, int xsz, int val, int *nValLocation)
06540 {
06541    int k, *tmp, *ValLocation;
06542    static char FuncName[]={"SUMA_Find_inIntVect"};
06543    SUMA_Boolean LocalHead = NOPE;
06544    
06545    SUMA_ENTRY;
06546 
06547    
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    
06568       ValLocation = (int *) SUMA_calloc(*nValLocation,sizeof(int));
06569    
06570       SUMA_SCALE_VEC(tmp,ValLocation,1,*nValLocation,int,int);
06571    
06572       SUMA_free(tmp);
06573 
06574    SUMA_RETURN (ValLocation);
06575 
06576 }
06577 
06578 
06579 
06580 
06581 
06582 
06583 
06584 
06585 
06586 
06587 
06588 
06589 
06590 
06591 
06592 
06593 
06594 
06595 
06596 
06597 
06598 
06599 
06600 
06601 
06602 
06603 
06604 
06605 
06606 
06607 
06608 
06609 
06610 
06611 
06612 
06613 
06614 
06615 int * SUMA_UniqueInt (int *y, int xsz, int *kunq, int Sorted )
06616 {
06617    int *xtmp, *xunq, k ,*x;
06618    SUMA_Boolean LocalHead = NOPE;
06619    static char FuncName[]={"SUMA_UniqueInt"};
06620 
06621    SUMA_ENTRY;
06622    *kunq = 0;
06623 
06624    if (!xsz)
06625     {
06626       SUMA_RETURN(NULL);
06627    }
06628    if (!Sorted)
06629     {
06630       x = (int *)SUMA_calloc(xsz, sizeof(int));
06631       if (!x)
06632          {
06633             fprintf (SUMA_STDERR,"Error %s: Failed to allocate for x.", FuncName);
06634             SUMA_RETURN (NULL);
06635          }
06636       for (k=0; k < xsz; ++k)
06637          x[k] = y[k];
06638       qsort(x,xsz,sizeof(int), (int(*) (const void *, const void *)) SUMA_compare_int);
06639    }
06640    else
06641       x = y;
06642 
06643    if (!xsz)   
06644     SUMA_RETURN (NULL);
06645 
06646    xtmp = (int *) SUMA_calloc(xsz,sizeof(int));
06647    if (xtmp == NULL)
06648     {
06649       fprintf (SUMA_STDERR,"Error %s: Could not allocate memory", FuncName);
06650       SUMA_RETURN (NULL);
06651    }
06652 
06653    *kunq = 0;
06654    xtmp[0] = x[0];
06655    for (k=1;k<xsz;++k)
06656     {
06657       if ((x[k] != x[k - 1]))
06658          {
06659             ++*kunq;
06660             xtmp[*kunq] = x[k];   
06661          }
06662    }
06663    ++*kunq;
06664    
06665    
06666    
06667    xunq = (int *) SUMA_calloc(*kunq,sizeof(int));
06668    SUMA_COPY_VEC(xtmp,xunq,*kunq,int,int);
06669 
06670    SUMA_free(xtmp); 
06671 
06672    if (!Sorted)
06673       SUMA_free (x);
06674 
06675    SUMA_RETURN (xunq);
06676 }
06677 
06678 
06679 
06680 
06681 
06682 
06683 
06684 
06685 
06686 
06687 
06688 
06689 
06690 
06691 
06692 
06693 
06694 
06695 
06696 
06697 
06698 int * SUMA_UniqueInt_ind (int *ys, int N_y, int *kunq, int **iup)
06699 {
06700    int *yu=NULL, k ,*iu=NULL;
06701    SUMA_Boolean LocalHead = NOPE;
06702    static char FuncName[]={"SUMA_UniqueInt_ind"};
06703 
06704    SUMA_ENTRY;
06705    
06706    *kunq = 0;
06707 
06708    if (!N_y)
06709     {
06710       SUMA_RETURN(NULL);
06711    }
06712 
06713    if (!N_y)   
06714     SUMA_RETURN (NULL);
06715 
06716    yu = (int *) SUMA_calloc(N_y,sizeof(int));
06717    iu = (int *) SUMA_calloc(N_y,sizeof(int));
06718    if (!yu || !iu)
06719     {
06720       fprintf (SUMA_STDERR,"Error %s: Could not allocate memory", FuncName);
06721       SUMA_RETURN (NULL);
06722    }
06723 
06724    *kunq = 0;
06725    yu[0] = ys[0];
06726    iu[0] = 0;
06727    for (k=1;k<N_y;++k)
06728     {
06729       if ((ys[k] != ys[k - 1]))
06730          {
06731             ++*kunq;
06732             yu[*kunq] = ys[k];   
06733             iu[*kunq] = k;
06734          }
06735    }
06736    ++*kunq;
06737    
06738    
06739    
06740    yu = (int *) SUMA_realloc(yu, *kunq*sizeof(int));
06741    iu = (int *) SUMA_realloc(iu, *kunq*sizeof(int));
06742 
06743    *iup = iu;
06744    SUMA_RETURN (yu);
06745 }
06746 
06747 
06748 
06749 
06750 
06751 
06752 
06753 
06754 
06755 
06756 
06757 
06758 
06759 
06760 
06761 
06762 int *SUMA_reorder(int *y, int *isort, int N_isort)
06763 {
06764    static char FuncName[]={"SUMA_reorder"};
06765    int i = 0, *yr = NULL;
06766    
06767    SUMA_ENTRY;
06768    
06769    if (!y || !isort || N_isort <= 0) SUMA_RETURN(yr);
06770    
06771    yr = (int *)SUMA_calloc( N_isort, sizeof(int));
06772    if (!yr) SUMA_RETURN(yr);
06773    
06774    for (i=0; i<N_isort; ++i) yr[i] = y[isort[i]];
06775    
06776    SUMA_RETURN(yr);
06777 }
06778 
06779 
06780 
06781 
06782 
06783 
06784 int SUMA_suck_file( char *fname , char **fbuf )
06785 {
06786    static char FuncName[]={"SUMA_suck_file"};
06787    int  fd , ii ;
06788    unsigned long len; 
06789    char * buf ;
06790 
06791    SUMA_ENTRY;
06792    
06793    if( fname == NULL || fname[0] == '\0' || fbuf == NULL ) SUMA_RETURN(0) ;
06794 
06795    len = THD_filesize( fname ) ;
06796    if( len <= 0 ) SUMA_RETURN(0) ;
06797 
06798    buf = (char *) SUMA_malloc( sizeof(char) * (len+4) ) ;
06799    if( buf == NULL ) SUMA_RETURN(0) ;
06800 
06801    fd = open( fname , O_RDONLY ) ;
06802    if( fd < 0 ) SUMA_RETURN(0) ;
06803 
06804    ii = read( fd , buf , len ) ;
06805    close( fd ) ;
06806    if( ii <= 0 ){ SUMA_free(buf) ; SUMA_RETURN(0); }
06807    *fbuf = buf ; SUMA_RETURN(ii) ;
06808 }
06809 
06810 
06811 
06812 
06813 
06814 char * SUMA_file_suck( char *fname , int *nread )
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 }