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 }