1 /*#define TEST*/
2 /*#define DEBUG_3*/
3 #ifdef DEBUG_1
4 #define DEBUG_2
5 #define DEBUG_3
6 #endif
7
8 /* Header FILES */
9
10 #include "SUMA_suma.h"
11
12 /* CODE */
13
14
15
16 /*!**
17 File : SUMA_SurfNorm.c
18 \author Ziad Saad
19 Date : Thu Jan 3 14:46:55 EST 2002
20
21 Purpose :
22 Calculate the polygon and node normals making up a surface
23 It is assumed that FaceSetList does NOT contain an index >= N_NodeList
24
25 Usage :
26 RetStrct = SUMA_SurfNorm (NodeList, N_NodeList, FaceSetList, N_FaceSetList );
27
28
29 Input paramters :
30 \param NodeList (float *): N_NodeList x 3 vector containing the XYZ coordinates of each node
31 \param N_NodeList (int): number of nodes in NodeList
32 \param FaceSetList (int *): [N_FaceSetList x 3] vector (matrix prior to SUMA 1.2) of node indices forming each triangular facet (FaceSets).
33 The indices must be relative to NodeList
34 \param N_FaceSetList (int): 1st dimension of FaceSetList
35
36
37 Returns :
38 \return RetStrct, a structure of the type SUMA_SURF_NORM with the following fields
39 N_Node (int) Number of nodes, 1st dim of NodeNormList
40 N_Face (int) Number of facesets, 1st dim of FaceNormList
41 FaceNormList (float *) N_Face x 3 vector (was matrix prior to SUMA 1.2) containing normalized normal vectors for each triangular faceset
42 NodeNormList (float *) N_Node x 3 vector (was matrix prior to SUMA 1.2) containing normalized normal vectors for each node
43
44 Support :
45 \sa SUMA_MemberFaceSets
46 \sa
47
48 Side effects :
49 to free memory from RetStrct
50 if (RetStrct.FaceNormList) SUMA_free(RetStrct.FaceNormList);
51 if (RetStrct.NodeNormList) SUMA_free(RetStrct.NodeNormList);
52
53 The node normals are obtained by averaging the normals of the surrounding facesets (ie facesets containing the node).
54 This is an approximation and in special cases like a cube with a particular tesselation, the normals could be biased
55 in direction if the node is a member of more triangles on one side of the object than the other sides.
56 See Lab-book NIH-2, page 142 for an illustration or this miserable ascii rendition.
57
58 -------------------
59 * ## * Here, Node I will have its normal biased towards the direction of N2 and N3
60 * ## * | because node I is part of two triangles from one side and one triangle in
61 * ## N1 * | each of the other 2 sides.
62 * ## * |
63 *## I * N3 |
64 ------------------ * |
65 |## |######|
66 | ## N4 | *
67 | ## | N2 *
68 | ## | *
69 | ## | *
70 | ## | *
71 ------------------
72 ***/
73 static int surf_norm_quiet = 0;
set_surf_norm_quiet(int v)74 void set_surf_norm_quiet(int v) { surf_norm_quiet=v; }
SUMA_SurfNorm(float * NodeList,int N_NodeList,int * FaceSetList,int N_FaceSetList)75 SUMA_SURF_NORM SUMA_SurfNorm (float *NodeList, int N_NodeList,
76 int *FaceSetList, int N_FaceSetList )
77 {/*SUMA_SurfNorm*/
78 static char stmp[200], FuncName[]={"SUMA_SurfNorm"};
79 float d1[3], d2[3], d, nrm;
80 SUMA_SURF_NORM RetStrct;
81 int *Index, *N_Memb, i, j, maxind, NotMember, id, id2, ND, ip, NP;
82 SUMA_Boolean LocalHead = NOPE;
83
84 SUMA_ENTRY;
85
86 ND = 3;
87 NP = 3;
88 RetStrct.N_Node = N_NodeList; /* redundant, for clean up outside function */
89 RetStrct.N_Face = N_FaceSetList;
90
91 /* allocate space */
92 if (LocalHead)
93 fprintf(SUMA_STDERR,"%s: %d %d\n", FuncName, N_NodeList, N_FaceSetList);
94 RetStrct.FaceNormList =
95 (float *)SUMA_calloc (N_FaceSetList * NP, sizeof(float));
96 RetStrct.NodeNormList = (float *)SUMA_calloc (N_NodeList * ND, sizeof(float));
97 Index = (int *)SUMA_calloc (N_NodeList, sizeof(int));
98 N_Memb = (int *)SUMA_calloc (N_NodeList, sizeof(int));
99 if (!RetStrct.FaceNormList || !RetStrct.NodeNormList || !Index || !N_Memb)
100 {
101 SUMA_alloc_problem (FuncName);
102 SUMA_RETURN (RetStrct);
103 }
104
105 /* calculate and normalize triangle normals */
106 maxind = N_NodeList -1;
107 for (i=0; i < N_FaceSetList; i++) {
108 ip = NP * i;
109 for (j=0; j < 3; j++) {
110 d1[j] = NodeList[(ND*FaceSetList[ip])+j] -
111 NodeList[(ND*FaceSetList[ip+1])+j];
112 d2[j] = NodeList[(ND*FaceSetList[ip+1])+j] -
113 NodeList[(ND*FaceSetList[ip+2])+j];
114 }
115 RetStrct.FaceNormList[ip ] = d1[1]*d2[2] - d1[2]*d2[1];
116 RetStrct.FaceNormList[ip+1] = d1[2]*d2[0] - d1[0]*d2[2];
117 RetStrct.FaceNormList[ip+2] = d1[0]*d2[1] - d1[1]*d2[0];
118 d = sqrt(RetStrct.FaceNormList[ip]*RetStrct.FaceNormList[ip]+
119 RetStrct.FaceNormList[ip+1]*RetStrct.FaceNormList[ip+1]+
120 RetStrct.FaceNormList[ip+2]*RetStrct.FaceNormList[ip+2]);
121 if (d == 0.0) {
122 /* I used to return here with a nasty message,
123 but it seems that FreeSurfer surfaces contain such conditions
124 So, now I just set the normal to 1.0 in that case */
125 /*SUMA_error_message (FuncName,"Zero length vector, returning",1);
126 if (RetStrct.FaceNormList) SUMA_free(RetStrct.FaceNormList);
127 if (RetStrct.NodeNormList) SUMA_free(RetStrct.NodeNormList);
128 if (Index) SUMA_free(Index);
129 if (N_Memb) SUMA_free(N_Memb);
130 SUMA_RETURN (RetStrct); */
131 RetStrct.FaceNormList[ip] = 1.0;
132 RetStrct.FaceNormList[ip+1] = 1.0;
133 RetStrct.FaceNormList[ip+2] = 1.0;
134
135 } else {
136 RetStrct.FaceNormList[ip] /= d;
137 RetStrct.FaceNormList[ip+1] /= d;
138 RetStrct.FaceNormList[ip+2] /= d;
139 }
140
141 #if 0
142 /* a test for the function SUMA_TriNorm */
143 {
144 float test_norm[3];
145 SUMA_TriNorm (&(NodeList[(ND*FaceSetList[ip])]), &(NodeList[(ND*FaceSetList[ip+1])]), &(NodeList[(ND*FaceSetList[ip+2])]), test_norm);
146 if (test_norm[0] != RetStrct.FaceNormList[ip] || test_norm[1] != RetStrct.FaceNormList[ip+1] || test_norm[2] != RetStrct.FaceNormList[ip+2]) {
147 fprintf (SUMA_STDERR, "Error %s: Test of SUMA_TriNorm failed, difference in norms. Exiting.\n", FuncName);
148 exit(1);
149 }
150 fprintf (SUMA_STDERR,".");
151 }
152
153 #endif
154
155 /*each node making up the FaceSet will get its normal vector updated*/
156 if (FaceSetList[ip] > maxind || FaceSetList[ip+1] > maxind || FaceSetList[ip+2] > maxind) {
157 SUMA_error_message (FuncName,"FaceSetList contains indices >= N_NodeList",1);
158 if (RetStrct.FaceNormList) SUMA_free(RetStrct.FaceNormList);
159 if (RetStrct.NodeNormList) SUMA_free(RetStrct.NodeNormList);
160 if (Index) SUMA_free(Index);
161 if (N_Memb) SUMA_free(N_Memb);
162 SUMA_RETURN (RetStrct);
163 }
164
165
166 id2 = ND * FaceSetList[ip];
167 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
168 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
169 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
170 ++N_Memb[FaceSetList[ip]];
171 id2 = ND * FaceSetList[ip+1];
172 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
173 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
174 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
175 ++N_Memb[FaceSetList[ip+1]];
176 id2 = ND * FaceSetList[ip+2];
177 RetStrct.NodeNormList[id2] += RetStrct.FaceNormList[ip];
178 RetStrct.NodeNormList[id2+1] += RetStrct.FaceNormList[ip+1];
179 RetStrct.NodeNormList[id2+2] += RetStrct.FaceNormList[ip+2];
180 ++N_Memb[FaceSetList[ip+2]];
181 }
182 SUMA_LH("Normalizing");
183 /*Now normalize NodeNormList*/
184 NotMember = 0;
185 for (i=0; i<N_NodeList; ++i)
186 {
187 id = ND * i;
188 if (N_Memb[i])
189 {
190 /* fprintf(SUMA_STDERR,"%s: Node %d, Normal (pre scale) %f %f %f\n", FuncName, i, RetStrct.NodeNormList[id], RetStrct.NodeNormList[id+1], RetStrct.NodeNormList[id+2]); */
191 RetStrct.NodeNormList[id] /= N_Memb[i];
192 RetStrct.NodeNormList[id+1] /= N_Memb[i];
193 RetStrct.NodeNormList[id+2] /= N_Memb[i];
194 /* normalize */
195 nrm = sqrt(RetStrct.NodeNormList[id]*RetStrct.NodeNormList[id] + RetStrct.NodeNormList[id+1]*RetStrct.NodeNormList[id+1] + RetStrct.NodeNormList[id+2]*RetStrct.NodeNormList[id+2]);
196 if (nrm) { /* at times nrm is 0. This happened once on a flat surface that was not consistently wound. Nodes that were
197 members of two triangles of opposed normals ended up with a normal (and nrm) of 0 */
198 RetStrct.NodeNormList[id] /= nrm;
199 RetStrct.NodeNormList[id+1] /= nrm;
200 RetStrct.NodeNormList[id+2] /= nrm;
201 }
202
203 /* fprintf(SUMA_STDERR,"%s: Node %d, N_Memb[i] = %d, nrm = %f\n"
204 , FuncName, i, N_Memb[i], nrm); */
205 }
206 else
207 {
208 ++NotMember;
209 /*
210 For some patches (Free Surfer's) such conditions
211 are frequent, spitting out that message becomes old too quick */
212 /*
213 fprintf(stdout,
214 "\n%s Warning: Node %d is not a member of any FaceSets,"
215 " returning unit vector as normal.\n", FuncName, i);
216 */
217 RetStrct.NodeNormList[id] =
218 RetStrct.NodeNormList[id+1] =
219 RetStrct.NodeNormList[id+2] = 1.0;
220 }
221 }
222 if (NotMember && !surf_norm_quiet) {
223 sprintf (stmp, "(IGNORE for surfaces with cuts\n"
224 "%d nodes (%f%% of total) are\n"
225 "not members of any FaceSets.\n"
226 "Their normals are set to the\n"
227 "unit vector.\n",
228 NotMember, (float)NotMember/(float)N_NodeList*100.0);
229 SUMA_SL_Note("%s",stmp);
230 }
231
232 if (N_Memb) SUMA_free(N_Memb);
233 if (Index) SUMA_free(Index);
234 SUMA_RETURN (RetStrct);
235 }/*SUMA_SurfNorm*/
236
237 /*!
238 Try to guess the direction of the surface normals
239 0 = dunno
240 1 = out
241 -1 = inwards
242
243 */
SUMA_SurfNormDir(SUMA_SurfaceObject * SO)244 int SUMA_SurfNormDir (SUMA_SurfaceObject *SO)
245 {
246 static char FuncName[]={"SUMA_SurfNormDir"};
247 int in, cntneg, cntpos;
248 float *a, *b, dot, d, U[3];
249 SUMA_Boolean LocalHead = NOPE;
250
251 SUMA_ENTRY;
252
253 if (!SO->N_Node) {
254 SUMA_S_Err("No Nodes!");
255 SUMA_RETURN(0);
256 }
257 if (!SO->NodeNormList) {
258 SUMA_RECOMPUTE_NORMALS(SO);
259 }
260
261 cntneg = 0;
262 cntpos = 0;
263
264 for (in=0; in<SO->N_Node; ++in) {
265 a = &(SO->NodeList[3*in]);
266 SUMA_UNIT_VEC(SO->Center, a , U, d); /* original distance from center */
267 b = &(SO->NodeNormList[3*in]);
268 SUMA_DOTP_VEC(U, b, dot, 3, float, float); /* dot product with normal at node*/
269 if (dot < 0) {
270 ++cntneg;
271 } else {
272 ++cntpos;
273 }
274 }
275
276 if (cntneg < cntpos) {
277 if (LocalHead) fprintf(SUMA_STDERR,"%s: %.2f%% sure that normals point outwards.\n", FuncName, (cntpos-cntneg)/(float)SO->N_Node*100.0);
278 SUMA_RETURN(1);
279 } else if (cntneg > cntpos) {
280 if (LocalHead) fprintf(SUMA_STDERR,"%s: %.2f%% sure that normals point inwards.\n", FuncName, (cntneg-cntpos)/(float)SO->N_Node*100.0);
281 SUMA_RETURN(-1);
282 } else {
283 if (LocalHead) fprintf(SUMA_STDERR,"%s: Not sure where node normals point!\n", FuncName);
284 SUMA_RETURN(0);
285 }
286
287 SUMA_RETURN(0);
288 }
289
290
291
292 /*!**
293 File : SUMA_MemberFaceSets.c
294 \author Ziad Saad
295 Date : Thu Jan 3 12:01:49 EST 2002
296
297 Purpose :
298 Determines for each Node Index, from 0 to Nind-1 the indices of the FaceSets to which the node belongs.
299
300 Usage :
301 RetStruct = SUMA_MemberFaceSets (int Nind, int * FaceSetList, int nFace, int FaceSetDim);
302
303 Input paramters :
304 \param Nind (int) : Total number of nodes in Index. ( No value (node index) in FaceSetList can be > Nind-1 )
305 \param FaceSetList (int *) : The FaceSetList vector (used to be matrix, prior to SUMA 1.2), [nFace x FaceSetDim]
306 \param nFace (int) : number of FaceSets in FaceSetList
307 \param FaceSetDim (int): column (2nd) dimension of FaceSetList (usually 3 or 4)
308
309 Returns :
310 \return RetStruct : a SUMA_MEMBER_FACE_SETS type structure with the following fields
311 NodeMemberOfFaceSet (int **) a matrix containing the indices of the FaceSets containing the nodes.
312 Each row i in NodeMemberOfFaceSet lists all FaceSets containing the
313 node i. The list is delimited by a -1 entry
314 for all but the node that has a max of N_Memb_max members.
315 the size of NodeMemberOfFaceSet is [Nind x N_Memb_max];
316 N_Memb (int *) the number of 1st order neighbors for each node in Index
317 N_Memb_max (int) the second dimension of NodeMemberOfFaceSet, also represents
318 the maximum number of FaceSets that
319 any node belongs to.
320 Nnode (int) : it's equal to Nind but it is kept here so that the structure can contain enough info to
321 free memory that is allocated to it.
322
323
324 if (RetStruct.NodeMemberOfFaceSet == NULL) then an error occurred in the function
325
326 Support :
327 \sa MemberFaceSets.c and .m for methods that work when the search is required for a few nodes only
328
329 Example:
330 FaceSetList is:
331 1 4 6
332 6 4 2
333 6 1 3
334
335 and Nind is:
336 7
337
338 Then RetStruct.NodeMemberOfFaceSet is:
339 -1 0 0 ( Node 0 is not used in FaceSetList )
340 0 2 -1 (Node 1 is used in Facesets 0 & 2 )
341 1 -1 0 (Node 2 is used in Faceset 1)
342 2 -1 0 (Node 3 is used in Faceset 2)
343 0 1 -1 (Node 4 is used in Facesets 0 & 1)
344 -1 0 0 (Node 5 in not used in FaceSetList)
345 0 1 2 (Node 6 is used in all Facesets)
346
347 and RetStruct.N_Memb is:
348 0 2 1 1 2 0 3
349
350 Side effects :
351 To free RetStruct, use:
352 if (RetStruct.NodeMemberOfFaceSet) SUMA_free2D((char **)RetStruct.NodeMemberOfFaceSet, Nind);
353 if (RetStruct.N_Memb) SUMA_free(RetStruct.N_Memb);
354 if (RetStruct) SUMA_free(RetStrct);
355
356
357 ***/
SUMA_MemberFaceSets(int Nind,int * FaceSetList,int nFr,int FaceDim,char * ownerid)358 SUMA_MEMBER_FACE_SETS *SUMA_MemberFaceSets (int Nind, int * FaceSetList,
359 int nFr , int FaceDim, char *ownerid)
360 {/*SUMA_MemberFaceSets*/
361 static char FuncName[]={"SUMA_MemberFaceSets"};
362 SUMA_MEMBER_FACE_SETS *RetStrct;
363 int **tmpMember;
364 int i, inode, iface, ip , NP;
365 byte *filled=NULL;
366 SUMA_Boolean LocalHead = NOPE;
367
368 SUMA_ENTRY;
369
370 SUMA_LHv("Nind=%d, nFr=%d, FaceDim=%d\n", Nind, nFr, FaceDim);
371
372 NP = FaceDim;
373 RetStrct = (SUMA_MEMBER_FACE_SETS *)
374 SUMA_malloc(sizeof(SUMA_MEMBER_FACE_SETS));
375 RetStrct->idcode_str = NULL;
376 SUMA_NEW_ID(RetStrct->idcode_str, NULL);
377 RetStrct->N_links = 0;
378 if (ownerid) sprintf(RetStrct->owner_id, "%s", ownerid);
379 else RetStrct->owner_id[0] = '\0';
380 RetStrct->LinkedPtrType = SUMA_LINKED_MEMB_FACE_TYPE;
381 RetStrct->do_type = not_DO_type;
382
383 RetStrct->N_Memb_max = RetStrct->Nnode = 0;
384 RetStrct->N_Memb = NULL;
385 RetStrct->NodeMemberOfFaceSet = NULL;
386
387 /* Allocate return variables */
388 tmpMember = (int **) SUMA_allocate2D (Nind,
389 SUMA_MAX_MEMBER_FACE_SETS ,sizeof(int));
390 RetStrct->N_Memb = (int *) SUMA_calloc (Nind, sizeof(int));
391 filled = (byte *)SUMA_calloc(Nind, sizeof(int));
392 if (!tmpMember || !RetStrct->N_Memb)
393 {
394 SUMA_S_Err("Failed to allocate for tmpMember or RetStrct->N_Memb");
395 SUMA_RETURN (RetStrct);
396 }
397
398 /* loop through all facesets and tag nodes that make up FaceSets*/
399 for (iface=0; iface<nFr; ++iface) {/*iface*/
400 i = 0;
401 ip = NP * iface;
402 do {
403 inode = FaceSetList[ip + i];
404 if (inode > Nind) {
405 SUMA_S_Err("FaceSetList contains node indices >= Nind");
406 SUMA_RETURN (RetStrct);
407 }
408 if (RetStrct->N_Memb[inode] >= SUMA_MAX_MEMBER_FACE_SETS) {
409 if (!filled[inode]){
410 SUMA_S_Errv("Node %d is member of (%d FaceSets) more than"
411 " SUMA_MAX_MEMBER_FACE_SETS (%d)\n"
412 "Skipping remaining facets\n",
413 inode, RetStrct->N_Memb[inode], SUMA_MAX_MEMBER_FACE_SETS);
414 }
415 filled[inode]=1;
416 } else {
417 tmpMember[inode][RetStrct->N_Memb[inode]] = iface;
418 ++RetStrct->N_Memb[inode];
419 }
420 if (RetStrct->N_Memb[inode] > RetStrct->N_Memb_max)
421 RetStrct->N_Memb_max = RetStrct->N_Memb[inode];
422 ++i;
423 } while (i < FaceDim);
424 }
425
426 /*allocate just enough for returning variables */
427 RetStrct->NodeMemberOfFaceSet =
428 (int **) SUMA_allocate2D (Nind, RetStrct->N_Memb_max ,sizeof(int));
429 if (!RetStrct->NodeMemberOfFaceSet)
430 {
431 SUMA_S_Err("Failed to allocate for RetStrct->NodeMemberOfFaceSet\n");
432 SUMA_RETURN (RetStrct);
433 }
434
435 /* loop through all nodes,
436 cp results into RetStrct->NodeMemberOfFaceSet and seal with -1 */
437 for (inode = 0; inode < Nind; ++inode) {
438 i = 0;
439 while (i < RetStrct->N_Memb[inode]) {
440 RetStrct->NodeMemberOfFaceSet[inode][i] = tmpMember[inode][i];
441 ++i;
442 }
443 /*seal with -1 */
444 if (RetStrct->N_Memb[inode] < RetStrct->N_Memb_max)
445 RetStrct->NodeMemberOfFaceSet[inode][i] = -1;
446 }
447
448 /* Clean up time */
449 if (tmpMember) SUMA_free2D((char **)tmpMember, Nind);
450 if (filled) SUMA_free(filled); filled = NULL;
451
452 RetStrct->Nnode = Nind;
453 SUMA_RETURN (RetStrct);
454
455 }/*SUMA_MemberFaceSets*/
456
457 /*! Free a SUMA_MEMBER_FACE_SETS structure */
458
SUMA_Free_MemberFaceSets(SUMA_MEMBER_FACE_SETS * MF)459 SUMA_Boolean SUMA_Free_MemberFaceSets (SUMA_MEMBER_FACE_SETS *MF)
460 {
461 static char FuncName[]={"SUMA_Free_MemberFaceSets"};
462 SUMA_Boolean LocalHead = NOPE;
463
464 SUMA_ENTRY;
465 if (!MF) { SUMA_RETURN (YUP); }
466 if (MF->N_links) {
467 SUMA_LH("Just a link release");
468 MF = (SUMA_MEMBER_FACE_SETS *)SUMA_UnlinkFromPointer((void *)MF);
469 SUMA_RETURN (YUP);
470 }
471
472 SUMA_LH("No more links, here we go");
473 if (MF->idcode_str) SUMA_free(MF->idcode_str);
474 if (MF->NodeMemberOfFaceSet) SUMA_free2D((char **)MF->NodeMemberOfFaceSet, MF->Nnode);
475 if (MF->N_Memb) SUMA_free(MF->N_Memb);
476 if (MF) SUMA_free(MF);
477 SUMA_RETURN (YUP);
478 }
479