00001
00002
00003 #ifdef DEBUG_1
00004 #define DEBUG_2
00005 #define DEBUG_3
00006 #endif
00007
00008
00009
00010 #include "SUMA_suma.h"
00011
00012 #ifdef SUMA_SurfNorm_STAND_ALONE
00013
00014 SUMA_SurfaceViewer *SUMAg_cSV;
00015 SUMA_SurfaceViewer *SUMAg_SVv = NULL;
00016
00017 int SUMAg_N_SVv = 0;
00018 SUMA_DO *SUMAg_DOv;
00019 int SUMAg_N_DOv = 0;
00020 SUMA_CommonFields *SUMAg_CF;
00021 #else
00022 extern SUMA_CommonFields *SUMAg_CF;
00023 #endif
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 SUMA_SURF_NORM SUMA_SurfNorm (float *NodeList, int N_NodeList, int *FaceSetList, int N_FaceSetList )
00088 {
00089 static char stmp[200], FuncName[]={"SUMA_SurfNorm"};
00090 float d1[3], d2[3], d, nrm;
00091 SUMA_SURF_NORM RetStrct;
00092 int *Index, *N_Memb, i, j, maxind, NotMember, id, id2, ND, ip, NP;
00093 SUMA_Boolean LocalHead = NOPE;
00094
00095 SUMA_ENTRY;
00096
00097 ND = 3;
00098 NP = 3;
00099 RetStrct.N_Node = N_NodeList;
00100 RetStrct.N_Face = N_FaceSetList;
00101
00102
00103 if (LocalHead) fprintf(SUMA_STDERR,"%s: %d %d\n", FuncName, N_NodeList, N_FaceSetList);
00104 RetStrct.FaceNormList = (float *)SUMA_calloc (N_FaceSetList * NP, sizeof(float));
00105 RetStrct.NodeNormList = (float *)SUMA_calloc (N_NodeList * ND, sizeof(float));
00106 Index = (int *)SUMA_calloc (N_NodeList, sizeof(int));
00107 N_Memb = (int *)SUMA_calloc (N_NodeList, sizeof(int));
00108 if (!RetStrct.FaceNormList || !RetStrct.NodeNormList || !Index || !N_Memb)
00109 {
00110 SUMA_alloc_problem (FuncName);
00111 SUMA_RETURN (RetStrct);
00112 }
00113
00114
00115 maxind = N_NodeList -1;
00116 for (i=0; i < N_FaceSetList; i++) {
00117 ip = NP * i;
00118 for (j=0; j < 3; j++) {
00119 d1[j] = NodeList[(ND*FaceSetList[ip])+j] - NodeList[(ND*FaceSetList[ip+1])+j];
00120 d2[j] = NodeList[(ND*FaceSetList[ip+1])+j] - NodeList[(ND*FaceSetList[ip+2])+j];
00121 }
00122 RetStrct.FaceNormList[ip] = d1[1]*d2[2] - d1[2]*d2[1];
00123 RetStrct.FaceNormList[ip+1] = d1[2]*d2[0] - d1[0]*d2[2];
00124 RetStrct.FaceNormList[ip+2] = d1[0]*d2[1] - d1[1]*d2[0];
00125 d = sqrt(RetStrct.FaceNormList[ip]*RetStrct.FaceNormList[ip]+RetStrct.FaceNormList[ip+1]*RetStrct.FaceNormList[ip+1]+RetStrct.FaceNormList[ip+2]*RetStrct.FaceNormList[ip+2]);
00126 if (d == 0.0) {
00127
00128
00129
00130
00131
00132
00133
00134
00135 RetStrct.FaceNormList[ip] = 1.0;
00136 RetStrct.FaceNormList[ip+1] = 1.0;
00137 RetStrct.FaceNormList[ip+2] = 1.0;
00138
00139 } else {
00140 RetStrct.FaceNormList[ip] /= d;
00141 RetStrct.FaceNormList[ip+1] /= d;
00142 RetStrct.FaceNormList[ip+2] /= d;
00143 }
00144
00145 #if 0
00146
00147 {
00148 float test_norm[3];
00149 SUMA_TriNorm (&(NodeList[(ND*FaceSetList[ip])]), &(NodeList[(ND*FaceSetList[ip+1])]), &(NodeList[(ND*FaceSetList[ip+2])]), test_norm);
00150 if (test_norm[0] != RetStrct.FaceNormList[ip] || test_norm[1] != RetStrct.FaceNormList[ip+1] || test_norm[2] != RetStrct.FaceNormList[ip+2]) {
00151 fprintf (SUMA_STDERR, "Error %s: Test of SUMA_TriNorm failed, difference in norms. Exiting.\n", FuncName);
00152 exit(1);
00153 }
00154 fprintf (SUMA_STDERR,".");
00155 }
00156
00157 #endif
00158
00159
00160 if (FaceSetList[ip] > maxind || FaceSetList[ip+1] > maxind || FaceSetList[ip+2] > maxind) {
00161 SUMA_error_message (FuncName,"FaceSetList contains indices >= N_NodeList",1);
00162 if (RetStrct.FaceNormList) SUMA_free(RetStrct.FaceNormList);
00163 if (RetStrct.NodeNormList) SUMA_free(RetStrct.NodeNormList);
00164 if (Index) SUMA_free(Index);
00165 if (N_Memb) SUMA_free(N_Memb);
00166 SUMA_RETURN (RetStrct);
00167 }
00168
00169
00170 id2 = ND * FaceSetList[ip];
00171 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
00172 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
00173 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
00174 ++N_Memb[FaceSetList[ip]];
00175 id2 = ND * FaceSetList[ip+1];
00176 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
00177 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
00178 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
00179 ++N_Memb[FaceSetList[ip+1]];
00180 id2 = ND * FaceSetList[ip+2];
00181 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
00182 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
00183 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
00184 ++N_Memb[FaceSetList[ip+2]];
00185 }
00186 SUMA_LH("Normalizing");
00187
00188 NotMember = 0;
00189 for (i=0; i<N_NodeList; ++i)
00190 {
00191 id = ND * i;
00192 if (N_Memb[i])
00193 {
00194
00195 RetStrct.NodeNormList[id] /= N_Memb[i];
00196 RetStrct.NodeNormList[id+1] /= N_Memb[i];
00197 RetStrct.NodeNormList[id+2] /= N_Memb[i];
00198
00199 nrm = sqrt(RetStrct.NodeNormList[id]*RetStrct.NodeNormList[id] + RetStrct.NodeNormList[id+1]*RetStrct.NodeNormList[id+1] + RetStrct.NodeNormList[id+2]*RetStrct.NodeNormList[id+2]);
00200 if (nrm) {
00201
00202 RetStrct.NodeNormList[id] /= nrm;
00203 RetStrct.NodeNormList[id+1] /= nrm;
00204 RetStrct.NodeNormList[id+2] /= nrm;
00205 }
00206
00207
00208 }
00209 else
00210 {
00211 ++NotMember;
00212
00213
00214
00215
00216
00217 RetStrct.NodeNormList[id] = RetStrct.NodeNormList[id+1] = RetStrct.NodeNormList[id+2] = 1.0;
00218 }
00219 }
00220 if (NotMember) {
00221 sprintf (stmp, "(IGNORE for surface patches\n"
00222 "%d nodes (%f%% of total) are\n"
00223 "not members of any FaceSets.\n"
00224 "Their normals are set to the\n"
00225 "unit vector.\n", NotMember, (float)NotMember/(float)N_NodeList*100.0);
00226 SUMA_SL_Warn(stmp);
00227 }
00228
00229 if (N_Memb) SUMA_free(N_Memb);
00230 if (Index) SUMA_free(Index);
00231 SUMA_RETURN (RetStrct);
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 SUMA_MEMBER_FACE_SETS *SUMA_MemberFaceSets (int Nind, int * FaceSetList, int nFr , int FaceDim, char *ownerid)
00304 {
00305 static char FuncName[]={"SUMA_MemberFaceSets"};
00306 SUMA_MEMBER_FACE_SETS *RetStrct;
00307 int **tmpMember;
00308 int i, inode, iface, ip , NP;
00309
00310 SUMA_ENTRY;
00311
00312 NP = FaceDim;
00313 RetStrct = (SUMA_MEMBER_FACE_SETS *)SUMA_malloc(sizeof(SUMA_MEMBER_FACE_SETS));
00314 RetStrct->idcode_str = NULL;
00315 SUMA_NEW_ID(RetStrct->idcode_str, NULL);
00316 RetStrct->N_links = 0;
00317 if (ownerid) sprintf(RetStrct->owner_id, "%s", ownerid);
00318 else RetStrct->owner_id[0] = '\0';
00319 RetStrct->LinkedPtrType = SUMA_LINKED_MEMB_FACE_TYPE;
00320
00321 RetStrct->N_Memb_max = RetStrct->Nnode = 0;
00322 RetStrct->N_Memb = NULL;
00323 RetStrct->NodeMemberOfFaceSet = NULL;
00324
00325
00326 tmpMember = (int **) SUMA_allocate2D (Nind, SUMA_MAX_MEMBER_FACE_SETS ,sizeof(int));
00327 RetStrct->N_Memb = (int *) SUMA_calloc (Nind, sizeof(int));
00328
00329 if (!tmpMember || !RetStrct->N_Memb)
00330 {
00331 fprintf (SUMA_STDERR,"Error %s: Failed to allocate for tmpMember or RetStrct->N_Memb\n", FuncName);
00332 SUMA_RETURN (RetStrct);
00333 }
00334
00335
00336 for (iface=0; iface<nFr; ++iface) {
00337 i = 0;
00338 ip = NP * iface;
00339 do {
00340 inode = FaceSetList[ip + i];
00341 if (inode > Nind) {
00342 fprintf (SUMA_STDERR,"Error %s: FaceSetList contains node indices >= Nind\n", FuncName);
00343 SUMA_RETURN (RetStrct);
00344 }
00345 tmpMember[inode][RetStrct->N_Memb[inode]] = iface;
00346 ++RetStrct->N_Memb[inode];
00347 if (RetStrct->N_Memb[inode] >= SUMA_MAX_MEMBER_FACE_SETS) {
00348 fprintf (SUMA_STDERR,"Error %s: Node %d is member of (%d FaceSets) more than SUMA_MAX_MEMBER_FACE_SETS (%d)\n",\
00349 FuncName, inode, RetStrct->N_Memb[inode], SUMA_MAX_MEMBER_FACE_SETS);
00350 SUMA_RETURN (RetStrct);
00351 }
00352 if (RetStrct->N_Memb[inode] > RetStrct->N_Memb_max) RetStrct->N_Memb_max = RetStrct->N_Memb[inode];
00353 ++i;
00354 } while (i < FaceDim);
00355 }
00356
00357
00358 RetStrct->NodeMemberOfFaceSet = (int **) SUMA_allocate2D (Nind, RetStrct->N_Memb_max ,sizeof(int));
00359 if (!RetStrct->NodeMemberOfFaceSet)
00360 {
00361 fprintf(SUMA_STDERR,"Error %s: Failed to allocate for RetStrct->NodeMemberOfFaceSet\n", FuncName);
00362 SUMA_RETURN (RetStrct);
00363 }
00364
00365
00366 for (inode = 0; inode < Nind; ++inode) {
00367 i = 0;
00368 while (i < RetStrct->N_Memb[inode]) {
00369 RetStrct->NodeMemberOfFaceSet[inode][i] = tmpMember[inode][i];
00370 ++i;
00371 }
00372
00373 if (RetStrct->N_Memb[inode] < RetStrct->N_Memb_max) RetStrct->NodeMemberOfFaceSet[inode][i] = -1;
00374 }
00375
00376
00377 if (tmpMember) SUMA_free2D((char **)tmpMember, Nind);
00378
00379 RetStrct->Nnode = Nind;
00380 SUMA_RETURN (RetStrct);
00381
00382 }
00383
00384
00385
00386 SUMA_Boolean SUMA_Free_MemberFaceSets (SUMA_MEMBER_FACE_SETS *MF)
00387 {
00388 static char FuncName[]={"SUMA_Free_MemberFaceSets"};
00389 SUMA_Boolean LocalHead = NOPE;
00390
00391 SUMA_ENTRY;
00392 if (!MF) { SUMA_RETURN (YUP); }
00393 if (MF->N_links) {
00394 SUMA_LH("Just a link release");
00395 MF = (SUMA_MEMBER_FACE_SETS *)SUMA_UnlinkFromPointer((void *)MF);
00396 SUMA_RETURN (YUP);
00397 }
00398
00399 SUMA_LH("No more links, here we go");
00400 if (MF->idcode_str) SUMA_free(MF->idcode_str);
00401 if (MF->NodeMemberOfFaceSet) SUMA_free2D((char **)MF->NodeMemberOfFaceSet, MF->Nnode);
00402 if (MF->N_Memb) SUMA_free(MF->N_Memb);
00403 if (MF) SUMA_free(MF);
00404 SUMA_RETURN (YUP);
00405 }
00406 #ifdef TEST_SUMA_MemberFaceSets
00407 void usage ()
00408
00409 {
00410 printf ("\nUsage: SUMA_MemberFaceSets <FaceSetList> <indexfile> \n");
00411 printf ("\t <FaceSetList> : file containing the facesetlist \n");
00412 printf ("\t <index file> : file containing the indices of the nodes \n");
00413 printf ("\t You're looking up wich node belongs to which FaceSets\n\n");
00414 printf ("\t\t\t Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \tThu Jan 3 12:01:49 EST 2002 \n");
00415 exit (0);
00416 }
00417
00418 int main (int argc,char *argv[])
00419 {
00420 char FuncName[100];
00421 int n, nind, **X;
00422 SUMA_MEMBER_FACE_SETS *RetStrct;
00423
00424 X = (int **) SUMA_allocate2D(3,3, sizeof(int));
00425 X[0][0] = 1; X[0][1]= 4; X[0][2]= 6;
00426 X[1][0] = 6; X[1][1]= 4; X[1][2]= 2;
00427 X[2][0] = 6; X[2][1]= 1; X[2][2]= 3;
00428
00429
00430
00431 sprintf (FuncName,"SUMA_MemberFaceSets-Main-");
00432
00433
00434 n = 3;
00435 nind = 7;
00436
00437 RetStrct = SUMA_MemberFaceSets (nind, X, n , n);
00438
00439 printf ("\nMember FaceSets :\n");
00440 SUMA_disp_dmat (RetStrct->NodeMemberOfFaceSet, nind, RetStrct->N_Memb_max, 1);
00441
00442 printf ("\n# of Member FaceSets :\n");
00443 SUMA_disp_dvect (RetStrct->N_Memb , RetStrct->Nnode);
00444
00445 if (RetStrct->N_Memb)
00446 {
00447 if (!SUMA_Free_MemberFaceSets (RetStrct)) {
00448 fprintf(SUMA_STDERR,"Error %s : Failed to free RetStrct", FuncName);
00449 exit(1);
00450 }
00451 }
00452 SUMA_free2D((char **)X, 3);
00453 SUMA_RETURN (0);
00454 }
00455 #endif
00456
00457 #ifdef SUMA_SurfNorm_STAND_ALONE
00458 void usage ()
00459
00460 {
00461 printf ("\nUsage: SUMA_SurfNorm <NodeList> <FaceSetList> \n");
00462 printf ("\t ..... \n\n");
00463 printf ("\t\t\t Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \tThu Jan 3 14:46:55 EST 2002 \n");
00464 exit (0);
00465 }
00466
00467 int main (int argc,char *argv[])
00468 {
00469 static char FuncName[]={"SUMA_SurfNorm-Main-"};
00470 SUMA_SURF_NORM RetStrct;
00471 int nface, nnode;
00472 int *FaceSetList;
00473 float *NodeList;
00474 SUMA_Boolean LocalHead = NOPE;
00475
00476
00477 if (LocalHead) fprintf (SUMA_STDERR,"%s: Calling SUMA_Create_CommonFields ...\n", FuncName);
00478
00479 SUMAg_CF = SUMA_Create_CommonFields ();
00480 if (SUMAg_CF == NULL) {
00481 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_Create_CommonFields\n", FuncName);
00482 exit(1);
00483 }
00484 if (LocalHead) fprintf (SUMA_STDERR,"%s: SUMA_Create_CommonFields Done.\n", FuncName);
00485
00486
00487
00488 if (argc < 3)
00489 {
00490 usage ();
00491 exit (1);
00492 }
00493
00494 nnode = SUMA_float_file_size(argv[1]);
00495 nnode /= 3;
00496 nface = SUMA_float_file_size(argv[2]);
00497 nface /= 3;
00498
00499 FaceSetList = (int *) SUMA_calloc (nface * 3, sizeof(int));
00500 NodeList = (float *) SUMA_calloc (nnode * 3, sizeof(float));
00501
00502 if (!FaceSetList || !NodeList)
00503 {
00504 SUMA_alloc_problem(FuncName);
00505 exit (0);
00506 }
00507
00508 SUMA_Read_dfile (FaceSetList, argv[2], nface * 3);
00509 printf ("\nFaceSets (%d):\n", nface);
00510 SUMA_disp_vecdmat (FaceSetList, nface, 3, 1);
00511
00512 SUMA_Read_file (NodeList, argv[1], 3*nnode);
00513
00514 printf ("\nNodes : (%d)\n", nnode);
00515 SUMA_disp_vecmat (NodeList, nnode, 3, 1, SUMA_ROW_MAJOR, NULL, 0);
00516
00517 RetStrct = SUMA_SurfNorm(NodeList, nnode, FaceSetList, nface );
00518 printf ("\nNode Norms:\n");
00519 SUMA_disp_vecmat (RetStrct.NodeNormList, RetStrct.N_Node, 3, 1, SUMA_ROW_MAJOR, NULL, 0);
00520 printf ("\nFaceSet Norms:\n");
00521 SUMA_disp_vecmat (RetStrct.FaceNormList, RetStrct.N_Face, 3, 1, SUMA_ROW_MAJOR, NULL, 0);
00522
00523
00524 if (RetStrct.FaceNormList) SUMA_free(RetStrct.FaceNormList);
00525 if (RetStrct.NodeNormList) SUMA_free(RetStrct.NodeNormList);
00526 if (FaceSetList) SUMA_free(FaceSetList);
00527 if (NodeList) SUMA_free(NodeList);
00528 SUMA_RETURN (0);
00529 }
00530 #endif