Bug Summary

File:src/limn/polymod.c
Location:line 1468, column 5
Description:Potential leak of memory pointed to by 'icntNew'

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6 Copyright (C) 2012, 2011, 2010, 2009, 2008 Thomas Schultz
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public License
10 (LGPL) as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12 The terms of redistributing and/or modifying this software also
13 include exceptions to the LGPL that facilitate static linking.
14
15 This library is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 Lesser General Public License for more details.
19
20 You should have received a copy of the GNU Lesser General Public License
21 along with this library; if not, write to Free Software Foundation, Inc.,
22 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23*/
24
25
26#include "limn.h"
27
28/*
29** determines intersection of elements of srcA and srcB.
30** assumes:
31** - there are no repeats in either list
32** - dstC is allocated for at least as long as the longer of srcA and srcB
33*/
34static unsigned int
35flipListIntx(unsigned int *dstC,
36 const unsigned int *_srcA, const unsigned int *_srcB) {
37 const unsigned int *srcA, *srcB;
38 unsigned int numA, numB, numC, idxA, idxB;
39
40 numA = _srcA[0];
41 srcA = _srcA + 1;
42 numB = _srcB[0];
43 srcB = _srcB + 1;
44 numC = 0;
45 for (idxA=0; idxA<numA; idxA++) {
46 for (idxB=0; idxB<numB; idxB++) {
47 if (srcA[idxA] == srcB[idxB]) {
48 dstC[numC++] = srcA[idxA];
49 }
50 }
51 }
52 return numC;
53}
54
55/*
56** given triangle identified by triIdx,
57** set neighGot[] and neighInfo[][]
58** neighbors are index 0,1,2;
59** neighbor ii is on edge between vert ii and (ii+1)%3
60** neighGot[ii] is non-zero iff there was such a neighbor
61** neighInfo[ii][0]: index of the (triangle) neighbor
62** neighInfo[ii][1], neighInfo[ii][2]: the two vertices shared with neighbor,
63** in the order that the *neighbor* should be traversing them
64*/
65static void
66flipNeighborsGet(Nrrd *nTriWithVert, Nrrd *nVertWithTri,
67 unsigned int neighGot[3], unsigned int neighInfo[3][3],
68 unsigned int *intxBuff, unsigned int triIdx) {
69 /* static const char me[]="flipNeighborsGet"; */
70 unsigned int intxNum, vertA, vertB, neighIdx, maxTriPerVert,
71 *vertWithTri, *triWithVert;
72 int ii;
73
74 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
75 triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data)((unsigned int*)(nTriWithVert->data));
76 maxTriPerVert = nTriWithVert->axis[0].size - 1;
77 for (ii=0; ii<3; ii++) {
78 vertA = (vertWithTri + 3*triIdx)[ii];
79 vertB = (vertWithTri + 3*triIdx)[AIR_MOD(ii+1, 3)((ii+1)%(3) >= 0 ? (ii+1)%(3) : 3 + (ii+1)%(3))];
80 /*
81 fprintf(stderr, "!%s: %u edge %u: vert{A,B} = %u %u\n", me,
82 triIdx, ii, vertA, vertB);
83 */
84 /* find the intersection of the sets of {triangles using vertA}
85 and {triangles using vertB}: for reasonable surfaces should
86 be either 0 or 2 triangles, and if its 2, then triIdx
87 should be one of them */
88 intxNum = flipListIntx(intxBuff,
89 triWithVert + (1+maxTriPerVert)*vertA,
90 triWithVert + (1+maxTriPerVert)*vertB);
91 if (2 == intxNum) {
92 neighIdx = intxBuff[0];
93 if (neighIdx == triIdx) {
94 neighIdx = intxBuff[1];
95 }
96 neighGot[ii] = AIR_TRUE1;
97 neighInfo[ii][0] = neighIdx;
98 neighInfo[ii][1] = vertB;
99 neighInfo[ii][2] = vertA;
100 } else {
101 neighGot[ii] = AIR_FALSE0;
102 }
103 }
104 return;
105}
106
107/*
108** determines if triIdx needs to be flipped, given that it should
109** be seeing vertices vertA and vertB in that order
110*/
111static int
112flipNeed(Nrrd *nVertWithTri, unsigned int triIdx,
113 unsigned int vertA, unsigned int vertB) {
114 unsigned int *vertWithTri, vert[3];
115 int ai, bi;
116
117 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
118 ELL_3V_COPY(vert, vertWithTri + 3*triIdx)((vert)[0] = (vertWithTri + 3*triIdx)[0], (vert)[1] = (vertWithTri
+ 3*triIdx)[1], (vert)[2] = (vertWithTri + 3*triIdx)[2]);
119 for (ai=0; vert[ai] != vertA; ai++)
120 ;
121 for (bi=0; vert[bi] != vertB; bi++)
122 ;
123 return (1 != AIR_MOD(bi - ai, 3)((bi - ai)%(3) >= 0 ? (bi - ai)%(3) : 3 + (bi - ai)%(3)));
124}
125
126/*
127** this is a weird dual-personality function that is the inner
128** loop of both vertex winding fixing, and the learning stage of
129** vertex splitting
130**
131** for flipping (!splitting)
132** assumes that triIdx was just popped from "okay" stack
133** (triIdx has just been fixed to have correct winding)
134** then goes through the not-yet-done neighbors of triIdx,
135** flipping them if needed, and
136** then adding those neighbors to the stack.
137** returns the number of tris added to stack
138**
139** NOTE: the "flipping" is done within the nVertWithTri representation,
140** but *not* in the limnPolyData itself.
141*/
142static unsigned int
143neighborsCheckPush(Nrrd *nTriWithVert, Nrrd *nVertWithTri,
144 unsigned char *triDone, airArray *okayArr,
145 unsigned int *intxBuff, airArray *splitArr,
146 unsigned int triIdx, int splitting) {
147 /* static const char me[]="neighborsCheckPush"; */
148 unsigned int neighGot[3], neighInfo[3][3], ii, *okay, okayIdx,
149 *vertWithTri, pushedNum;
150
151 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
152 flipNeighborsGet(nTriWithVert, nVertWithTri,
153 neighGot, neighInfo,
154 intxBuff, triIdx);
155 /*
156 for (ii=0; ii<3; ii++) {
157 fprintf(stderr, "!%s: %u neigh[%u]: ", me, triIdx, ii);
158 if (neighGot[ii]) {
159 fprintf(stderr, "%u (%u %u) (done %u)\n",
160 neighInfo[ii][0], neighInfo[ii][1], neighInfo[ii][2],
161 triDone[neighInfo[ii][0]]);
162 } else {
163 fprintf(stderr, "nope\n");
164 }
165 }
166 */
167 pushedNum = 0;
168 for (ii=0; ii<3; ii++) {
169 /* WARNING: complicated logic WRT triDone, splitting, and need */
170 if (neighGot[ii]) {
171 unsigned int tmp, *idxLine, need;
172 if (!splitting) {
173 if (!triDone[neighInfo[ii][0]]) {
174 /* we only take time to learn need if as yet undone */
175 need = flipNeed(nVertWithTri, neighInfo[ii][0],
176 neighInfo[ii][1], neighInfo[ii][2]);
177 if (need) {
178 idxLine = vertWithTri + 3*neighInfo[ii][0];
179 /* here is the vertex winding flip */
180 ELL_SWAP2(idxLine[0], idxLine[1], tmp)((tmp)=(idxLine[0]),(idxLine[0])=(idxLine[1]),(idxLine[1])=(tmp
));
181 }
182 }
183 } else {
184 /* we're here for splitting */
185 /* we have to learn need regardless of done-ness */
186 need = flipNeed(nVertWithTri, neighInfo[ii][0],
187 neighInfo[ii][1], neighInfo[ii][2]);
188 if (need && triDone[neighInfo[ii][0]]) {
189 /* we "need" to flip and yet we've already visited that triangle
190 ==> edge between triIdx and neighInfo[ii][0] needs splitting.
191 See if its a new split, and add it if so */
192 unsigned int *split, splitIdx, splitNum, vert0, vert1;
193 vert0 = AIR_MIN(neighInfo[ii][1], neighInfo[ii][2])((neighInfo[ii][1]) < (neighInfo[ii][2]) ? (neighInfo[ii][
1]) : (neighInfo[ii][2]));
194 vert1 = AIR_MAX(neighInfo[ii][1], neighInfo[ii][2])((neighInfo[ii][1]) > (neighInfo[ii][2]) ? (neighInfo[ii][
1]) : (neighInfo[ii][2]));
195 splitNum = splitArr->len;
196 split = AIR_CAST(unsigned int*, splitArr->data)((unsigned int*)(splitArr->data));
197 for (splitIdx=0; splitIdx<splitNum; splitIdx++) {
198 if (split[2 + 5*splitIdx] == vert0
199 && split[3 + 5*splitIdx] == vert1) {
200 break;
201 }
202 }
203 if (splitIdx == splitNum) {
204 /* this is a new split, add it */
205 /*
206 fprintf(stderr, "!%s: new split(%u,%u) (have %u)\n",
207 me, vert0, vert1, splitArr->len);
208 */
209 splitIdx = airArrayLenIncr(splitArr, 1);
210 split = AIR_CAST(unsigned int*, splitArr->data)((unsigned int*)(splitArr->data));
211 split[0 + 5*splitIdx] = triIdx;
212 split[1 + 5*splitIdx] = neighInfo[ii][0];
213 split[2 + 5*splitIdx] = vert0;
214 split[3 + 5*splitIdx] = vert1;
215 split[4 + 5*splitIdx] = AIR_FALSE0;
216 }
217 }
218 }
219 /* regardless of splitting, we push onto the okay stack all
220 the un-done neighbors that we just processed */
221 if (!triDone[neighInfo[ii][0]]) {
222 triDone[neighInfo[ii][0]] = AIR_TRUE1;
223 okayIdx = airArrayLenIncr(okayArr, 1);
224 okay = AIR_CAST(unsigned int*, okayArr->data)((unsigned int*)(okayArr->data));
225 okay[okayIdx] = neighInfo[ii][0];
226 ++pushedNum;
227 }
228 } /* if (neighGot[ii]) */
229 } /* for ii */
230 return pushedNum;
231}
232
233/*
234** ONLY GOOD FOR limnPrimitiveTriangles!!
235*/
236static unsigned int
237maxTrianglePerPrimitive(limnPolyData *pld) {
238 unsigned int ret, primIdx;
239
240 ret = 0;
241 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
242 ret = AIR_MAX(ret, pld->icnt[primIdx]/3)((ret) > (pld->icnt[primIdx]/3) ? (ret) : (pld->icnt
[primIdx]/3));
243 }
244 return ret;
245}
246
247/*
248** fills nTriWithVert with 2D array about which triangles use which vertices
249*/
250static int
251triangleWithVertex(Nrrd *nTriWithVert, limnPolyData *pld) {
252 static const char me[]="triangleWithVertex";
253 unsigned int *triWithVertNum, /* vert ii has triWithVertNum[ii] tris */
254 *triWithVert, baseVertIdx, primIdx, vertIdx,
255 maxTriPerVert, totTriIdx;
256 airArray *mop;
257
258 if (!(nTriWithVert && pld)) {
259 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
260 return 1;
261 }
262 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
263 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
264 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
265 return 1;
266 }
267
268 triWithVertNum = AIR_CAST(unsigned int*,((unsigned int*)(calloc(pld->xyzwNum, sizeof(unsigned int)
)))
269 calloc(pld->xyzwNum, sizeof(unsigned int)))((unsigned int*)(calloc(pld->xyzwNum, sizeof(unsigned int)
)));
270 if (!triWithVertNum) {
271 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate temp array", me);
272 return 1;
273 }
274 mop = airMopNew();
275 airMopAdd(mop, triWithVertNum, airFree, airMopAlways);
276
277 /* fill in triWithVertNum */
278 baseVertIdx = 0;
279 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
280 unsigned int triNum, triIdx, *indxLine, ii;
281 triNum = pld->icnt[primIdx]/3;
282 for (triIdx=0; triIdx<triNum; triIdx++) {
283 indxLine = pld->indx + baseVertIdx + 3*triIdx;
284 for (ii=0; ii<3; ii++) {
285 triWithVertNum[indxLine[ii]]++;
286 }
287 }
288 baseVertIdx += pld->icnt[primIdx];
289 }
290
291 /* find max # tris per vert, allocate output */
292 maxTriPerVert = 0;
293 for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
294 maxTriPerVert = AIR_MAX(maxTriPerVert, triWithVertNum[vertIdx])((maxTriPerVert) > (triWithVertNum[vertIdx]) ? (maxTriPerVert
) : (triWithVertNum[vertIdx]));
295 }
296 if (nrrdMaybeAlloc_va(nTriWithVert, nrrdTypeUInt, 2,
297 AIR_CAST(size_t, 1 + maxTriPerVert)((size_t)(1 + maxTriPerVert)),
298 AIR_CAST(size_t, pld->xyzwNum)((size_t)(pld->xyzwNum)))) {
299 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
300 airMopError(mop); return 1;
301 }
302 triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data)((unsigned int*)(nTriWithVert->data));
303
304 baseVertIdx = 0;
305 totTriIdx = 0;
306 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
307 unsigned int triNum, *indxLine, *twvLine, ii, triIdx;
308 triNum = pld->icnt[primIdx]/3;
309 for (triIdx=0; triIdx<triNum; triIdx++) {
310 indxLine = pld->indx + baseVertIdx + 3*triIdx;
311 for (ii=0; ii<3; ii++) {
312 twvLine = triWithVert + (1+maxTriPerVert)*indxLine[ii];
313 twvLine[1+twvLine[0]] = totTriIdx;
314 twvLine[0]++;
315 }
316 ++totTriIdx;
317 }
318 baseVertIdx += pld->icnt[primIdx];
319 }
320
321 airMopOkay(mop);
322 return 0;
323}
324
325/*
326** learns which (three vertices) are with which triangle
327*/
328static int
329vertexWithTriangle(Nrrd *nVertWithTri, limnPolyData *pld) {
330 static const char me[]="vertexWithTriangle";
331 unsigned int baseVertIdx, primIdx, *vertWithTri, triNum;
332
333 if (!(nVertWithTri && pld)) {
334 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
335 return 1;
336 }
337 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
338 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
339 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
340 return 1;
341 }
342
343 triNum = limnPolyDataPolygonNumber(pld);
344 if (nrrdMaybeAlloc_va(nVertWithTri, nrrdTypeUInt, 2,
345 AIR_CAST(size_t, 3)((size_t)(3)),
346 AIR_CAST(size_t, triNum)((size_t)(triNum)))) {
347 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
348 return 1;
349 }
350 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
351
352 baseVertIdx = 0;
353 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
354 unsigned int triIdx, *indxLine, totTriIdx, ii;
355 triNum = pld->icnt[primIdx]/3;
356 for (triIdx=0; triIdx<triNum; triIdx++) {
357 totTriIdx = triIdx + baseVertIdx/3;
358 indxLine = pld->indx + baseVertIdx + 3*triIdx;
359 for (ii=0; ii<3; ii++) {
360 (vertWithTri + 3*totTriIdx)[ii] = indxLine[ii];
361 }
362 }
363 baseVertIdx += pld->icnt[primIdx];
364 }
365
366 return 0;
367}
368
369static int
370splitListExtract(unsigned int *listLenP,
371 airArray *edgeArr, unsigned char *hitCount,
372 unsigned int firstVertIdx, unsigned int edgeDoneNum) {
373 static const char me[]="splitListExtract";
374 unsigned int *edgeData, edgeNum, *edgeLine, edgeIdx, edgeTmp[5],
375 tmp, nextVertIdx, listLen;
376
377 edgeNum = edgeArr->len;
378 edgeData = AIR_CAST(unsigned int*, edgeArr->data)((unsigned int*)(edgeArr->data));
379 edgeNum -= edgeDoneNum;
380 edgeData += 5*edgeDoneNum;
381
382 /* put first edge in first position */
383 for (edgeIdx=0; edgeIdx<edgeNum; edgeIdx++) {
384 edgeLine = edgeData + 5*edgeIdx;
385 if (edgeLine[2] == firstVertIdx || edgeLine[3] == firstVertIdx) {
386 break;
387 }
388 }
389 if (edgeIdx == edgeNum) {
390 biffAddf(LIMNlimnBiffKey, "%s: never found first vertex %u", me, firstVertIdx);
391 return 1;
392 }
393 if (edgeLine[3] == firstVertIdx) {
394 ELL_SWAP2(edgeLine[2], edgeLine[3], tmp)((tmp)=(edgeLine[2]),(edgeLine[2])=(edgeLine[3]),(edgeLine[3]
)=(tmp));
395 }
396 ELL_5V_COPY(edgeTmp, edgeData)((edgeTmp)[0]=(edgeData)[0], (edgeTmp)[1]=(edgeData)[1], (edgeTmp
)[2]=(edgeData)[2], (edgeTmp)[3]=(edgeData)[3], (edgeTmp)[4]=
(edgeData)[4]);
397 ELL_5V_COPY(edgeData, edgeLine)((edgeData)[0]=(edgeLine)[0], (edgeData)[1]=(edgeLine)[1], (edgeData
)[2]=(edgeLine)[2], (edgeData)[3]=(edgeLine)[3], (edgeData)[4
]=(edgeLine)[4]);
398 ELL_5V_COPY(edgeLine, edgeTmp)((edgeLine)[0]=(edgeTmp)[0], (edgeLine)[1]=(edgeTmp)[1], (edgeLine
)[2]=(edgeTmp)[2], (edgeLine)[3]=(edgeTmp)[3], (edgeLine)[4]=
(edgeTmp)[4]);
399
400 /* start looking for the rest */
401 listLen = 1;
402 hitCount[firstVertIdx]--;
403 nextVertIdx = edgeData[3];
404 hitCount[nextVertIdx]--;
405 /*
406 fprintf(stderr, "!%s: found first %u --> %u (tris %u %u)\n", me,
407 firstVertIdx, nextVertIdx,
408 edgeData[0], edgeData[1]);
409 */
410
411 /* the search start progresses so that we don't see the same edge twice */
412#define SEARCH \
413 for (edgeIdx=listLen; edgeIdx<edgeNum; edgeIdx++) { \
414 edgeLine = edgeData + 5*edgeIdx; \
415 if (edgeLine[2] == nextVertIdx || edgeLine[3] == nextVertIdx) { \
416 break; \
417 } \
418 }
419 SEARCH;
420 while (edgeIdx < edgeNum) {
421 if (edgeLine[3] == nextVertIdx) {
422 ELL_SWAP2(edgeLine[2], edgeLine[3], tmp)((tmp)=(edgeLine[2]),(edgeLine[2])=(edgeLine[3]),(edgeLine[3]
)=(tmp));
423 }
424 ELL_5V_COPY(edgeTmp, edgeData + 5*listLen)((edgeTmp)[0]=(edgeData + 5*listLen)[0], (edgeTmp)[1]=(edgeData
+ 5*listLen)[1], (edgeTmp)[2]=(edgeData + 5*listLen)[2], (edgeTmp
)[3]=(edgeData + 5*listLen)[3], (edgeTmp)[4]=(edgeData + 5*listLen
)[4]);
425 ELL_5V_COPY(edgeData + 5*listLen, edgeLine)((edgeData + 5*listLen)[0]=(edgeLine)[0], (edgeData + 5*listLen
)[1]=(edgeLine)[1], (edgeData + 5*listLen)[2]=(edgeLine)[2], (
edgeData + 5*listLen)[3]=(edgeLine)[3], (edgeData + 5*listLen
)[4]=(edgeLine)[4]);
426 ELL_5V_COPY(edgeLine, edgeTmp)((edgeLine)[0]=(edgeTmp)[0], (edgeLine)[1]=(edgeTmp)[1], (edgeLine
)[2]=(edgeTmp)[2], (edgeLine)[3]=(edgeTmp)[3], (edgeLine)[4]=
(edgeTmp)[4]);
427 hitCount[nextVertIdx]--;
428 /*
429 fprintf(stderr, "!%s: (len %u) found %u --> %u (tris %u %u)\n", me,
430 listLen, nextVertIdx,
431 (edgeData + 5*listLen)[3],
432 (edgeData + 5*listLen)[0],
433 (edgeData + 5*listLen)[1]);
434 */
435 nextVertIdx = (edgeData + 5*listLen)[3];
436 hitCount[nextVertIdx]--;
437 listLen++;
438 SEARCH;
439 }
440 /*
441 fprintf(stderr, "!%s: finishing with Len %u, ended at %u\n", me,
442 listLen, nextVertIdx);
443 */
444 *listLenP = listLen;
445 return 0;
446#undef SEARCH
447}
448
449/*
450** returns the element of vert[] that is not v0 or v1
451*/
452static unsigned int
453sweepVertNext(unsigned int *vert, unsigned int v0, unsigned int v1) {
454 unsigned int v2;
455
456 v2 = vert[0];
457 if (v2 == v0 || v2 == v1) {
458 v2 = vert[1];
459 }
460 if (v2 == v0 || v2 == v1) {
461 v2 = vert[2];
462 }
463 return v2;
464}
465
466/*
467** returns non-zero iff A and B are in {v[0],v[1],v[2]}
468*/
469static int
470sweepHave2(unsigned int v[3], unsigned int A, unsigned B) {
471 int haveA, haveB;
472
473 haveA = (A == v[0] || A == v[1] || A == v[2]);
474 haveB = (B == v[0] || B == v[1] || B == v[2]);
475 return (haveA && haveB);
476}
477
478/*
479** returns UINT_MAX if there is no other triangle
480*/
481static unsigned int
482sweepTriNext(unsigned int *triLine, unsigned int v0, unsigned int v1,
483 unsigned int triNot, Nrrd *nVertWithTri) {
484 unsigned int triIdx, ret, *vertLine, *vertWithTri;
485
486 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
487
488 for (triIdx=0; triIdx<triLine[0]; triIdx++) {
489 if (triLine[1+triIdx] == triNot) {
490 continue;
491 }
492 vertLine = vertWithTri + 3*triLine[1+triIdx];
493 if (sweepHave2(vertLine, v0, v1)) {
494 break;
495 }
496 }
497 if (triIdx == triLine[0]) {
498 ret = UINT_MAX(2147483647 *2U +1U);
499 } else {
500 ret = triLine[1+triIdx];
501 }
502 return ret;
503}
504
505/*
506** the sweep does NOT include triStart, but it does include whichever
507** triStop it hit (if any)
508** returns: length of sweep
509** sweep: output (does not include triStart)
510** triStartIdx: what triangle to start at
511** vertPivotIdx, vertStartIdx: two vertices of start triangle; sweep
512** proceeds around the pivot index
513** triStop{0,1}Idx: triangles to stop sweeping at
514*/
515static unsigned int
516splitTriSweep(unsigned int *sweep,
517 unsigned int triStart,
518 unsigned int vertPivot, unsigned int vertStart,
519 unsigned int triStop0, unsigned int triStop1,
520 Nrrd *nTriWithVert, Nrrd *nVertWithTri) {
521 /* static const char me[]="splitTriSweep"; */
522 unsigned int sweepLen;
523 unsigned int maxTriPerVert, *triWithVert,
524 *vertWithTri, *triLine, *vertLine, triCurr, vertLast, vertNext;
525
526 maxTriPerVert = AIR_CAST(unsigned int, nTriWithVert->axis[0].size-1)((unsigned int)(nTriWithVert->axis[0].size-1));
527 triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data)((unsigned int*)(nTriWithVert->data));
528 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
529
530 /*
531 fprintf(stderr, "!%s: hi, triStart %u, pivot %u, start %u, "
532 "stop = %u, %u\n", me,
533 triStart, vertPivot, vertStart, triStop0, triStop1);
534 */
535 if (triStart == triStop0 || triStart == triStop1) {
536 /* nowhere to go */
537 return 0;
538 }
539
540 triLine = triWithVert + (1+maxTriPerVert)*vertPivot;
541 vertLast = vertStart;
542 triCurr = triStart;
543 sweepLen = 0;
544 do {
545 if (!(triCurr == triStart)) {
546 sweep[sweepLen++] = triCurr;
547 /*
548 fprintf(stderr, "!%s: saving sweep[%u] = %u\n", me,
549 sweepLen-1, triCurr);
550 */
551 }
552 vertLine = vertWithTri + 3*triCurr;
553 vertNext = sweepVertNext(vertLine, vertPivot, vertLast);
554 /*
555 fprintf(stderr, "!%s: vertNext(%u,%u) = %u\n", me,
556 vertPivot, vertLast, vertNext);
557 */
558 triCurr = sweepTriNext(triLine, vertPivot, vertNext,
559 triCurr, nVertWithTri);
560 /*
561 fprintf(stderr, "!%s: triNext(%u,%u) = %u\n", me,
562 vertPivot, vertNext, triCurr);
563 */
564 vertLast = vertNext;
565 } while (!( UINT_MAX(2147483647 *2U +1U) == triCurr
566 || triStart == triCurr
567 || triStop0 == triCurr
568 || triStop1 == triCurr ));
569 if (!( UINT_MAX(2147483647 *2U +1U) == triCurr )) {
570 sweep[sweepLen++] = triCurr;
571 /*
572 fprintf(stderr, "!%s: saving sweep[%u] = %u\n", me,
573 sweepLen-1, triCurr);
574 */
575 }
576
577 return sweepLen;
578}
579
580/*
581** track0: first triangle track, length *track0LenP
582** track1: first triangle track, length *track1LenP
583** sweep: buffer for sweep
584**
585** NOTE: triangles may be internally repeated in a track
586**
587** when vert path a loop on a non-orientable surface (e.g. mobius strip),
588** then track0 will NOT include the endpoint triangles
589** (or its not supposed to), and track1 will include them.
590*/
591static int
592splitTriTrack(unsigned int *track0, unsigned int *track0LenP,
593 unsigned int *track1, unsigned int *track1LenP,
594 unsigned int *sweep,
595 Nrrd *nTriWithVert, Nrrd *nVertWithTri,
596 airArray *edgeArr, unsigned startIdx, unsigned int listLen,
597 int looping) {
598 static const char me[]="splitTriTrack";
599 unsigned int len0, len1, *edgeData, *edgeLine, edgeIdx, triIdx,
600 /* maxTriPerVert, *triWithVert, *vertWithTri, */
601 sweepLen, loopEnd0, loopEnd1, loopStart0, loopStart1;
602 int doBack0, doBack1;
603
604 len0 = len1 = 0;
605 edgeData = AIR_CAST(unsigned int*, edgeArr->data)((unsigned int*)(edgeArr->data));
606 edgeData += 5*startIdx;
607 /* maxTriPerVert = AIR_CAST(unsigned int, nTriWithVert->axis[0].size-1); */
608 /* triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data); */
609 /* vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data); */
610
611 if (looping) {
612 loopStart0 = (edgeData)[0];
613 loopStart1 = (edgeData)[1];
614 loopEnd0 = (edgeData + 5*(listLen - 1))[0];
615 loopEnd1 = (edgeData + 5*(listLen - 1))[1];
616 /*
617 fprintf(stderr, "!%s: loop start = %u, %u, end = %u,%u\n", me,
618 loopStart0, loopStart1, loopEnd0, loopEnd1);
619 */
620 } else {
621 loopStart0 = loopStart1 = UINT_MAX(2147483647 *2U +1U);
622 loopEnd0 = loopEnd1 = UINT_MAX(2147483647 *2U +1U);
623 }
624
625 /* ,,,,,,,,,,,,,,,,,,,,,
626 fprintf(stderr, "!%s: 1st 2 tris %u %u, verts %u %u\n", me,
627 edgeData[0], edgeData[1], edgeData[2], edgeData[3]);
628 fprintf(stderr, "!%s: triangles at start vert %u:\n", me, edgeData[2]);
629 triLine = triWithVert + edgeData[2]*(1+maxTriPerVert);
630 for (triIdx=0; triIdx<triLine[0]; triIdx++) {
631 unsigned int *vertLine;
632 vertLine = vertWithTri + 3*triLine[1+triIdx];
633 fprintf(stderr, "!%s: %u: %u (verts %u %u %u)\n",
634 me, triIdx, triLine[1+triIdx],
635 vertLine[0], vertLine[1], vertLine[2]);
636 }
637 ````````````````````` */
638
639 /* we turn on backward sweeping for the initial edge;
640 doBack{0,1} will be set explicitly at each edge thereafter */
641 doBack0 = doBack1 = AIR_TRUE1;
642 for (edgeIdx=0; edgeIdx<(looping
643 ? listLen-1
644 : listLen); edgeIdx++) {
645 unsigned int stop0, stop1;
646 edgeLine = edgeData + 5*edgeIdx;
647 /* ,,,,,,,,,,,,,,,,,,,,,
648 fprintf(stderr, "!%s: edge %u: vert %u->%u, tris %u, %u\n", me,
649 edgeIdx, edgeLine[2], edgeLine[3],
650 edgeLine[0], edgeLine[1]);
651 fprintf(stderr, "!%s: triangles at next vert %u:\n", me, edgeLine[3]);
652 triLine = triWithVert + edgeLine[3]*(1+maxTriPerVert);
653 for (triIdx=0; triIdx<triLine[0]; triIdx++) {
654 vertLine = vertWithTri + 3*triLine[1+triIdx];
655 fprintf(stderr, "!%s: %u: %u (verts %u %u %u)\n",
656 me, triIdx, triLine[1+triIdx],
657 vertLine[0], vertLine[1], vertLine[2]);
658 }
659 ````````````````````` */
660 if (0 == edgeIdx && looping) {
661 /* sweeps from 1st link on loop are stopped by a tris on last edge */
662 stop0 = loopEnd0;
663 stop1 = loopEnd1;
664 } else {
665 stop0 = UINT_MAX(2147483647 *2U +1U);
666 stop1 = UINT_MAX(2147483647 *2U +1U);
667 }
668 if (doBack0) {
669 sweepLen = splitTriSweep(sweep, edgeLine[0], edgeLine[2], edgeLine[3],
670 stop0, stop1, nTriWithVert, nVertWithTri);
671 if (0 == edgeIdx && looping && sweepLen > 0) {
672 /* don't include either stop triangle on track 0 */
673 for (triIdx=0; triIdx<sweepLen-1; triIdx++) {
674 track0[len0++] = sweep[sweepLen-2-triIdx];
675 }
676 } else {
677 for (triIdx=0; triIdx<sweepLen; triIdx++) {
678 track0[len0++] = sweep[sweepLen-1-triIdx];
679 }
680 }
681 track0[len0++] = edgeLine[0];
682 }
683 if (doBack1) {
684 sweepLen = splitTriSweep(sweep, edgeLine[1], edgeLine[2], edgeLine[3],
685 stop0, stop1, nTriWithVert, nVertWithTri);
686 /* on this side we *do* include the stop triangle */
687 for (triIdx=0; triIdx<sweepLen; triIdx++) {
688 track1[len1++] = sweep[sweepLen-1-triIdx];
689 }
690 track1[len1++] = edgeLine[1];
691 }
692
693 if (edgeIdx<listLen-1) {
694 stop0 = (edgeLine + 5)[0];
695 stop1 = (edgeLine + 5)[1];
696 } else {
697 if (looping) {
698 stop0 = loopStart0;
699 stop1 = loopStart1;
700 } else {
701 stop0 = UINT_MAX(2147483647 *2U +1U);
702 stop1 = UINT_MAX(2147483647 *2U +1U);
703 }
704 }
705 sweepLen = splitTriSweep(sweep, edgeLine[0], edgeLine[3], edgeLine[2],
706 stop0, stop1, nTriWithVert, nVertWithTri);
707 for (triIdx=0; triIdx<sweepLen; triIdx++) {
708 track0[len0++] = sweep[triIdx];
709 }
710 sweepLen = splitTriSweep(sweep, edgeLine[1], edgeLine[3], edgeLine[2],
711 stop0, stop1, nTriWithVert, nVertWithTri);
712 for (triIdx=0; triIdx<sweepLen; triIdx++) {
713 track1[len1++] = sweep[triIdx];
714 }
715 if (edgeIdx<listLen-1) {
716 unsigned int *nextLine, tmp;
717 /* re-arrange the next edgeLine according to sweep results */
718 nextLine = edgeData + 5*(1 + edgeIdx);
719 if (track0[len0-1] == nextLine[0]
720 && track1[len1-1] == nextLine[1]) {
721 /* fprintf(stderr, "!%s: tracking went 0->0, 1->1\n", me); */
722 doBack0 = doBack1 = AIR_FALSE0;
723 } else if (track0[len0-1] == nextLine[1]
724 && track1[len1-1] == nextLine[0]) {
725 /* fprintf(stderr, "!%s: tracking went 0->1, 0->1\n", me); */
726 ELL_SWAP2(nextLine[0], nextLine[1], tmp)((tmp)=(nextLine[0]),(nextLine[0])=(nextLine[1]),(nextLine[1]
)=(tmp));
727 doBack0 = doBack1 = AIR_FALSE0;
728 } else if (track0[len0-1] == nextLine[0]) {
729 /* fprintf(stderr, "!%s: tracking went 0->0, 1->x\n", me); */
730 doBack0 = AIR_FALSE0;
731 doBack1 = AIR_TRUE1;
732 } else if (track1[len1-1] == nextLine[1]) {
733 /* fprintf(stderr, "!%s: tracking went 0->x, 1->1\n", me); */
734 doBack0 = AIR_TRUE1;
735 doBack1 = AIR_FALSE0;
736 } else if (track0[len0-1] == nextLine[1]) {
737 /* fprintf(stderr, "!%s: tracking went 0->1, 1->x\n", me); */
738 ELL_SWAP2(nextLine[0], nextLine[1], tmp)((tmp)=(nextLine[0]),(nextLine[0])=(nextLine[1]),(nextLine[1]
)=(tmp));
739 doBack0 = AIR_FALSE0;
740 doBack1 = AIR_TRUE1;
741 } else if (track1[len1-1] == nextLine[0]) {
742 /* fprintf(stderr, "!%s: tracking went 0->x, 1->0\n", me); */
743 ELL_SWAP2(nextLine[0], nextLine[1], tmp)((tmp)=(nextLine[0]),(nextLine[0])=(nextLine[1]),(nextLine[1]
)=(tmp));
744 doBack0 = AIR_TRUE1;
745 doBack1 = AIR_FALSE0;
746 } else {
747 biffAddf(LIMNlimnBiffKey, "%s: edge %u/%u, sweep ends %u,%u != want %u,%u", me,
748 edgeIdx, listLen, track0[len0-1], track1[len1-1],
749 nextLine[0], nextLine[1]);
750 return 1;
751 }
752 } else {
753 doBack0 = doBack1 = AIR_FALSE0;
754 }
755 }
756 if (looping) {
757 /* the end of track0 shouldn't include the stop */
758 len0--;
759 }
760
761 *track0LenP = len0;
762 *track1LenP = len1;
763 return 0;
764}
765
766static int
767splitVertDup(limnPolyData *pld, airArray *edgeArr,
768 unsigned int edgeDoneNum, unsigned int listLen,
769 unsigned int *track, unsigned int trackLen,
770 int looping) {
771 static const char me[]="splitVertDup";
772 unsigned int *vixLut, ii, vixLutLen, oldVertNum, newVertNum, *edgeData,
773 bitflag, trackIdx, vert0, vert1;
774 airArray *mop;
775 limnPolyData pldTmp;
776
777 mop = airMopNew();
778 edgeData = AIR_CAST(unsigned int*, edgeArr->data)((unsigned int*)(edgeArr->data));
779 edgeData += 5*edgeDoneNum;
780 oldVertNum = pld->xyzwNum;
781 vixLutLen = looping ? listLen : listLen+1;
782 newVertNum = oldVertNum + vixLutLen;
783
784 /* quiet compiler warnings */
785 pldTmp.rgba = NULL((void*)0);
786 pldTmp.norm = NULL((void*)0);
787 pldTmp.tex2 = NULL((void*)0);
788 pldTmp.tang = NULL((void*)0);
789
790 if (looping) {
791 vert0 = edgeData[2]; /* don't use dupe of this on first triangle */
792 vert1 = edgeData[3]; /* don't use dupe of this on last triangle */
793 } else {
794 vert0 = vert1 = UINT_MAX(2147483647 *2U +1U);
795 }
796
797 /* HEY: sneakily preserve the old per-vertex arrays; we own them now */
798 pldTmp.xyzw = pld->xyzw;
799 airMopAdd(mop, pldTmp.xyzw, airFree, airMopAlways);
800 pld->xyzw = NULL((void*)0);
801 pld->xyzwNum = 0;
802 bitflag = limnPolyDataInfoBitFlag(pld);
803 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
804 pldTmp.rgba = pld->rgba;
805 airMopAdd(mop, pldTmp.rgba, airFree, airMopAlways);
806 pld->rgba = NULL((void*)0);
807 pld->rgbaNum = 0;
808 }
809 if ((1 << limnPolyDataInfoNorm) & bitflag) {
810 pldTmp.norm = pld->norm;
811 airMopAdd(mop, pldTmp.norm, airFree, airMopAlways);
812 pld->norm = NULL((void*)0);
813 pld->normNum = 0;
814 }
815 if ((1 << limnPolyDataInfoTex2) & bitflag) {
816 pldTmp.tex2 = pld->tex2;
817 airMopAdd(mop, pldTmp.tex2, airFree, airMopAlways);
818 pld->tex2 = NULL((void*)0);
819 pld->tex2Num = 0;
820 }
821 if ((1 << limnPolyDataInfoTang) & bitflag) {
822 pldTmp.tang = pld->tang;
823 airMopAdd(mop, pldTmp.tang, airFree, airMopAlways);
824 pld->tang = NULL((void*)0);
825 pld->tangNum = 0;
826 }
827 if (limnPolyDataAlloc(pld, bitflag, newVertNum,
828 pld->indxNum, pld->primNum)) {
829 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate new vert # %u", me, newVertNum);
830 airMopError(mop); return 1;
831 }
832
833 /* copy old data */
834 memcpy(pld->xyzw, pldTmp.xyzw, oldVertNum*4*sizeof(float))__builtin___memcpy_chk (pld->xyzw, pldTmp.xyzw, oldVertNum
*4*sizeof(float), __builtin_object_size (pld->xyzw, 0));
835 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
836 memcpy(pld->rgba, pldTmp.rgba, oldVertNum*4*sizeof(unsigned char))__builtin___memcpy_chk (pld->rgba, pldTmp.rgba, oldVertNum
*4*sizeof(unsigned char), __builtin_object_size (pld->rgba
, 0));
837 }
838 if ((1 << limnPolyDataInfoNorm) & bitflag) {
839 memcpy(pld->norm, pldTmp.norm, oldVertNum*3*sizeof(float))__builtin___memcpy_chk (pld->norm, pldTmp.norm, oldVertNum
*3*sizeof(float), __builtin_object_size (pld->norm, 0));
840 }
841 if ((1 << limnPolyDataInfoTex2) & bitflag) {
842 memcpy(pld->tex2, pldTmp.tex2, oldVertNum*2*sizeof(float))__builtin___memcpy_chk (pld->tex2, pldTmp.tex2, oldVertNum
*2*sizeof(float), __builtin_object_size (pld->tex2, 0));
843 }
844 if ((1 << limnPolyDataInfoTang) & bitflag) {
845 memcpy(pld->tang, pldTmp.tang, oldVertNum*3*sizeof(float))__builtin___memcpy_chk (pld->tang, pldTmp.tang, oldVertNum
*3*sizeof(float), __builtin_object_size (pld->tang, 0));
846 }
847
848 vixLut = AIR_CAST(unsigned int *, calloc(2*vixLutLen,((unsigned int *)(calloc(2*vixLutLen, sizeof(unsigned int))))
849 sizeof(unsigned int)))((unsigned int *)(calloc(2*vixLutLen, sizeof(unsigned int))));
850 airMopAdd(mop, vixLut, airFree, airMopAlways);
851 if (looping) {
852 for (ii=0; ii<vixLutLen; ii++) {
853 vixLut[0 + 2*ii] = edgeData[2 + 5*ii];
854 vixLut[1 + 2*ii] = oldVertNum + ii;
855 }
856 } else {
857 for (ii=0; ii<vixLutLen-1; ii++) {
858 vixLut[0 + 2*ii] = edgeData[2 + 5*ii];
859 vixLut[1 + 2*ii] = oldVertNum + ii;
860 }
861 /* now ii == vixLutLen-1 == listLen */
862 vixLut[0 + 2*ii] = edgeData[3 + 5*(ii-1)];
863 vixLut[1 + 2*ii] = oldVertNum + ii;
864 }
865
866 /* copy pld's vertex information to duped vertices */
867 for (ii=0; ii<vixLutLen; ii++) {
868 ELL_4V_COPY(pld->xyzw + 4*vixLut[1 + 2*ii],((pld->xyzw + 4*vixLut[1 + 2*ii])[0] = (pld->xyzw + 4*vixLut
[0 + 2*ii])[0], (pld->xyzw + 4*vixLut[1 + 2*ii])[1] = (pld
->xyzw + 4*vixLut[0 + 2*ii])[1], (pld->xyzw + 4*vixLut[
1 + 2*ii])[2] = (pld->xyzw + 4*vixLut[0 + 2*ii])[2], (pld->
xyzw + 4*vixLut[1 + 2*ii])[3] = (pld->xyzw + 4*vixLut[0 + 2
*ii])[3])
869 pld->xyzw + 4*vixLut[0 + 2*ii])((pld->xyzw + 4*vixLut[1 + 2*ii])[0] = (pld->xyzw + 4*vixLut
[0 + 2*ii])[0], (pld->xyzw + 4*vixLut[1 + 2*ii])[1] = (pld
->xyzw + 4*vixLut[0 + 2*ii])[1], (pld->xyzw + 4*vixLut[
1 + 2*ii])[2] = (pld->xyzw + 4*vixLut[0 + 2*ii])[2], (pld->
xyzw + 4*vixLut[1 + 2*ii])[3] = (pld->xyzw + 4*vixLut[0 + 2
*ii])[3]);
870 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
871 ELL_4V_COPY(pld->rgba + 4*vixLut[1 + 2*ii],((pld->rgba + 4*vixLut[1 + 2*ii])[0] = (pld->rgba + 4*vixLut
[0 + 2*ii])[0], (pld->rgba + 4*vixLut[1 + 2*ii])[1] = (pld
->rgba + 4*vixLut[0 + 2*ii])[1], (pld->rgba + 4*vixLut[
1 + 2*ii])[2] = (pld->rgba + 4*vixLut[0 + 2*ii])[2], (pld->
rgba + 4*vixLut[1 + 2*ii])[3] = (pld->rgba + 4*vixLut[0 + 2
*ii])[3])
872 pld->rgba + 4*vixLut[0 + 2*ii])((pld->rgba + 4*vixLut[1 + 2*ii])[0] = (pld->rgba + 4*vixLut
[0 + 2*ii])[0], (pld->rgba + 4*vixLut[1 + 2*ii])[1] = (pld
->rgba + 4*vixLut[0 + 2*ii])[1], (pld->rgba + 4*vixLut[
1 + 2*ii])[2] = (pld->rgba + 4*vixLut[0 + 2*ii])[2], (pld->
rgba + 4*vixLut[1 + 2*ii])[3] = (pld->rgba + 4*vixLut[0 + 2
*ii])[3]);
873 }
874 if ((1 << limnPolyDataInfoNorm) & bitflag) {
875 ELL_3V_COPY(pld->norm + 3*vixLut[1 + 2*ii],((pld->norm + 3*vixLut[1 + 2*ii])[0] = (pld->norm + 3*vixLut
[0 + 2*ii])[0], (pld->norm + 3*vixLut[1 + 2*ii])[1] = (pld
->norm + 3*vixLut[0 + 2*ii])[1], (pld->norm + 3*vixLut[
1 + 2*ii])[2] = (pld->norm + 3*vixLut[0 + 2*ii])[2])
876 pld->norm + 3*vixLut[0 + 2*ii])((pld->norm + 3*vixLut[1 + 2*ii])[0] = (pld->norm + 3*vixLut
[0 + 2*ii])[0], (pld->norm + 3*vixLut[1 + 2*ii])[1] = (pld
->norm + 3*vixLut[0 + 2*ii])[1], (pld->norm + 3*vixLut[
1 + 2*ii])[2] = (pld->norm + 3*vixLut[0 + 2*ii])[2]);
877 }
878 if ((1 << limnPolyDataInfoTex2) & bitflag) {
879 ELL_2V_COPY(pld->tex2 + 2*vixLut[1 + 2*ii],((pld->tex2 + 2*vixLut[1 + 2*ii])[0] = (pld->tex2 + 2*vixLut
[0 + 2*ii])[0], (pld->tex2 + 2*vixLut[1 + 2*ii])[1] = (pld
->tex2 + 2*vixLut[0 + 2*ii])[1])
880 pld->tex2 + 2*vixLut[0 + 2*ii])((pld->tex2 + 2*vixLut[1 + 2*ii])[0] = (pld->tex2 + 2*vixLut
[0 + 2*ii])[0], (pld->tex2 + 2*vixLut[1 + 2*ii])[1] = (pld
->tex2 + 2*vixLut[0 + 2*ii])[1]);
881 }
882 if ((1 << limnPolyDataInfoTang) & bitflag) {
883 ELL_3V_COPY(pld->tang + 3*vixLut[1 + 2*ii],((pld->tang + 3*vixLut[1 + 2*ii])[0] = (pld->tang + 3*vixLut
[0 + 2*ii])[0], (pld->tang + 3*vixLut[1 + 2*ii])[1] = (pld
->tang + 3*vixLut[0 + 2*ii])[1], (pld->tang + 3*vixLut[
1 + 2*ii])[2] = (pld->tang + 3*vixLut[0 + 2*ii])[2])
884 pld->tang + 3*vixLut[0 + 2*ii])((pld->tang + 3*vixLut[1 + 2*ii])[0] = (pld->tang + 3*vixLut
[0 + 2*ii])[0], (pld->tang + 3*vixLut[1 + 2*ii])[1] = (pld
->tang + 3*vixLut[0 + 2*ii])[1], (pld->tang + 3*vixLut[
1 + 2*ii])[2] = (pld->tang + 3*vixLut[0 + 2*ii])[2]);
885 }
886 }
887
888 /* for triangles in track, update indices of duped vertices */
889 /* we do this updating ONLY in the limnPolyData, and that's okay:
890 the split information is computed entirely from nVertWithTri
891 and nTriWithVert (which were based on the original polydata),
892 but not the current polydata */
893 /* HEY: this is one place where we really exploit the fact that we only
894 have triangles: it makes the indxLine computation much much easier */
895 for (trackIdx=0; trackIdx<trackLen; trackIdx++) {
896 unsigned int *indxLine, jj;
897 indxLine = pld->indx + 3*track[trackIdx];
898 for (ii=0; ii<vixLutLen; ii++) {
899 for (jj=0; jj<3; jj++) {
900 if (indxLine[jj] == vixLut[0 + 2*ii]
901 && !((0 == trackIdx && indxLine[jj] == vert0)
902 || (trackLen-1 == trackIdx && indxLine[jj] == vert1))) {
903 indxLine[jj] = vixLut[1 + 2*ii];
904 }
905 }
906 }
907 }
908
909 airMopOkay(mop);
910 return 0;
911}
912
913/*
914** edge[0, 1]: two neighboring triangles,
915** edge[2, 3]: their shared vertices
916** edge[4]: non-zero if this split has been processed
917**
918** this really makes no effort to be fast (or comprehensible)
919**
920** HEY should we be returning some statistics (e.g. how many points added)?
921*/
922static int
923doSplitting(limnPolyData *pld, Nrrd *nTriWithVert, Nrrd *nVertWithTri,
924 airArray *edgeArr) {
925 static const char me[]="doSplitting";
926 unsigned int edgeIdx, *edgeData,
927 *edgeLine=NULL((void*)0), vertIdx, vertNum, splitNum, edgeDoneNum, listLen=0,
928 *track0, track0Len=0, *track1, *sweep, track1Len=0, maxTriPerVert;
929 unsigned char *hitCount;
930 airArray *mop;
931 int passIdx;
932
933 if (!edgeArr->len) {
934 /* actually, no splitting was required! */
935 return 0;
936 }
937
938 mop = airMopNew();
939 /* NOTE: It is necessary to save out the number of (initial)
940 number of vertices here, because as we do the splitting
941 (which is done once per track, as tracks are computed),
942 pld->xyzwNum will increase ... */
943 vertNum = pld->xyzwNum;
944 hitCount = AIR_CAST(unsigned char *, calloc(vertNum,((unsigned char *)(calloc(vertNum, sizeof(unsigned char))))
945 sizeof(unsigned char)))((unsigned char *)(calloc(vertNum, sizeof(unsigned char))));
946 maxTriPerVert = AIR_CAST(unsigned int, nTriWithVert->axis[0].size - 1)((unsigned int)(nTriWithVert->axis[0].size - 1));
947 track0 = AIR_CAST(unsigned int *, calloc(maxTriPerVert*edgeArr->len,((unsigned int *)(calloc(maxTriPerVert*edgeArr->len, sizeof
(unsigned int))))
948 sizeof(unsigned int)))((unsigned int *)(calloc(maxTriPerVert*edgeArr->len, sizeof
(unsigned int))));
949 track1 = AIR_CAST(unsigned int *, calloc(maxTriPerVert*edgeArr->len,((unsigned int *)(calloc(maxTriPerVert*edgeArr->len, sizeof
(unsigned int))))
950 sizeof(unsigned int)))((unsigned int *)(calloc(maxTriPerVert*edgeArr->len, sizeof
(unsigned int))));
951 sweep = AIR_CAST(unsigned int *, calloc(maxTriPerVert,((unsigned int *)(calloc(maxTriPerVert, sizeof(unsigned int))
))
952 sizeof(unsigned int)))((unsigned int *)(calloc(maxTriPerVert, sizeof(unsigned int))
));
953 if (!(hitCount && track0 && track1 && sweep)) {
954 biffAddf(LIMNlimnBiffKey, "%s: couldn't alloc buffers", me);
955 airMopError(mop); return 1;
956 }
957 airMopAdd(mop, hitCount, airFree, airMopAlways);
958 airMopAdd(mop, track0, airFree, airMopAlways);
959 airMopAdd(mop, track1, airFree, airMopAlways);
960 airMopAdd(mop, sweep, airFree, airMopAlways);
961
962 edgeData = AIR_CAST(unsigned int*, edgeArr->data)((unsigned int*)(edgeArr->data));
963
964 /* initialize hitCount */
965 for (edgeIdx=0; edgeIdx<edgeArr->len; edgeIdx++) {
966 unsigned int ha, hb;
967 edgeLine = edgeData + 5*edgeIdx;
968 ha = hitCount[edgeLine[2]]++;
969 hb = hitCount[edgeLine[3]]++;
970 if (ha > 2 || hb > 2) {
971 biffAddf(LIMNlimnBiffKey, "%s: edge %u (vert %u %u) created hit counts %u %u", me,
972 edgeIdx, edgeLine[2], edgeLine[3], ha, hb);
973 airMopError(mop); return 1;
974 }
975 }
976
977 /* scan hitCount */
978#define SEARCH(x) \
979 for (vertIdx=0; vertIdx<vertNum; vertIdx++) { \
980 if ((x) == hitCount[vertIdx]) { \
981 break; \
982 } \
983 }
984
985 splitNum = 0;
986 edgeDoneNum = 0;
987 /* pass 0: look for singleton hits ==> non-loop tracks
988 pass 1: look for hitCount[2] ==> loop tracks
989 */
990 for (passIdx=0; passIdx<2; passIdx++) {
991 if (0 == passIdx) {
992 SEARCH(1);
993 } else {
994 SEARCH(2);
995 }
996 while (vertIdx < vertNum) {
997 unsigned int E;
998 E = 0;
999 if (1) {
1000 unsigned int hitIdx, hitSum;
1001 hitSum = 0;
1002 for (hitIdx=0; hitIdx<vertNum; hitIdx++) {
1003 hitSum += hitCount[hitIdx];
1004 }
1005 /*
1006 fprintf(stderr, "!%s: PRE hitSum = %u (pass %u)\n", me,
1007 hitSum, passIdx);
1008 */
1009 }
1010 if (!E) E |= splitListExtract(&listLen, edgeArr, hitCount,
1011 vertIdx, edgeDoneNum);
1012 /* HEY: should do a splitListShorten() that cuts across repeated
1013 triangles, and then shifting downward the rest of the list.
1014 take care with loops. iterate until there is no shortening */
1015 /*
1016 if (1) {
1017 unsigned int hitIdx, hitSum;
1018 hitSum = 0;
1019 for (hitIdx=0; hitIdx<vertNum; hitIdx++) {
1020 hitSum += hitCount[hitIdx];
1021 }
1022 fprintf(stderr, "!%s: (%d) POST hitSum = %u (pass %u)\n", me, E,
1023 hitSum, passIdx);
1024 }
1025 if (1 == passIdx) {
1026 fprintf(stderr, "!%s: loop len %u, verts %u,%u --- %u,%u\n"
1027 " tris %u,%u --- %u,%u\n", me,
1028 listLen,
1029 (edgeData + 5*(edgeDoneNum + listLen - 1))[2],
1030 (edgeData + 5*(edgeDoneNum + listLen - 1))[3],
1031 (edgeData + 5*edgeDoneNum)[2],
1032 (edgeData + 5*edgeDoneNum)[3],
1033 (edgeData + 5*(edgeDoneNum + listLen - 1))[0],
1034 (edgeData + 5*(edgeDoneNum + listLen - 1))[1],
1035 (edgeData + 5*edgeDoneNum)[0],
1036 (edgeData + 5*edgeDoneNum)[1]);
1037 }
1038 */
1039 if (!E) E |= splitTriTrack(track0, &track0Len, track1, &track1Len,
1040 sweep, nTriWithVert, nVertWithTri,
1041 edgeArr, edgeDoneNum, listLen, passIdx);
1042 /* ,,,,,,,,,,,,,,,,,,,,,
1043 if (!E) {
1044 fprintf(stderr, "!%s: track0:\n", me);
1045 for (triIdx=0; triIdx<track0Len; triIdx++) {
1046 fprintf(stderr, "!%s: %u: %u\n", me, triIdx, track0[triIdx]);
1047 }
1048 fprintf(stderr, "!%s: track1:\n", me);
1049 for (triIdx=0; triIdx<track1Len; triIdx++) {
1050 fprintf(stderr, "!%s: %u: %u\n", me, triIdx, track1[triIdx]);
1051 }
1052 }
1053 ````````````````````` */
1054 /* see- this is the only time pld is used (so it can be modified) */
1055 /* HEY: we should be using track1, since that's the one that includes
1056 the endpoint triangles, but on a mobius strip demo it looked worse...
1057 this still needs debugging */
1058 if (!E) E |= splitVertDup(pld, edgeArr, edgeDoneNum, listLen,
1059 track0, track0Len, passIdx);
1060 if (E) {
1061 biffAddf(LIMNlimnBiffKey, "%s: trouble on split %u (done %u/%u)", me,
1062 splitNum, edgeDoneNum, AIR_CAST(unsigned int, edgeArr->len)((unsigned int)(edgeArr->len)));
1063 return 1;
1064 }
1065 edgeDoneNum += listLen;
1066 /*
1067 fprintf(stderr, "!%s: edgeDoneNum now %u (%u)\n", me,
1068 edgeDoneNum, AIR_CAST(unsigned int, edgeArr->len));
1069 */
1070 if (0 == passIdx) {
1071 SEARCH(1);
1072 } else {
1073 SEARCH(2);
1074 }
1075 }
1076 }
1077#undef SEARCH
1078 airMopOkay(mop);
1079 return 0;
1080}
1081
1082int
1083_limnPolyDataVertexWindingProcess(limnPolyData *pld, int splitting) {
1084 static const char me[]="limnPolyDataVertexWindingProcess";
1085 unsigned int
1086 primIdx, /* for indexing through primitives */
1087 triIdx, /* for indexing through triangles in each primitive */
1088 maxTriPerPrim, /* max # triangles per primitive, which is essential for
1089 the indexing of each triangle (in each primitive)
1090 into a single triangle index */
1091 totTriIdx, /* another triangle index */
1092 totTriNum, /* total # triangles */
1093 trueTriNum, /* correct total # triangles in all primitives */
1094 baseVertIdx, /* first vertex for current primitive */
1095 maxTriPerVert, /* max # of tris on single vertex */
1096 /* *triWithVert, 2D array ((1+maxTriPerVert) x pld->xyzwNum)
1097 of per-vertex triangles */
1098 *vertWithTri, /* 3D array (3 x maxTriPerPrim x pld->primNum)
1099 of per-tri vertices (vertex indices), which is just
1100 a repackaging of the information in the lpld */
1101 doneTriNum, /* # triangles finished so far */
1102 *intxBuff, /* stupid buffer */
1103 *okay, /* the stack of triangles with okay (possibly fixed)
1104 winding, but with some neighbors that may as yet
1105 need fixing */
1106 *split; /* stack of 5-tuples about edges needing vertex splits:
1107 split[0, 1]: two neighboring triangles,
1108 split[2, 3]: their shared vertices
1109 split[4]: non-zero if this split has been processed */
1110 unsigned char
1111 *triDone; /* 1D array (len totTriNum) record of done-ness */
1112 Nrrd *nTriWithVert, *nVertWithTri;
1113 airArray *mop, /* house-keeping */
1114 *okayArr, /* airArray around "okay" */
1115 *splitArr; /* airArray around "split" */
1116 airPtrPtrUnion appu;
1117 /*
1118 fprintf(stderr, "!%s: hi\n", me);
1119 */
1120 if (!pld) {
1121 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1122 return 1;
1123 }
1124
1125 if (!(pld->xyzwNum && pld->primNum)) {
1126 /* this is empty? */
1127 return 0;
1128 }
1129
1130 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
1131 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
1132 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
1133 return 1;
1134 }
1135
1136 maxTriPerPrim = maxTrianglePerPrimitive(pld);
1137 totTriNum = limnPolyDataPolygonNumber(pld);
1138
1139 mop = airMopNew();
1140 triDone = AIR_CAST(unsigned char *, calloc(totTriNum,((unsigned char *)(calloc(totTriNum, sizeof(unsigned char))))
1141 sizeof(unsigned char)))((unsigned char *)(calloc(totTriNum, sizeof(unsigned char))));
1142 airMopAdd(mop, triDone, airFree, airMopAlways);
1143 if (!triDone) {
1144 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate temp array", me);
1145 airMopError(mop); return 1;
1146 }
1147
1148 /* allocate TriWithVert, VertWithTri, intxBuff */
1149 nTriWithVert = nrrdNew();
1150 airMopAdd(mop, nTriWithVert, (airMopper)nrrdNuke, airMopAlways);
1151 nVertWithTri = nrrdNew();
1152 airMopAdd(mop, nVertWithTri, (airMopper)nrrdNuke, airMopAlways);
1153 if (triangleWithVertex(nTriWithVert, pld)
1154 || vertexWithTriangle(nVertWithTri, pld)) {
1155 biffAddf(LIMNlimnBiffKey, "%s: couldn't set nTriWithVert or nVertWithTri", me);
1156 airMopError(mop); return 1;
1157 }
1158 vertWithTri = AIR_CAST(unsigned int*, nVertWithTri->data)((unsigned int*)(nVertWithTri->data));
1159 /* triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data); */
1160
1161 maxTriPerVert = nTriWithVert->axis[0].size - 1;
1162 intxBuff = AIR_CAST(unsigned int*, calloc(maxTriPerVert,((unsigned int*)(calloc(maxTriPerVert, sizeof(unsigned int)))
)
1163 sizeof(unsigned int)))((unsigned int*)(calloc(maxTriPerVert, sizeof(unsigned int)))
);
1164 if (!intxBuff) {
1165 biffAddf(LIMNlimnBiffKey, "%s: failed to alloc an itty bitty buffer", me);
1166 airMopError(mop); return 1;
1167 }
1168 airMopAdd(mop, intxBuff, airFree, airMopAlways);
1169
1170 /*
1171 nrrdSave("triWithVert.nrrd", nTriWithVert, NULL);
1172 nrrdSave("vertWithTri.nrrd", nVertWithTri, NULL);
1173 */
1174
1175 /* create the stack of recently fixed triangles */
1176 appu.ui = &okay;
1177 okayArr = airArrayNew(appu.v, NULL((void*)0), sizeof(unsigned int),
1178 maxTriPerPrim);
1179 airMopAdd(mop, okayArr, (airMopper)airArrayNuke, airMopAlways);
1180 if (splitting) {
1181 appu.ui = &split;
1182 splitArr = airArrayNew(appu.v, NULL((void*)0), 5*sizeof(unsigned int),
1183 maxTriPerPrim);
1184 /* split set as it is used */
1185 } else {
1186 splitArr = NULL((void*)0);
1187 split = NULL((void*)0);
1188 }
1189
1190 /* the skinny */
1191 doneTriNum = 0;
1192 trueTriNum = 0;
1193 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
1194 trueTriNum += pld->icnt[primIdx]/3;
1195 }
1196 /*
1197 fprintf(stderr, "!%s: trueTriNum %u; other tri num %u\n", me,
1198 trueTriNum, limnPolyDataPolygonNumber(pld));
1199 */
1200 while (doneTriNum < trueTriNum) {
1201 /* find first undone triangle, which should be on a different
1202 connected component than any processed so far */
1203 for (totTriIdx=0; triDone[totTriIdx]; totTriIdx++)
1204 ;
1205 /* we use the winding of this triangle to determine the correct
1206 winding of all neighboring trianges, so this one is now done */
1207 triDone[totTriIdx] = AIR_TRUE1;
1208 ++doneTriNum;
1209 /*
1210 fprintf(stderr, "!%s: considering tri %u done (%u)\n",
1211 me, totTriIdx, doneTriNum);
1212 */
1213 doneTriNum += neighborsCheckPush(nTriWithVert, nVertWithTri,
1214 triDone, okayArr, intxBuff, splitArr,
1215 totTriIdx, splitting);
1216 while (okayArr->len) {
1217 unsigned int popped;
1218 popped = okay[okayArr->len-1];
1219 airArrayLenIncr(okayArr, -1);
1220 /*
1221 fprintf(stderr, "!%s: popped %u\n", me, popped);
1222 */
1223 doneTriNum += neighborsCheckPush(nTriWithVert, nVertWithTri,
1224 triDone, okayArr, intxBuff, splitArr,
1225 popped, splitting);
1226 }
1227 }
1228
1229 if (splitting) {
1230 if (doSplitting(pld, nTriWithVert, nVertWithTri, splitArr)) {
1231 biffAddf(LIMNlimnBiffKey, "%s: problem doing vertex splitting", me);
1232 return 1;
1233 }
1234 } else {
1235 /* Copy from nVertWithTri back into polydata */
1236 baseVertIdx = 0;
1237 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
1238 unsigned int triNum, *indxLine, ii;
1239 triNum = pld->icnt[primIdx]/3;
1240 for (triIdx=0; triIdx<triNum; triIdx++) {
1241 totTriIdx = triIdx + baseVertIdx/3;
1242 indxLine = pld->indx + baseVertIdx + 3*triIdx;
1243 for (ii=0; ii<3; ii++) {
1244 indxLine[ii] = (vertWithTri + 3*totTriIdx)[ii];
1245 }
1246 }
1247 baseVertIdx += pld->icnt[primIdx];
1248 }
1249 }
1250
1251 airMopOkay(mop);
1252 return 0;
1253}
1254
1255/*
1256** with non-zero splitting, this does vertex splitting so that
1257** non-orientable surfaces can be rendered without seams. Took longer
1258** to implement than intended.
1259**
1260** HEY: still has a bug in handling which triangles get which
1261** (new) vertices when the seam in the non-orientable surface
1262** is a closed loop. Can be debugged later...
1263*/
1264int
1265limnPolyDataVertexWindingFix(limnPolyData *pld, int splitting) {
1266 static const char me[]="limnPolyDataVertexWindingFix";
1267
1268 if (!splitting) {
1269 if (_limnPolyDataVertexWindingProcess(pld, AIR_FALSE0)) {
1270 biffAddf(LIMNlimnBiffKey, "%s: trouble", me);
1271 return 1;
1272 }
1273 } else {
1274 if (_limnPolyDataVertexWindingProcess(pld, AIR_FALSE0)
1275 || _limnPolyDataVertexWindingProcess(pld, AIR_TRUE1)) {
1276 biffAddf(LIMNlimnBiffKey, "%s: trouble", me);
1277 return 1;
1278 }
1279 }
1280 return 0;
1281}
1282
1283int
1284limnPolyDataCCFind(limnPolyData *pld) {
1285 static const char me[]="limnPolyDataCCFind";
1286 unsigned int realTriNum, *triMap, *triWithVert, vertIdx,
1287 *indxOld, *indxNew, primNumOld, *icntOld, *icntNew, *baseIndx,
1288 primIdxNew, primNumNew, passIdx, eqvNum=0;
1289 unsigned char *typeOld, *typeNew;
1290 Nrrd *nTriWithVert, *nccSize, *nTriMap;
1291 airArray *mop, *eqvArr;
1292
1293 if (!pld) {
1294 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1295 return 1;
1296 }
1297 if (!(pld->xyzwNum && pld->primNum)) {
1298 /* this is empty? */
1299 return 0;
1300 }
1301
1302 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
1303 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
1304 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
1305 return 1;
1306 }
1307
1308 mop = airMopNew();
1309
1310 realTriNum = limnPolyDataPolygonNumber(pld);
1311
1312 eqvArr = airArrayNew(NULL((void*)0), NULL((void*)0), 2*sizeof(unsigned int),
1313 /* this is only a heuristic */ pld->xyzwNum);
1314 airMopAdd(mop, eqvArr, (airMopper)airArrayNuke, airMopAlways);
1315
1316 nTriWithVert = nrrdNew();
1317 airMopAdd(mop, nTriWithVert, (airMopper)nrrdNuke, airMopAlways);
1318 if (triangleWithVertex(nTriWithVert, pld)) {
1319 biffAddf(LIMNlimnBiffKey, "%s: couldn't set nTriWithVert", me);
1320 airMopError(mop); return 1;
1321 }
1322
1323 /* simple profiling showed that stupid amount of time was spent
1324 adding the equivalences. So we go in two passes- first two see
1325 how many equivalences are needed, and then actually adding them */
1326 /* yea, so, its like you don't really even need an airArray ... */
1327 triWithVert = AIR_CAST(unsigned int*, nTriWithVert->data)((unsigned int*)(nTriWithVert->data));
1328 for (passIdx=0; passIdx<2; passIdx++) {
1329 if (0 == passIdx) {
1330 eqvNum = 0;
1331 } else {
1332 airArrayLenPreSet(eqvArr, eqvNum);
1333 }
1334 for (vertIdx=0; vertIdx<nTriWithVert->axis[1].size; vertIdx++) {
1335 unsigned int *triLine, triIdx;
1336 triLine = triWithVert + vertIdx*(nTriWithVert->axis[0].size);
1337 for (triIdx=1; triIdx<triLine[0]; triIdx++) {
1338 if (0 == passIdx) {
1339 ++eqvNum;
1340 } else {
1341 airEqvAdd(eqvArr, triLine[1], triLine[1+triIdx]);
1342 }
1343 }
1344 }
1345 }
1346
1347 nTriMap = nrrdNew();
1348 airMopAdd(mop, nTriMap, (airMopper)nrrdNuke, airMopAlways);
1349 nccSize = nrrdNew();
1350 airMopAdd(mop, nccSize, (airMopper)nrrdNuke, airMopAlways);
1351 if (nrrdMaybeAlloc_va(nTriMap, nrrdTypeUInt, 1,
1352 AIR_CAST(size_t, realTriNum)((size_t)(realTriNum)))) {
1353 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate equivalence map", me);
1354 airMopError(mop); return 1;
1355 }
1356 triMap = AIR_CAST(unsigned int*, nTriMap->data)((unsigned int*)(nTriMap->data));
1357 primNumNew = airEqvMap(eqvArr, triMap, realTriNum);
1358 if (nrrdHisto(nccSize, nTriMap, NULL((void*)0), NULL((void*)0), primNumNew, nrrdTypeUInt)) {
1359 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: couldn't histogram CC map", me);
1360 airMopError(mop); return 1;
1361 }
1362
1363 /* indxNumOld == indxNumNew */
1364 indxOld = pld->indx;
1365 primNumOld = pld->primNum;
1366 if (1 != primNumOld) {
1367 biffAddf(LIMNlimnBiffKey, "%s: sorry! stupid implementation can't "
1368 "do primNum %u (only 1)",
1369 me, primNumOld);
1370 airMopError(mop); return 1;
1371 }
1372 typeOld = pld->type;
1373 icntOld = pld->icnt;
1374 indxNew = AIR_CAST(unsigned int*,((unsigned int*)(calloc(pld->indxNum, sizeof(unsigned int)
)))
1375 calloc(pld->indxNum, sizeof(unsigned int)))((unsigned int*)(calloc(pld->indxNum, sizeof(unsigned int)
)));
1376 typeNew = AIR_CAST(unsigned char*,((unsigned char*)(calloc(primNumNew, sizeof(unsigned char))))
1377 calloc(primNumNew, sizeof(unsigned char)))((unsigned char*)(calloc(primNumNew, sizeof(unsigned char))));
1378 icntNew = AIR_CAST(unsigned int*,((unsigned int*)(calloc(primNumNew, sizeof(unsigned int))))
1379 calloc(primNumNew, sizeof(unsigned int)))((unsigned int*)(calloc(primNumNew, sizeof(unsigned int))));
1380 if (!(indxNew && typeNew && icntNew)) {
1381 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate new polydata arrays", me);
1382 airMopError(mop); return 1;
1383 }
1384 pld->indx = indxNew;
1385 pld->primNum = primNumNew;
1386 pld->type = typeNew;
1387 pld->icnt = icntNew;
1388 airMopAdd(mop, indxOld, airFree, airMopAlways);
1389 airMopAdd(mop, typeOld, airFree, airMopAlways);
1390 airMopAdd(mop, icntOld, airFree, airMopAlways);
1391
1392 /* this multi-pass thing is really stupid
1393 (and assumes stupid primNumOld = 1) */
1394 baseIndx = pld->indx;
1395 for (primIdxNew=0; primIdxNew<pld->primNum; primIdxNew++) {
1396 unsigned int realTriIdx;
1397 pld->type[primIdxNew] = limnPrimitiveTriangles;
1398 pld->icnt[primIdxNew] = 0;
1399 for (realTriIdx=0; realTriIdx<realTriNum; realTriIdx++) {
1400 if (triMap[realTriIdx] == primIdxNew) {
1401 ELL_3V_COPY(baseIndx, indxOld + 3*realTriIdx)((baseIndx)[0] = (indxOld + 3*realTriIdx)[0], (baseIndx)[1] =
(indxOld + 3*realTriIdx)[1], (baseIndx)[2] = (indxOld + 3*realTriIdx
)[2]);
1402 baseIndx += 3;
1403 pld->icnt[primIdxNew] += 3;
1404 }
1405 }
1406 }
1407
1408 airMopOkay(mop);
1409 return 0;
1410}
1411
1412int
1413limnPolyDataPrimitiveSort(limnPolyData *pld, const Nrrd *_nval) {
1414 static const char me[]="limnPolyDataPrimitiveSort";
1415 Nrrd *nval, *nrec;
1416 const Nrrd *ntwo[2];
1417 airArray *mop;
1418 double *rec;
1419 unsigned int primIdx, **startIndx, *indxNew, *baseIndx, *icntNew;
1420 unsigned char *typeNew;
1421 int E;
1422
1423 if (!(pld && _nval)) {
1
Taking false branch
→
1424 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1425 return 1;
1426 }
1427 if (!(1 == _nval->dim
2
←
Taking false branch
→
1428 && nrrdTypeBlock != _nval->type
1429 && _nval->axis[0].size == pld->primNum)) {
1430 biffAddf(LIMNlimnBiffKey, "%s: need 1-D %u-len scalar nrrd "
1431 "(not %u-D type %s, axis[0].size %u)", me,
1432 pld->primNum,
1433 _nval->dim, airEnumStr(nrrdType, _nval->type),
1434 AIR_CAST(unsigned int, _nval->axis[0].size)((unsigned int)(_nval->axis[0].size)));
1435 return 1;
1436 }
1437
1438 mop = airMopNew();
1439 nval = nrrdNew();
1440 airMopAdd(mop, nval, (airMopper)nrrdNuke, airMopAlways);
1441 nrec = nrrdNew();
1442 airMopAdd(mop, nrec, (airMopper)nrrdNuke, airMopAlways);
1443 E = 0;
1444 if (!E) E |= nrrdConvert(nval, _nval, nrrdTypeDouble);
3
←
Taking true branch
→
1445 ntwo[0] = nval;
1446 ntwo[1] = nval;
1447 if (!E) E |= nrrdJoin(nrec, ntwo, 2, 0, AIR_TRUE1);
4
←
Assuming 'E' is 0
→
5
←
Taking true branch
→
1448 if (E) {
6
←
Assuming 'E' is 0
→
7
←
Taking false branch
→
1449 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: problem creating records", me);
1450 airMopError(mop); return 1;
1451 }
1452 rec = AIR_CAST(double *, nrec->data)((double *)(nrec->data));
1453 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
8
←
Loop condition is true. Entering loop body
→
9
←
Loop condition is false. Execution continues on line 1456
→
1454 rec[1 + 2*primIdx] = primIdx;
1455 }
1456 qsort(rec, pld->primNum, 2*sizeof(double),
1457 nrrdValCompareInv[nrrdTypeDouble]);
1458
1459 startIndx = AIR_CAST(unsigned int**, calloc(pld->primNum,((unsigned int**)(calloc(pld->primNum, sizeof(unsigned int
*))))
1460 sizeof(unsigned int*)))((unsigned int**)(calloc(pld->primNum, sizeof(unsigned int
*))));
1461 indxNew = AIR_CAST(unsigned int*, calloc(pld->indxNum,((unsigned int*)(calloc(pld->indxNum, sizeof(unsigned int)
)))
1462 sizeof(unsigned int)))((unsigned int*)(calloc(pld->indxNum, sizeof(unsigned int)
)));
1463 icntNew = AIR_CAST(unsigned int*, calloc(pld->primNum,((unsigned int*)(calloc(pld->primNum, sizeof(unsigned int)
)))
10
←
Within the expansion of the macro 'AIR_CAST':
→
a
Memory is allocated
1464 sizeof(unsigned int)))((unsigned int*)(calloc(pld->primNum, sizeof(unsigned int)
)));
1465 typeNew = AIR_CAST(unsigned char*, calloc(pld->primNum,((unsigned char*)(calloc(pld->primNum, sizeof(unsigned char
))))
1466 sizeof(unsigned char)))((unsigned char*)(calloc(pld->primNum, sizeof(unsigned char
))));
1467 if (!(startIndx && indxNew && icntNew && typeNew)) {
11
←
Taking true branch
→
1468 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocated temp buffers", me);
12
←
Potential leak of memory pointed to by 'icntNew'
1469 airMopError(mop); return 1;
1470 }
1471 airMopAdd(mop, startIndx, airFree, airMopAlways);
1472
1473 baseIndx = pld->indx;
1474 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
1475 startIndx[primIdx] = baseIndx;
1476 baseIndx += pld->icnt[primIdx];
1477 }
1478 baseIndx = indxNew;
1479 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
1480 unsigned int sortIdx;
1481 sortIdx = AIR_CAST(unsigned int, rec[1 + 2*primIdx])((unsigned int)(rec[1 + 2*primIdx]));
1482 memcpy(baseIndx, startIndx[sortIdx],__builtin___memcpy_chk (baseIndx, startIndx[sortIdx], pld->
icnt[sortIdx]*sizeof(unsigned int), __builtin_object_size (baseIndx
, 0))
1483 pld->icnt[sortIdx]*sizeof(unsigned int))__builtin___memcpy_chk (baseIndx, startIndx[sortIdx], pld->
icnt[sortIdx]*sizeof(unsigned int), __builtin_object_size (baseIndx
, 0));
1484 icntNew[primIdx] = pld->icnt[sortIdx];
1485 typeNew[primIdx] = pld->type[sortIdx];
1486 baseIndx += pld->icnt[sortIdx];
1487 }
1488
1489 airFree(pld->indx);
1490 pld->indx = indxNew;
1491 airFree(pld->type);
1492 pld->type = typeNew;
1493 airFree(pld->icnt);
1494 pld->icnt = icntNew;
1495
1496 airMopOkay(mop);
1497 return 0;
1498}
1499
1500int
1501limnPolyDataVertexWindingFlip(limnPolyData *pld) {
1502 static const char me[]="limnPolyDataVertexWindingFlip";
1503 unsigned int baseVertIdx, primIdx;
1504
1505 if (!pld) {
1506 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1507 return 1;
1508 }
1509 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
1510 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
1511 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
1512 return 1;
1513 }
1514
1515 baseVertIdx = 0;
1516 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
1517 unsigned int triNum, triIdx, *indxLine, tmpIdx;
1518 triNum = pld->icnt[primIdx]/3;
1519 for (triIdx=0; triIdx<triNum; triIdx++) {
1520 indxLine = pld->indx + baseVertIdx + 3*triIdx;
1521 ELL_SWAP2(indxLine[0], indxLine[2], tmpIdx)((tmpIdx)=(indxLine[0]),(indxLine[0])=(indxLine[2]),(indxLine
[2])=(tmpIdx));
1522 }
1523 baseVertIdx += pld->icnt[primIdx];
1524 }
1525
1526 return 0;
1527}
1528
1529int
1530limnPolyDataPrimitiveSelect(limnPolyData *pldOut,
1531 const limnPolyData *pldIn,
1532 const Nrrd *_nmask) {
1533 static const char me[]="limnPolyDataPrimitiveSelect";
1534 Nrrd *nmask;
1535 double *mask;
1536 unsigned int oldBaseVertIdx, oldPrimIdx, oldVertIdx, bitflag,
1537 *old2newMap, *new2oldMap,
1538 newPrimNum, newBaseVertIdx, newPrimIdx, newIndxNum, newVertIdx, newVertNum;
1539 unsigned char *vertUsed;
1540 airArray *mop;
1541
1542 if (!(pldOut && pldIn && _nmask)) {
1543 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1544 return 1;
1545 }
1546 if (!(1 == _nmask->dim
1547 && nrrdTypeBlock != _nmask->type
1548 && _nmask->axis[0].size == pldIn->primNum)) {
1549 biffAddf(LIMNlimnBiffKey, "%s: need 1-D %u-len scalar nrrd "
1550 "(not %u-D type %s, axis[0].size %u)", me,
1551 pldIn->primNum, _nmask->dim, airEnumStr(nrrdType, _nmask->type),
1552 AIR_CAST(unsigned int, _nmask->axis[0].size)((unsigned int)(_nmask->axis[0].size)));
1553 return 1;
1554 }
1555
1556 mop = airMopNew();
1557 nmask = nrrdNew();
1558 airMopAdd(mop, nmask, (airMopper)nrrdNuke, airMopAlways);
1559 if (nrrdConvert(nmask, _nmask, nrrdTypeDouble)) {
1560 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: trouble converting mask to %s", me,
1561 airEnumStr(nrrdType, nrrdTypeDouble));
1562 return 1;
1563 }
1564 mask = AIR_CAST(double *, nmask->data)((double *)(nmask->data));
1565
1566 old2newMap = AIR_CAST(unsigned int *, calloc(pldIn->xyzwNum,((unsigned int *)(calloc(pldIn->xyzwNum, sizeof(unsigned int
))))
1567 sizeof(unsigned int)))((unsigned int *)(calloc(pldIn->xyzwNum, sizeof(unsigned int
))));
1568 airMopAdd(mop, old2newMap, airFree, airMopAlways);
1569 vertUsed = AIR_CAST(unsigned char *, calloc(pldIn->xyzwNum,((unsigned char *)(calloc(pldIn->xyzwNum, sizeof(unsigned char
))))
1570 sizeof(unsigned char)))((unsigned char *)(calloc(pldIn->xyzwNum, sizeof(unsigned char
))));
1571 airMopAdd(mop, vertUsed, airFree, airMopAlways);
1572
1573 /* initialize all verts as unused */
1574 for (oldVertIdx=0; oldVertIdx<pldIn->xyzwNum; oldVertIdx++) {
1575 vertUsed[oldVertIdx] = AIR_FALSE0;
1576 }
1577 /* mark the used verts, and count # new indices and primitives */
1578 oldBaseVertIdx = 0;
1579 newPrimNum = 0;
1580 newIndxNum = 0;
1581 for (oldPrimIdx=0; oldPrimIdx<pldIn->primNum; oldPrimIdx++) {
1582 unsigned indxIdx;
1583 if (mask[oldPrimIdx]) {
1584 for (indxIdx=0; indxIdx<pldIn->icnt[oldPrimIdx]; indxIdx++) {
1585 vertUsed[(pldIn->indx + oldBaseVertIdx)[indxIdx]] = AIR_TRUE1;
1586 }
1587 newIndxNum += pldIn->icnt[oldPrimIdx];
1588 newPrimNum++;
1589 }
1590 oldBaseVertIdx += pldIn->icnt[oldPrimIdx];
1591 }
1592 /* count the used verts, and set up map from old to new indices */
1593 newVertNum = 0;
1594 for (oldVertIdx=0; oldVertIdx<pldIn->xyzwNum; oldVertIdx++) {
1595 if (vertUsed[oldVertIdx]) {
1596 old2newMap[oldVertIdx] = newVertNum++;
1597 }
1598 }
1599 /* allocate and fill reverse map */
1600 new2oldMap = AIR_CAST(unsigned int *, calloc(newVertNum,((unsigned int *)(calloc(newVertNum, sizeof(unsigned int))))
1601 sizeof(unsigned int)))((unsigned int *)(calloc(newVertNum, sizeof(unsigned int))));
1602 airMopAdd(mop, new2oldMap, airFree, airMopAlways);
1603 newVertIdx = 0;
1604 for (oldVertIdx=0; oldVertIdx<pldIn->xyzwNum; oldVertIdx++) {
1605 if (vertUsed[oldVertIdx]) {
1606 new2oldMap[newVertIdx++] = oldVertIdx;
1607 }
1608 }
1609
1610 /* allocate output polydata */
1611 bitflag = limnPolyDataInfoBitFlag(pldIn);
1612 if (limnPolyDataAlloc(pldOut, bitflag, newVertNum, newIndxNum, newPrimNum)) {
1613 biffAddf(LIMNlimnBiffKey, "%s: trouble allocating output", me);
1614 return 1;
1615 }
1616
1617 /* transfer per-primitive information from old to new */
1618 oldBaseVertIdx = 0;
1619 newBaseVertIdx = 0;
1620 newPrimIdx = 0;
1621 for (oldPrimIdx=0; oldPrimIdx<pldIn->primNum; oldPrimIdx++) {
1622 if (mask[oldPrimIdx]) {
1623 unsigned indxIdx;
1624 pldOut->icnt[newPrimIdx] = pldIn->icnt[oldPrimIdx];
1625 pldOut->type[newPrimIdx] = pldIn->type[oldPrimIdx];
1626 for (indxIdx=0; indxIdx<pldIn->icnt[oldPrimIdx]; indxIdx++) {
1627 oldVertIdx = (pldIn->indx + oldBaseVertIdx)[indxIdx];
1628 (pldOut->indx + newBaseVertIdx)[indxIdx] = old2newMap[oldVertIdx];
1629 }
1630 newBaseVertIdx += pldIn->icnt[oldPrimIdx];
1631 newPrimIdx++;
1632 }
1633 oldBaseVertIdx += pldIn->icnt[oldPrimIdx];
1634 }
1635 /* transfer per-vertex info */
1636 for (newVertIdx=0; newVertIdx<newVertNum; newVertIdx++) {
1637 oldVertIdx = new2oldMap[newVertIdx];
1638 ELL_4V_COPY(pldOut->xyzw + 4*newVertIdx, pldIn->xyzw + 4*oldVertIdx)((pldOut->xyzw + 4*newVertIdx)[0] = (pldIn->xyzw + 4*oldVertIdx
)[0], (pldOut->xyzw + 4*newVertIdx)[1] = (pldIn->xyzw +
4*oldVertIdx)[1], (pldOut->xyzw + 4*newVertIdx)[2] = (pldIn
->xyzw + 4*oldVertIdx)[2], (pldOut->xyzw + 4*newVertIdx
)[3] = (pldIn->xyzw + 4*oldVertIdx)[3]);
1639 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1640 ELL_4V_COPY(pldOut->rgba + 4*newVertIdx, pldIn->rgba + 4*oldVertIdx)((pldOut->rgba + 4*newVertIdx)[0] = (pldIn->rgba + 4*oldVertIdx
)[0], (pldOut->rgba + 4*newVertIdx)[1] = (pldIn->rgba +
4*oldVertIdx)[1], (pldOut->rgba + 4*newVertIdx)[2] = (pldIn
->rgba + 4*oldVertIdx)[2], (pldOut->rgba + 4*newVertIdx
)[3] = (pldIn->rgba + 4*oldVertIdx)[3]);
1641 }
1642 if ((1 << limnPolyDataInfoNorm) & bitflag) {
1643 ELL_3V_COPY(pldOut->norm + 3*newVertIdx, pldIn->norm + 3*oldVertIdx)((pldOut->norm + 3*newVertIdx)[0] = (pldIn->norm + 3*oldVertIdx
)[0], (pldOut->norm + 3*newVertIdx)[1] = (pldIn->norm +
3*oldVertIdx)[1], (pldOut->norm + 3*newVertIdx)[2] = (pldIn
->norm + 3*oldVertIdx)[2]);
1644 }
1645 if ((1 << limnPolyDataInfoTex2) & bitflag) {
1646 ELL_3V_COPY(pldOut->tex2 + 2*newVertIdx, pldIn->tex2 + 2*oldVertIdx)((pldOut->tex2 + 2*newVertIdx)[0] = (pldIn->tex2 + 2*oldVertIdx
)[0], (pldOut->tex2 + 2*newVertIdx)[1] = (pldIn->tex2 +
2*oldVertIdx)[1], (pldOut->tex2 + 2*newVertIdx)[2] = (pldIn
->tex2 + 2*oldVertIdx)[2]);
1647 }
1648 if ((1 << limnPolyDataInfoTang) & bitflag) {
1649 ELL_3V_COPY(pldOut->tang + 3*newVertIdx, pldIn->tang + 3*oldVertIdx)((pldOut->tang + 3*newVertIdx)[0] = (pldIn->tang + 3*oldVertIdx
)[0], (pldOut->tang + 3*newVertIdx)[1] = (pldIn->tang +
3*oldVertIdx)[1], (pldOut->tang + 3*newVertIdx)[2] = (pldIn
->tang + 3*oldVertIdx)[2]);
1650 }
1651 }
1652
1653 airMopOkay(mop);
1654 return 0;
1655}
1656
1657/* Helper function for limnPolyDataClipMulti - clips the edge between
1658 * disc and kept that partially fulfills the thresholds and maintains
1659 * a data structure that keeps track of edges we have clipped already,
1660 * to avoid creating duplicate vertices.
1661 */
1662static int
1663clipEdge(int disc, int kept, Nrrd *nval, double *thresh, int *newIdx,
1664 airArray *llistArr, limnPolyData *pld, unsigned int bitflag,
1665 limnPolyData *newpld, airArray *xyzwArr, airArray *rgbaArr,
1666 airArray *normArr, airArray *tex2Arr, airArray *tangArr) {
1667 int ref=-1, *llist=(int*)llistArr->data;
1668 int next=newIdx[disc];
1669 double alpha=0;
1670 unsigned int i,q,p,nk;
1671 double (*lup)(const void *v, size_t I)=nrrdDLookup[nval->type];
1672 /* check if we clipped the edge previously */
1673 while (next!=-1) {
1674 if (llist[next]==kept) /* found the desired vertex */
1675 return llist[next+1];
1676 ref=next+2;
1677 next=llist[next+2];
1678 }
1679 /* we need to interpolate - find the weight */
1680 nk=(nval->dim==1)?1:nval->axis[0].size;
1681 for (i=0; i<nk; i++) {
1682 double discval = lup(nval->data, nk*disc+i);
1683 double keptval = lup(nval->data, nk*kept+i);
1684 double thisalpha = AIR_AFFINE(discval,thresh[i],keptval,0.0,1.0)( ((double)(1.0)-(0.0))*((double)(thresh[i])-(discval)) / ((double
)(keptval)-(discval)) + (0.0));
1685 if (thisalpha<1.0 && thisalpha>alpha)
1686 alpha=thisalpha;
1687 }
1688 /* add interpolated vertex */
1689 q=airArrayLenIncr(xyzwArr, 1);
1690 ELL_4V_LERP_TT(newpld->xyzw+4*q, float, alpha, pld->xyzw+4*disc, pld->xyzw+4*kept)((newpld->xyzw+4*q)[0] = ((float)((((alpha))*(((pld->xyzw
+4*kept)[0]) - ((pld->xyzw+4*disc)[0])) + ((pld->xyzw+4
*disc)[0])))), (newpld->xyzw+4*q)[1] = ((float)((((alpha))
*(((pld->xyzw+4*kept)[1]) - ((pld->xyzw+4*disc)[1])) + (
(pld->xyzw+4*disc)[1])))), (newpld->xyzw+4*q)[2] = ((float
)((((alpha))*(((pld->xyzw+4*kept)[2]) - ((pld->xyzw+4*disc
)[2])) + ((pld->xyzw+4*disc)[2])))), (newpld->xyzw+4*q)
[3] = ((float)((((alpha))*(((pld->xyzw+4*kept)[3]) - ((pld
->xyzw+4*disc)[3])) + ((pld->xyzw+4*disc)[3])))));
1691 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1692 airArrayLenIncr(rgbaArr, 1);
1693 ELL_4V_LERP_TT(newpld->rgba+4*q, unsigned char, alpha, pld->rgba+4*disc, pld->rgba+4*kept)((newpld->rgba+4*q)[0] = ((unsigned char)((((alpha))*(((pld
->rgba+4*kept)[0]) - ((pld->rgba+4*disc)[0])) + ((pld->
rgba+4*disc)[0])))), (newpld->rgba+4*q)[1] = ((unsigned char
)((((alpha))*(((pld->rgba+4*kept)[1]) - ((pld->rgba+4*disc
)[1])) + ((pld->rgba+4*disc)[1])))), (newpld->rgba+4*q)
[2] = ((unsigned char)((((alpha))*(((pld->rgba+4*kept)[2])
- ((pld->rgba+4*disc)[2])) + ((pld->rgba+4*disc)[2])))
), (newpld->rgba+4*q)[3] = ((unsigned char)((((alpha))*(((
pld->rgba+4*kept)[3]) - ((pld->rgba+4*disc)[3])) + ((pld
->rgba+4*disc)[3])))));
1694 }
1695 if ((1 << limnPolyDataInfoNorm) & bitflag) {
1696 float fnorm[3];
1697 double len;
1698 /* take special care to treat non-orientable surface normals correctly */
1699 if (ELL_3V_DOT(pld->norm+3*disc, pld->norm+3*kept)((pld->norm+3*disc)[0]*(pld->norm+3*kept)[0] + (pld->
norm+3*disc)[1]*(pld->norm+3*kept)[1] + (pld->norm+3*disc
)[2]*(pld->norm+3*kept)[2])<0) {
1700 ELL_3V_SCALE_TT(fnorm, float, -1.0, pld->norm+3*kept)((fnorm)[0] = ((float)((-1.0)*(pld->norm+3*kept)[0])), (fnorm
)[1] = ((float)((-1.0)*(pld->norm+3*kept)[1])), (fnorm)[2]
= ((float)((-1.0)*(pld->norm+3*kept)[2])));
1701 } else {
1702 ELL_3V_COPY(fnorm, pld->norm+3*kept)((fnorm)[0] = (pld->norm+3*kept)[0], (fnorm)[1] = (pld->
norm+3*kept)[1], (fnorm)[2] = (pld->norm+3*kept)[2]);
1703 }
1704 airArrayLenIncr(normArr, 1);
1705 ELL_3V_LERP_TT(newpld->norm+3*q, float, alpha, pld->norm+3*disc, fnorm)((newpld->norm+3*q)[0] = ((float)((((alpha))*(((fnorm)[0])
- ((pld->norm+3*disc)[0])) + ((pld->norm+3*disc)[0])))
), (newpld->norm+3*q)[1] = ((float)((((alpha))*(((fnorm)[1
]) - ((pld->norm+3*disc)[1])) + ((pld->norm+3*disc)[1])
))), (newpld->norm+3*q)[2] = ((float)((((alpha))*(((fnorm)
[2]) - ((pld->norm+3*disc)[2])) + ((pld->norm+3*disc)[2
])))));
1706 /* re-normalize */
1707 len=ELL_3V_LEN(newpld->norm+3*q)(sqrt((((newpld->norm+3*q))[0]*((newpld->norm+3*q))[0] +
((newpld->norm+3*q))[1]*((newpld->norm+3*q))[1] + ((newpld
->norm+3*q))[2]*((newpld->norm+3*q))[2])));
1708 if (len>1e-20) {
1709 ELL_3V_SCALE_TT(newpld->norm+3*q, float, 1.0/len, newpld->norm+3*q)((newpld->norm+3*q)[0] = ((float)((1.0/len)*(newpld->norm
+3*q)[0])), (newpld->norm+3*q)[1] = ((float)((1.0/len)*(newpld
->norm+3*q)[1])), (newpld->norm+3*q)[2] = ((float)((1.0
/len)*(newpld->norm+3*q)[2])));
1710 }
1711 }
1712 if ((1 << limnPolyDataInfoTex2) & bitflag) {
1713 airArrayLenIncr(tex2Arr, 1);
1714 ELL_2V_LERP_TT(newpld->tex2+2*q, float, alpha, pld->tex2+2*disc, pld->tex2+2*kept)((newpld->tex2+2*q)[0] = ((float)((((alpha))*(((pld->tex2
+2*kept)[0]) - ((pld->tex2+2*disc)[0])) + ((pld->tex2+2
*disc)[0])))), (newpld->tex2+2*q)[1] = ((float)((((alpha))
*(((pld->tex2+2*kept)[1]) - ((pld->tex2+2*disc)[1])) + (
(pld->tex2+2*disc)[1])))));
1715 }
1716 if ((1 << limnPolyDataInfoTang) & bitflag) {
1717 airArrayLenIncr(tangArr, 1);
1718 ELL_3V_LERP_TT(newpld->tang+3*q, float, alpha, pld->tang+3*disc, pld->tang+3*kept)((newpld->tang+3*q)[0] = ((float)((((alpha))*(((pld->tang
+3*kept)[0]) - ((pld->tang+3*disc)[0])) + ((pld->tang+3
*disc)[0])))), (newpld->tang+3*q)[1] = ((float)((((alpha))
*(((pld->tang+3*kept)[1]) - ((pld->tang+3*disc)[1])) + (
(pld->tang+3*disc)[1])))), (newpld->tang+3*q)[2] = ((float
)((((alpha))*(((pld->tang+3*kept)[2]) - ((pld->tang+3*disc
)[2])) + ((pld->tang+3*disc)[2])))));
1719 }
1720 /* add new vertex to linked list */
1721 p=airArrayLenIncr(llistArr, 1);
1722 llist=(int*)llistArr->data; /* update in case of re-allocation */
1723 llist[3*p]=kept;
1724 llist[3*p+1]=q;
1725 llist[3*p+2]=-1;
1726 if (ref==-1) newIdx[disc]=3*p;
1727 else llist[ref]=3*p;
1728 return q;
1729}
1730
1731/*
1732 * Clips the given triangles (limnPrimitiveTriangles) according to the
1733 * input matrix nval and the threshold array thresh. First axis of
1734 * nval are different clipping criteria, second axis are vertex
1735 * indices. The length of thresh has to equal the size of the first
1736 * axis, the vertex count in pld has to equal the size of the second
1737 * axis. If nval is 1D, it is assumed to have a single criterion.
1738 *
1739 * A vertex is preserved if all values are >= the respective
1740 * threshold; triangles with partially discarded vertices are clipped,
1741 * potentially generating a quad that is then triangulated arbitrarily.
1742 */
1743int
1744limnPolyDataClipMulti(limnPolyData *pld, Nrrd *nval, double *thresh) {
1745 static const char me[]="limnPolyDataClipMulti";
1746 unsigned char *keepVert=NULL((void*)0);
1747 airArray *mop;
1748 unsigned int E, i, idx=0;
1749 double (*lup)(const void *v, size_t I);
1750 airArray *xyzwArr, *rgbaArr=NULL((void*)0), *normArr=NULL((void*)0), *tex2Arr=NULL((void*)0),
1751 *tangArr=NULL((void*)0), *indxArr, *typeArr, *icntArr, *llistArr=NULL((void*)0);
1752 limnPolyData *newpld=NULL((void*)0);
1753 int *newIdx=NULL((void*)0), *llist=NULL((void*)0);
1754 unsigned int bitflag, nk, nvert;
1755 airPtrPtrUnion appu;
1756
1757 if (!(pld && nval)) {
1758 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1759 return 1;
1760 }
1761
1762 if (nrrdTypeBlock == nval->type) {
1763 biffAddf(LIMNlimnBiffKey, "%s: need scalar type (not %s)", me,
1764 airEnumStr(nrrdType, nval->type));
1765 return 1;
1766 }
1767
1768 if (nval->dim==1) {
1769 nk=1; nvert=nval->axis[0].size;
1770 } else if (nval->dim==2) {
1771 nk=nval->axis[0].size; nvert=nval->axis[1].size;
1772 } else {
1773 biffAddf(LIMNlimnBiffKey, "%s: need 1D or 2D input array, got %uD", me, nval->dim);
1774 return 1;
1775 }
1776 if (nvert!=pld->xyzwNum) {
1777 biffAddf(LIMNlimnBiffKey, "%s: # verts %u != # values %u", me,
1778 pld->xyzwNum, nvert);
1779 return 1;
1780 }
1781 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pld)) {
1782 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
1783 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
1784 return 1;
1785 }
1786
1787 /* Memory allocation in C IS a headache */
1788 mop=airMopNew();
1789 E = AIR_FALSE0;
1790 if (!E) {
1791 E|=!(keepVert = AIR_CAST(unsigned char *,((unsigned char *)(calloc(pld->xyzwNum, sizeof(char))))
1792 calloc(pld->xyzwNum, sizeof(char)))((unsigned char *)(calloc(pld->xyzwNum, sizeof(char)))));
1793 }
1794 if (!E) {
1795 airMopAdd(mop, keepVert, airFree, airMopAlways);
1796 E|=!(newIdx = AIR_CAST(int *, malloc(pld->xyzwNum*sizeof(int)))((int *)(malloc(pld->xyzwNum*sizeof(int)))));
1797 }
1798 if (!E) {
1799 unsigned int incr;
1800 airMopAdd(mop, newIdx, airFree, airMopAlways);
1801 memset(newIdx, -1, sizeof(int)*pld->xyzwNum)__builtin___memset_chk (newIdx, -1, sizeof(int)*pld->xyzwNum
, __builtin_object_size (newIdx, 0));
1802 /* This setting of incr is arbitrary and was not optimized in any way: */
1803 incr = pld->xyzwNum/10; /* 10% of previous vertex count... */
1804 if (incr<50) incr=50; /* ...but at least 50. */
1805 appu.i = &llist;
1806 E|=!(llistArr=airArrayNew(appu.v, NULL((void*)0), 3*sizeof(int), incr));
1807 }
1808 if (!E) {
1809 airMopAdd(mop, llistArr, (airMopper)airArrayNuke, airMopAlways);
1810 E|=!(newpld = limnPolyDataNew());
1811 }
1812 bitflag = limnPolyDataInfoBitFlag(pld);
1813 if (!E) {
1814 unsigned int incr;
1815 airMopAdd(mop, newpld, airFree, airMopAlways); /* "shallow" free */
1816 incr = pld->xyzwNum/20; /* 5% of previous vertex count... */
1817 if (incr<10) incr=10; /* ...but at least 10. */
1818 appu.f = &(newpld->xyzw);
1819 E|=!(xyzwArr=airArrayNew(appu.v, &(newpld->xyzwNum),
1820 4*sizeof(float), incr));
1821 if (!E) {
1822 airMopAdd(mop, xyzwArr, (airMopper)airArrayNuke, airMopOnError);
1823 airMopAdd(mop, xyzwArr, (airMopper)airArrayNix, airMopOnOkay);
1824 }
1825 if (!E && (1 << limnPolyDataInfoRGBA) & bitflag) {
1826 appu.uc = &(newpld->rgba);
1827 E|=!(rgbaArr=airArrayNew(appu.v, &(newpld->rgbaNum),
1828 4*sizeof(unsigned char), incr));
1829 if (!E) {
1830 airMopAdd(mop, rgbaArr, (airMopper)airArrayNuke, airMopOnError);
1831 airMopAdd(mop, rgbaArr, (airMopper)airArrayNix, airMopOnOkay);
1832 }
1833 }
1834 if (!E && (1 << limnPolyDataInfoNorm) & bitflag) {
1835 appu.f = &(newpld->norm);
1836 E|=!(normArr=airArrayNew(appu.v, &(newpld->normNum),
1837 3*sizeof(float), incr));
1838 if (!E) {
1839 airMopAdd(mop, normArr, (airMopper)airArrayNuke, airMopOnError);
1840 airMopAdd(mop, normArr, (airMopper)airArrayNix, airMopOnOkay);
1841 }
1842 }
1843 if (!E && (1 << limnPolyDataInfoTex2) & bitflag) {
1844 appu.f = &(newpld->tex2);
1845 E|=!(tex2Arr=airArrayNew(appu.v, &(newpld->tex2Num),
1846 2*sizeof(float), incr));
1847 if (!E) {
1848 airMopAdd(mop, tex2Arr, (airMopper)airArrayNuke, airMopOnError);
1849 airMopAdd(mop, tex2Arr, (airMopper)airArrayNix, airMopOnOkay);
1850 }
1851 }
1852 if (!E && (1 << limnPolyDataInfoTang) & bitflag) {
1853 appu.f = &(newpld->tang);
1854 E|=!(tangArr=airArrayNew(appu.v, &(newpld->tangNum),
1855 3*sizeof(float), incr));
1856 if (!E) {
1857 airMopAdd(mop, tangArr, (airMopper)airArrayNuke, airMopOnError);
1858 airMopAdd(mop, tangArr, (airMopper)airArrayNix, airMopOnOkay);
1859 }
1860 }
1861 if (!E) {
1862 incr = pld->indxNum/20; /* 5% of previous index count... */
1863 if (incr<10) incr=10; /* ...but at least 10. */
1864 appu.ui = &(newpld->indx);
1865 E|=!(indxArr=airArrayNew(appu.v, &(newpld->indxNum),
1866 sizeof(unsigned int), incr));
1867 if (!E) {
1868 airMopAdd(mop, indxArr, (airMopper)airArrayNuke, airMopOnError);
1869 airMopAdd(mop, indxArr, (airMopper)airArrayNix, airMopOnOkay);
1870 }
1871 }
1872 if (!E) {
1873 incr = pld->primNum/10; /* 10% of previous primNum... */
1874 if (incr<1) incr=1; /* ...but at least 1. */
1875 appu.uc = &(newpld->type);
1876 E|=!(typeArr=airArrayNew(appu.v, &(newpld->primNum),
1877 sizeof(unsigned char), incr));
1878 if (!E) {
1879 airMopAdd(mop, typeArr, (airMopper)airArrayNuke, airMopOnError);
1880 airMopAdd(mop, typeArr, (airMopper)airArrayNix, airMopOnOkay);
1881 }
1882 appu.ui = &(newpld->icnt);
1883 E|=!(icntArr=airArrayNew(appu.v, NULL((void*)0),
1884 sizeof(unsigned int), incr));
1885 if (!E) {
1886 airMopAdd(mop, icntArr, (airMopper)airArrayNuke, airMopOnError);
1887 airMopAdd(mop, icntArr, (airMopper)airArrayNix, airMopOnOkay);
1888 }
1889 }
1890 }
1891 if (E) {
1892 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate buffers", me);
1893 airMopError(mop); return 1;
1894 }
1895
1896 /* mark vertices, leaving at 0 means "discard" */
1897 lup = nrrdDLookup[nval->type];
1898 for (i=0; i<pld->xyzwNum; i++) {
1899 unsigned int j, keep = AIR_TRUE1;
1900 for (j=0; j<nk; j++, idx++) {
1901 if (lup(nval->data, idx) < thresh[j])
1902 keep = AIR_FALSE0;
1903 }
1904 if (keep) {
1905 keepVert[i]=AIR_TRUE1;
1906 }
1907 }
1908
1909 /* now, iterate over all primitives and triangles */
1910
1911 /* Note: If keepVert[i]==AIR_TRUE, newIdx[i] is its new index; else, it is
1912 * an index j into llist, which is a linked list:
1913 * llist[j] == other (kept) end of the edge
1914 * llist[j+1] == index of new vertex for that edge
1915 * llist[j+2] == next index into llist
1916 */
1917
1918 /* TODO: All the airArray stuff should have allocation error checking */
1919 idx=0;
1920 for (i=0; i<pld->primNum; i++) {
1921 int j, oldTriNum=pld->icnt[i]/3, newTriNum=0;
1922 unsigned int kept=0, disck=0; /* index of last kept / discarded vertex */
1923 for (j=0; j<oldTriNum; j++, idx+=3) {
1924 unsigned int p, quad[4];
1925 int k, keepN=0;
1926 for (k=0; k<3; k++) {
1927 unsigned int oldidx=pld->indx[idx+k];
1928 if (keepVert[oldidx]) {
1929 keepN++; kept=oldidx;
1930 /* make sure the vertex is copied over */
1931 if (newIdx[oldidx]==-1) {
1932 unsigned int q=newIdx[oldidx]=airArrayLenIncr(xyzwArr, 1);
1933 ELL_4V_COPY(newpld->xyzw+4*q, pld->xyzw+4*oldidx)((newpld->xyzw+4*q)[0] = (pld->xyzw+4*oldidx)[0], (newpld
->xyzw+4*q)[1] = (pld->xyzw+4*oldidx)[1], (newpld->xyzw
+4*q)[2] = (pld->xyzw+4*oldidx)[2], (newpld->xyzw+4*q)[
3] = (pld->xyzw+4*oldidx)[3]);
1934 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1935 airArrayLenIncr(rgbaArr, 1);
1936 ELL_4V_COPY(newpld->rgba+4*q, pld->rgba+4*oldidx)((newpld->rgba+4*q)[0] = (pld->rgba+4*oldidx)[0], (newpld
->rgba+4*q)[1] = (pld->rgba+4*oldidx)[1], (newpld->rgba
+4*q)[2] = (pld->rgba+4*oldidx)[2], (newpld->rgba+4*q)[
3] = (pld->rgba+4*oldidx)[3]);
1937 }
1938 if ((1 << limnPolyDataInfoNorm) & bitflag) {
1939 airArrayLenIncr(normArr, 1);
1940 ELL_3V_COPY(newpld->norm+3*q, pld->norm+3*oldidx)((newpld->norm+3*q)[0] = (pld->norm+3*oldidx)[0], (newpld
->norm+3*q)[1] = (pld->norm+3*oldidx)[1], (newpld->norm
+3*q)[2] = (pld->norm+3*oldidx)[2]);
1941 }
1942 if ((1 << limnPolyDataInfoTex2) & bitflag) {
1943 airArrayLenIncr(tex2Arr, 1);
1944 ELL_2V_COPY(newpld->tex2+2*q, pld->tex2+2*oldidx)((newpld->tex2+2*q)[0] = (pld->tex2+2*oldidx)[0], (newpld
->tex2+2*q)[1] = (pld->tex2+2*oldidx)[1]);
1945 }
1946 if ((1 << limnPolyDataInfoTang) & bitflag) {
1947 airArrayLenIncr(tangArr, 1);
1948 ELL_3V_COPY(newpld->tang+3*q, pld->tang+3*oldidx)((newpld->tang+3*q)[0] = (pld->tang+3*oldidx)[0], (newpld
->tang+3*q)[1] = (pld->tang+3*oldidx)[1], (newpld->tang
+3*q)[2] = (pld->tang+3*oldidx)[2]);
1949 }
1950 }
1951 } else {
1952 disck=k;
1953 }
1954 }
1955 switch (keepN) {
1956 case 0: /* nothing to be done; discard this triangle */
1957 break;
1958 case 1: /* result of clipping is a single triangle */
1959 newTriNum++;
1960 p=airArrayLenIncr(indxArr, 3);
1961 for (k=0; k<3; k++) {
1962 if (keepVert[pld->indx[idx+k]])
1963 newpld->indx[p+k]=newIdx[pld->indx[idx+k]];
1964 else
1965 newpld->indx[p+k]=clipEdge(pld->indx[idx+k], kept, nval, thresh,
1966 newIdx, llistArr, pld, bitflag, newpld,
1967 xyzwArr, rgbaArr, normArr,
1968 tex2Arr, tangArr);
1969 }
1970 break;
1971 case 2: /* result of clipping is a quad, triangulate */
1972 newTriNum+=2;
1973 p=0;
1974 for (k=0; k<3; k++) {
1975 if (keepVert[pld->indx[idx+k]]) quad[p++]=newIdx[pld->indx[idx+k]];
1976 else {
1977 quad[p++]=clipEdge(pld->indx[idx+k], pld->indx[idx+(disck+2)%3],
1978 nval, thresh, newIdx, llistArr,
1979 pld, bitflag, newpld, xyzwArr, rgbaArr,
1980 normArr, tex2Arr, tangArr);
1981 quad[p++]=clipEdge(pld->indx[idx+k], pld->indx[idx+(disck+1)%3],
1982 nval, thresh, newIdx, llistArr,
1983 pld, bitflag, newpld, xyzwArr, rgbaArr,
1984 normArr, tex2Arr, tangArr);
1985 }
1986 }
1987 p=airArrayLenIncr(indxArr, 6);
1988 ELL_3V_SET(newpld->indx+p, quad[0], quad[1], quad[3])((newpld->indx+p)[0] = (quad[0]), (newpld->indx+p)[1] =
(quad[1]), (newpld->indx+p)[2] = (quad[3]));
1989 ELL_3V_SET(newpld->indx+p+3, quad[1], quad[2], quad[3])((newpld->indx+p+3)[0] = (quad[1]), (newpld->indx+p+3)[
1] = (quad[2]), (newpld->indx+p+3)[2] = (quad[3]));
1990 break;
1991 case 3: /* simply copy the existing triangle */
1992 newTriNum++;
1993 p=airArrayLenIncr(indxArr, 3);
1994 for (k=0; k<3; k++) {
1995 newpld->indx[p+k]=newIdx[pld->indx[idx+k]];
1996 }
1997 break;
1998 }
1999 }
2000 if (newTriNum>0) {
2001 unsigned int p=airArrayLenIncr(typeArr, 1);
2002 airArrayLenIncr(icntArr, 1);
2003 newpld->type[p]=limnPrimitiveTriangles;
2004 newpld->icnt[p]=newTriNum*3;
2005 }
2006 }
2007
2008 /* finally, replace contents of pld with new data */
2009 airFree(pld->xyzw);
2010 airFree(pld->rgba);
2011 airFree(pld->norm);
2012 airFree(pld->tex2);
2013 airFree(pld->tang);
2014 airFree(pld->indx);
2015 airFree(pld->type);
2016 airFree(pld->icnt);
2017 memcpy(pld, newpld, sizeof(limnPolyData))__builtin___memcpy_chk (pld, newpld, sizeof(limnPolyData), __builtin_object_size
(pld, 0));
2018
2019 airMopOkay(mop);
2020 return 0;
2021}
2022
2023/* Simple wrapper around limnPolyDataClipMulti, in case of only one
2024 * clipping criterion.
2025 */
2026int
2027limnPolyDataClip(limnPolyData *pld, Nrrd *nval, double thresh) {
2028 return limnPolyDataClipMulti(pld, nval, &thresh);
2029}
2030
2031/* limnPolyDataCompress:
2032 * returns a "compressed" copy of the given limnPolyData pld that only
2033 * contains vertices that are referenced by some primitive
2034 * returns NULL and adds a message to biff upon error
2035 */
2036limnPolyData *
2037limnPolyDataCompress(const limnPolyData *pld) {
2038 static const char me[]="limnPolyDataCompress";
2039 limnPolyData *ret = NULL((void*)0);
2040 unsigned int infoBitFlag=0, vertNum=0, i, used_indxNum=0;
2041 int *vertMap;
2042 if (pld==NULL((void*)0)) {
2043 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
2044 return NULL((void*)0);
2045 }
2046 infoBitFlag=limnPolyDataInfoBitFlag(pld);
2047 vertMap=(int*) calloc(pld->xyzwNum,sizeof(int));
2048 if (vertMap==NULL((void*)0)) {
2049 biffAddf(LIMNlimnBiffKey, "%s: could not allocate memory", me);
2050 return NULL((void*)0);
2051 }
2052 /* how many indices are actually used? */
2053 for (i=0; i<pld->primNum; i++) {
2054 used_indxNum+=pld->icnt[i];
2055 }
2056 /* loop over all indices and mark referenced vertices in vertMap */
2057 for (i=0; i<used_indxNum; i++) {
2058 vertMap[pld->indx[i]]=1;
2059 }
2060 /* turn vertMap into an index map */
2061 for (i=0; i<pld->xyzwNum; i++) {
2062 if (vertMap[i]==0)
2063 vertMap[i]=-1;
2064 else
2065 vertMap[i]=vertNum++;
2066 }
2067 /* allocate new limnPolyData */
2068 if (NULL((void*)0) == (ret=limnPolyDataNew()) ||
2069 0!=limnPolyDataAlloc(ret, infoBitFlag, vertNum,
2070 used_indxNum, pld->primNum)) {
2071 biffAddf(LIMNlimnBiffKey, "%s: Could not allocate result", me);
2072 free(vertMap);
2073 return NULL((void*)0);
2074 }
2075 /* fill the newly allocated structure */
2076 for (i=0; i<pld->xyzwNum; i++) {
2077 if (vertMap[i]>=0) {
2078 ELL_4V_COPY(ret->xyzw+4*vertMap[i], pld->xyzw+4*i)((ret->xyzw+4*vertMap[i])[0] = (pld->xyzw+4*i)[0], (ret
->xyzw+4*vertMap[i])[1] = (pld->xyzw+4*i)[1], (ret->
xyzw+4*vertMap[i])[2] = (pld->xyzw+4*i)[2], (ret->xyzw+
4*vertMap[i])[3] = (pld->xyzw+4*i)[3]);
2079 }
2080 }
2081 if (ret->rgba!=NULL((void*)0)) {
2082 for (i=0; i<pld->xyzwNum; i++) {
2083 if (vertMap[i]>=0) {
2084 ELL_4V_COPY(ret->rgba+4*vertMap[i], pld->rgba+4*i)((ret->rgba+4*vertMap[i])[0] = (pld->rgba+4*i)[0], (ret
->rgba+4*vertMap[i])[1] = (pld->rgba+4*i)[1], (ret->
rgba+4*vertMap[i])[2] = (pld->rgba+4*i)[2], (ret->rgba+
4*vertMap[i])[3] = (pld->rgba+4*i)[3]);
2085 }
2086 }
2087 }
2088 if (ret->norm!=NULL((void*)0)) {
2089 for (i=0; i<pld->xyzwNum; i++) {
2090 if (vertMap[i]>=0) {
2091 ELL_3V_COPY(ret->norm+3*vertMap[i], pld->norm+3*i)((ret->norm+3*vertMap[i])[0] = (pld->norm+3*i)[0], (ret
->norm+3*vertMap[i])[1] = (pld->norm+3*i)[1], (ret->
norm+3*vertMap[i])[2] = (pld->norm+3*i)[2]);
2092 }
2093 }
2094 }
2095 if (ret->tex2!=NULL((void*)0)) {
2096 for (i=0; i<pld->xyzwNum; i++) {
2097 if (vertMap[i]>=0) {
2098 ELL_2V_COPY(ret->tex2+2*vertMap[i], pld->tex2+2*i)((ret->tex2+2*vertMap[i])[0] = (pld->tex2+2*i)[0], (ret
->tex2+2*vertMap[i])[1] = (pld->tex2+2*i)[1]);
2099 }
2100 }
2101 }
2102 if (ret->tang!=NULL((void*)0)) {
2103 for (i=0; i<pld->xyzwNum; i++) {
2104 if (vertMap[i]>=0) {
2105 ELL_3V_COPY(ret->tang+3*vertMap[i], pld->tang+3*i)((ret->tang+3*vertMap[i])[0] = (pld->tang+3*i)[0], (ret
->tang+3*vertMap[i])[1] = (pld->tang+3*i)[1], (ret->
tang+3*vertMap[i])[2] = (pld->tang+3*i)[2]);
2106 }
2107 }
2108 }
2109 for (i=0; i<used_indxNum; i++) {
2110 ret->indx[i]=vertMap[pld->indx[i]];
2111 }
2112 memcpy(ret->type, pld->type, sizeof(char)*pld->primNum)__builtin___memcpy_chk (ret->type, pld->type, sizeof(char
)*pld->primNum, __builtin_object_size (ret->type, 0));
2113 memcpy(ret->icnt, pld->icnt, sizeof(int)*pld->primNum)__builtin___memcpy_chk (ret->icnt, pld->icnt, sizeof(int
)*pld->primNum, __builtin_object_size (ret->icnt, 0));
2114
2115 free(vertMap);
2116
2117 return ret;
2118}
2119
2120/* limnPolyDataJoin:
2121 * concatenates the primitives in all num limnPolyDatas given in plds
2122 * and returns the result as a newly allocated limnPolyData
2123 * the new limnPolyData will only have color/normals/texture coordinates
2124 * if _all_ input limnPolyDatas had the respective attribute
2125 * returns NULL and adds a message to biff upon error
2126 */
2127limnPolyData *limnPolyDataJoin(const limnPolyData **plds,
2128 unsigned int num) {
2129 static const char me[]="limnPolyDataJoin";
2130 limnPolyData *ret = NULL((void*)0);
2131 unsigned int infoBitFlag=(1 << limnPolyDataInfoRGBA) |
2132 (1 << limnPolyDataInfoNorm) |
2133 (1 << limnPolyDataInfoTex2) |
2134 (1 << limnPolyDataInfoTang); /* by default, assume we have all these */
2135 unsigned int vertNum=0, indxNum=0, primNum=0;
2136 unsigned int i;
2137 if (plds==NULL((void*)0)) {
2138 biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
2139 return NULL((void*)0);
2140 }
2141 /* loop over all input plds to find infoBitFlag and the total number of
2142 * vertices / indices / primitives */
2143 for (i=0; i<num; i++) {
2144 if (plds[i]==NULL((void*)0)) {
2145 biffAddf(LIMNlimnBiffKey, "%s: plds[%d] is a NULL pointer", me, i);
2146 return NULL((void*)0);
2147 }
2148 infoBitFlag &= limnPolyDataInfoBitFlag(plds[i]);
2149 vertNum += plds[i]->xyzwNum;
2150 indxNum += plds[i]->indxNum;
2151 primNum += plds[i]->primNum;
2152 }
2153 if (NULL((void*)0) == (ret=limnPolyDataNew()) ||
2154 0!=limnPolyDataAlloc(ret, infoBitFlag, vertNum, indxNum, primNum)) {
2155 biffAddf(LIMNlimnBiffKey, "%s: Could not allocate result", me);
2156 return NULL((void*)0);
2157 }
2158 /* loop again over all input plds and fill the newly allocated structure */
2159 vertNum=indxNum=primNum=0;
2160 for (i=0; i<num; i++) {
2161 unsigned int j, used_indxNum=0;
2162 memcpy(ret->xyzw+4*vertNum, plds[i]->xyzw,__builtin___memcpy_chk (ret->xyzw+4*vertNum, plds[i]->xyzw
, sizeof(float)*4*plds[i]->xyzwNum, __builtin_object_size (
ret->xyzw+4*vertNum, 0))
2163 sizeof(float)*4*plds[i]->xyzwNum)__builtin___memcpy_chk (ret->xyzw+4*vertNum, plds[i]->xyzw
, sizeof(float)*4*plds[i]->xyzwNum, __builtin_object_size (
ret->xyzw+4*vertNum, 0));
2164 if (ret->rgba!=NULL((void*)0)) {
2165 memcpy(ret->rgba+4*vertNum, plds[i]->rgba,__builtin___memcpy_chk (ret->rgba+4*vertNum, plds[i]->rgba
, sizeof(unsigned char)*4*plds[i]->xyzwNum, __builtin_object_size
(ret->rgba+4*vertNum, 0))
2166 sizeof(unsigned char)*4*plds[i]->xyzwNum)__builtin___memcpy_chk (ret->rgba+4*vertNum, plds[i]->rgba
, sizeof(unsigned char)*4*plds[i]->xyzwNum, __builtin_object_size
(ret->rgba+4*vertNum, 0));
2167 }
2168 if (ret->norm!=NULL((void*)0)) {
2169 memcpy(ret->norm+3*vertNum, plds[i]->norm,__builtin___memcpy_chk (ret->norm+3*vertNum, plds[i]->norm
, sizeof(float)*3*plds[i]->xyzwNum, __builtin_object_size (
ret->norm+3*vertNum, 0))
2170 sizeof(float)*3*plds[i]->xyzwNum)__builtin___memcpy_chk (ret->norm+3*vertNum, plds[i]->norm
, sizeof(float)*3*plds[i]->xyzwNum, __builtin_object_size (
ret->norm+3*vertNum, 0));
2171 }
2172 if (ret->tex2!=NULL((void*)0)) {
2173 memcpy(ret->tex2+2*vertNum, plds[i]->tex2,__builtin___memcpy_chk (ret->tex2+2*vertNum, plds[i]->tex2
, sizeof(float)*2*plds[i]->xyzwNum, __builtin_object_size (
ret->tex2+2*vertNum, 0))
2174 sizeof(float)*2*plds[i]->xyzwNum)__builtin___memcpy_chk (ret->tex2+2*vertNum, plds[i]->tex2
, sizeof(float)*2*plds[i]->xyzwNum, __builtin_object_size (
ret->tex2+2*vertNum, 0));
2175 }
2176 if (ret->tang!=NULL((void*)0)) {
2177 memcpy(ret->tang+3*vertNum, plds[i]->tang,__builtin___memcpy_chk (ret->tang+3*vertNum, plds[i]->tang
, sizeof(float)*3*plds[i]->xyzwNum, __builtin_object_size (
ret->tang+3*vertNum, 0))
2178 sizeof(float)*3*plds[i]->xyzwNum)__builtin___memcpy_chk (ret->tang+3*vertNum, plds[i]->tang
, sizeof(float)*3*plds[i]->xyzwNum, __builtin_object_size (
ret->tang+3*vertNum, 0));
2179 }
2180 for (j=0; j<plds[i]->indxNum; j++) {
2181 ret->indx[indxNum+j]=vertNum+plds[i]->indx[j];
2182 }
2183 for (j=0; j<plds[i]->primNum; j++) {
2184 ret->type[primNum+j]=plds[i]->type[j];
2185 ret->icnt[primNum+j]=plds[i]->icnt[j];
2186 /* need to keep track of how many indices are actually used */
2187 used_indxNum+=plds[i]->icnt[j];
2188 }
2189 vertNum+=plds[i]->xyzwNum;
2190 indxNum+=used_indxNum;
2191 primNum+=plds[i]->primNum;
2192 }
2193 return ret;
2194}
2195
2196int
2197limnPolyDataEdgeHalve(limnPolyData *pldOut,
2198 const limnPolyData *pldIn) {
2199 static const char me[]="limnPolyDataEdgeHalve";
2200 Nrrd *nnewvert;
2201 unsigned int *newvert, nvold, nvidx, triidx, trinum, vlo, vhi, bitflag;
2202 airArray *mop;
2203
2204 if ((1 << limnPrimitiveTriangles) != limnPolyDataPrimitiveTypes(pldIn)) {
2205 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle %s primitives", me,
2206 airEnumStr(limnPrimitive, limnPrimitiveTriangles));
2207 return 1;
2208 }
2209 if (1 != pldIn->primNum) {
2210 biffAddf(LIMNlimnBiffKey, "%s: sorry, can only handle a single primitive", me);
2211 return 1;
2212 }
2213 mop = airMopNew();
2214 nnewvert = nrrdNew();
2215 airMopAdd(mop, nnewvert, AIR_CAST(airMopper, nrrdNuke)((airMopper)(nrrdNuke)), airMopAlways);
2216 nvold = pldIn->xyzwNum;
2217 if (nrrdMaybeAlloc_va(nnewvert, nrrdTypeUInt, 2, nvold, nvold)) {
2218 biffMovef(LIMNlimnBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate buffer", me);
2219 airMopError(mop); return 1;
2220 }
2221 newvert = AIR_CAST(unsigned int*, nnewvert->data)((unsigned int*)(nnewvert->data));
2222
2223 /* run through triangles, recording edges with the new vertex index */
2224 nvidx = nvold;
2225 trinum = pldIn->indxNum/3;
2226 for (triidx=0; triidx<trinum; triidx++) {
2227 vlo = pldIn->indx[0 + 3*triidx];
2228 vhi = pldIn->indx[1 + 3*triidx];
2229 if (!newvert[vlo + nvold*vhi]) {
2230 newvert[vlo + nvold*vhi] = newvert[vhi + nvold*vlo] = nvidx++;
2231 }
2232 vlo = pldIn->indx[1 + 3*triidx];
2233 vhi = pldIn->indx[2 + 3*triidx];
2234 if (!newvert[vlo + nvold*vhi]) {
2235 newvert[vlo + nvold*vhi] = newvert[vhi + nvold*vlo] = nvidx++;
2236 }
2237 vlo = pldIn->indx[2 + 3*triidx];
2238 vhi = pldIn->indx[0 + 3*triidx];
2239 if (!newvert[vlo + nvold*vhi]) {
2240 newvert[vlo + nvold*vhi] = newvert[vhi + nvold*vlo] = nvidx++;
2241 }
2242 }
2243
2244 /* allocate output */
2245 bitflag = limnPolyDataInfoBitFlag(pldIn);
2246 if (limnPolyDataAlloc(pldOut, bitflag, nvidx, 3*4*trinum, 1)) {
2247 biffAddf(LIMNlimnBiffKey, "%s: trouble allocating output", me);
2248 airMopError(mop); return 1;
2249 }
2250 pldOut->type[0] = limnPrimitiveTriangles;
2251 pldOut->icnt[0] = 3*4*trinum;
2252
2253 /* set output indx */
2254 for (triidx=0; triidx<trinum; triidx++) {
2255 unsigned int aa, ab, bb, bc, cc, ac;
2256 aa = pldIn->indx[0 + 3*triidx];
2257 bb = pldIn->indx[1 + 3*triidx];
2258 cc = pldIn->indx[2 + 3*triidx];
2259 ab = newvert[aa + nvold*bb];
2260 bc = newvert[bb + nvold*cc];
2261 ac = newvert[aa + nvold*cc];
2262 ELL_3V_SET(pldOut->indx + 3*(0 + 4*triidx), aa, ab, ac)((pldOut->indx + 3*(0 + 4*triidx))[0] = (aa), (pldOut->
indx + 3*(0 + 4*triidx))[1] = (ab), (pldOut->indx + 3*(0 +
4*triidx))[2] = (ac));
2263 ELL_3V_SET(pldOut->indx + 3*(1 + 4*triidx), ab, bc, ac)((pldOut->indx + 3*(1 + 4*triidx))[0] = (ab), (pldOut->
indx + 3*(1 + 4*triidx))[1] = (bc), (pldOut->indx + 3*(1 +
4*triidx))[2] = (ac));
2264 ELL_3V_SET(pldOut->indx + 3*(2 + 4*triidx), ab, bb, bc)((pldOut->indx + 3*(2 + 4*triidx))[0] = (ab), (pldOut->
indx + 3*(2 + 4*triidx))[1] = (bb), (pldOut->indx + 3*(2 +
4*triidx))[2] = (bc));
2265 ELL_3V_SET(pldOut->indx + 3*(3 + 4*triidx), ac, bc, cc)((pldOut->indx + 3*(3 + 4*triidx))[0] = (ac), (pldOut->
indx + 3*(3 + 4*triidx))[1] = (bc), (pldOut->indx + 3*(3 +
4*triidx))[2] = (cc));
2266 }
2267
2268 /* set output vertex info */
2269 for (vlo=0; vlo<nvold; vlo++) {
2270 ELL_4V_COPY(pldOut->xyzw + 4*vlo, pldIn->xyzw + 4*vlo)((pldOut->xyzw + 4*vlo)[0] = (pldIn->xyzw + 4*vlo)[0], (
pldOut->xyzw + 4*vlo)[1] = (pldIn->xyzw + 4*vlo)[1], (pldOut
->xyzw + 4*vlo)[2] = (pldIn->xyzw + 4*vlo)[2], (pldOut->
xyzw + 4*vlo)[3] = (pldIn->xyzw + 4*vlo)[3]);
2271 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
2272 ELL_4V_COPY(pldOut->rgba + 4*vlo, pldIn->rgba + 4*vlo)((pldOut->rgba + 4*vlo)[0] = (pldIn->rgba + 4*vlo)[0], (
pldOut->rgba + 4*vlo)[1] = (pldIn->rgba + 4*vlo)[1], (pldOut
->rgba + 4*vlo)[2] = (pldIn->rgba + 4*vlo)[2], (pldOut->
rgba + 4*vlo)[3] = (pldIn->rgba + 4*vlo)[3]);
2273 }
2274 if ((1 << limnPolyDataInfoNorm) & bitflag) {
2275 ELL_3V_COPY(pldOut->norm + 3*vlo, pldIn->norm + 3*vlo)((pldOut->norm + 3*vlo)[0] = (pldIn->norm + 3*vlo)[0], (
pldOut->norm + 3*vlo)[1] = (pldIn->norm + 3*vlo)[1], (pldOut
->norm + 3*vlo)[2] = (pldIn->norm + 3*vlo)[2]);
2276 }
2277 if ((1 << limnPolyDataInfoTex2) & bitflag) {
2278 ELL_2V_COPY(pldOut->tex2 + 2*vlo, pldIn->tex2 + 2*vlo)((pldOut->tex2 + 2*vlo)[0] = (pldIn->tex2 + 2*vlo)[0], (
pldOut->tex2 + 2*vlo)[1] = (pldIn->tex2 + 2*vlo)[1]);
2279 }
2280 if ((1 << limnPolyDataInfoTang) & bitflag) {
2281 ELL_3V_COPY(pldOut->tang + 3*vlo, pldIn->tang + 3*vlo)((pldOut->tang + 3*vlo)[0] = (pldIn->tang + 3*vlo)[0], (
pldOut->tang + 3*vlo)[1] = (pldIn->tang + 3*vlo)[1], (pldOut
->tang + 3*vlo)[2] = (pldIn->tang + 3*vlo)[2]);
2282 }
2283 for (vhi=vlo+1; vhi<nvold; vhi++) {
2284 unsigned int mid;
2285 mid = newvert[vlo + nvold*vhi];
2286 if (!mid) {
2287 continue;
2288 }
2289 ELL_4V_LERP(pldOut->xyzw + 4*mid, 0.5f,((pldOut->xyzw + 4*mid)[0] = (((0.5f))*(((pldIn->xyzw +
4*vhi)[0]) - ((pldIn->xyzw + 4*vlo)[0])) + ((pldIn->xyzw
+ 4*vlo)[0])), (pldOut->xyzw + 4*mid)[1] = (((0.5f))*(((pldIn
->xyzw + 4*vhi)[1]) - ((pldIn->xyzw + 4*vlo)[1])) + ((pldIn
->xyzw + 4*vlo)[1])), (pldOut->xyzw + 4*mid)[2] = (((0.5f
))*(((pldIn->xyzw + 4*vhi)[2]) - ((pldIn->xyzw + 4*vlo)
[2])) + ((pldIn->xyzw + 4*vlo)[2])), (pldOut->xyzw + 4*
mid)[3] = (((0.5f))*(((pldIn->xyzw + 4*vhi)[3]) - ((pldIn->
xyzw + 4*vlo)[3])) + ((pldIn->xyzw + 4*vlo)[3])))
2290 pldIn->xyzw + 4*vlo,((pldOut->xyzw + 4*mid)[0] = (((0.5f))*(((pldIn->xyzw +
4*vhi)[0]) - ((pldIn->xyzw + 4*vlo)[0])) + ((pldIn->xyzw
+ 4*vlo)[0])), (pldOut->xyzw + 4*mid)[1] = (((0.5f))*(((pldIn
->xyzw + 4*vhi)[1]) - ((pldIn->xyzw + 4*vlo)[1])) + ((pldIn
->xyzw + 4*vlo)[1])), (pldOut->xyzw + 4*mid)[2] = (((0.5f
))*(((pldIn->xyzw + 4*vhi)[2]) - ((pldIn->xyzw + 4*vlo)
[2])) + ((pldIn->xyzw + 4*vlo)[2])), (pldOut->xyzw + 4*
mid)[3] = (((0.5f))*(((pldIn->xyzw + 4*vhi)[3]) - ((pldIn->
xyzw + 4*vlo)[3])) + ((pldIn->xyzw + 4*vlo)[3])))
2291 pldIn->xyzw + 4*vhi)((pldOut->xyzw + 4*mid)[0] = (((0.5f))*(((pldIn->xyzw +
4*vhi)[0]) - ((pldIn->xyzw + 4*vlo)[0])) + ((pldIn->xyzw
+ 4*vlo)[0])), (pldOut->xyzw + 4*mid)[1] = (((0.5f))*(((pldIn
->xyzw + 4*vhi)[1]) - ((pldIn->xyzw + 4*vlo)[1])) + ((pldIn
->xyzw + 4*vlo)[1])), (pldOut->xyzw + 4*mid)[2] = (((0.5f
))*(((pldIn->xyzw + 4*vhi)[2]) - ((pldIn->xyzw + 4*vlo)
[2])) + ((pldIn->xyzw + 4*vlo)[2])), (pldOut->xyzw + 4*
mid)[3] = (((0.5f))*(((pldIn->xyzw + 4*vhi)[3]) - ((pldIn->
xyzw + 4*vlo)[3])) + ((pldIn->xyzw + 4*vlo)[3])));
2292 if ((1 << limnPolyDataInfoRGBA) & bitflag) {
2293 ELL_4V_LERP_TT(pldOut->rgba + 4*mid, unsigned char, 0.5f,((pldOut->rgba + 4*mid)[0] = ((unsigned char)((((0.5f))*((
(pldIn->rgba + 4*vhi)[0]) - ((pldIn->rgba + 4*vlo)[0]))
+ ((pldIn->rgba + 4*vlo)[0])))), (pldOut->rgba + 4*mid
)[1] = ((unsigned char)((((0.5f))*(((pldIn->rgba + 4*vhi)[
1]) - ((pldIn->rgba + 4*vlo)[1])) + ((pldIn->rgba + 4*vlo
)[1])))), (pldOut->rgba + 4*mid)[2] = ((unsigned char)((((
0.5f))*(((pldIn->rgba + 4*vhi)[2]) - ((pldIn->rgba + 4*
vlo)[2])) + ((pldIn->rgba + 4*vlo)[2])))), (pldOut->rgba
+ 4*mid)[3] = ((unsigned char)((((0.5f))*(((pldIn->rgba +
4*vhi)[3]) - ((pldIn->rgba + 4*vlo)[3])) + ((pldIn->rgba
+ 4*vlo)[3])))))
2294 pldIn->rgba + 4*vlo,((pldOut->rgba + 4*mid)[0] = ((unsigned char)((((0.5f))*((
(pldIn->rgba + 4*vhi)[0]) - ((pldIn->rgba + 4*vlo)[0]))
+ ((pldIn->rgba + 4*vlo)[0])))), (pldOut->rgba + 4*mid
)[1] = ((unsigned char)((((0.5f))*(((pldIn->rgba + 4*vhi)[
1]) - ((pldIn->rgba + 4*vlo)[1])) + ((pldIn->rgba + 4*vlo
)[1])))), (pldOut->rgba + 4*mid)[2] = ((unsigned char)((((
0.5f))*(((pldIn->rgba + 4*vhi)[2]) - ((pldIn->rgba + 4*
vlo)[2])) + ((pldIn->rgba + 4*vlo)[2])))), (pldOut->rgba
+ 4*mid)[3] = ((unsigned char)((((0.5f))*(((pldIn->rgba +
4*vhi)[3]) - ((pldIn->rgba + 4*vlo)[3])) + ((pldIn->rgba
+ 4*vlo)[3])))))
2295 pldIn->rgba + 4*vhi)((pldOut->rgba + 4*mid)[0] = ((unsigned char)((((0.5f))*((
(pldIn->rgba + 4*vhi)[0]) - ((pldIn->rgba + 4*vlo)[0]))
+ ((pldIn->rgba + 4*vlo)[0])))), (pldOut->rgba + 4*mid
)[1] = ((unsigned char)((((0.5f))*(((pldIn->rgba + 4*vhi)[
1]) - ((pldIn->rgba + 4*vlo)[1])) + ((pldIn->rgba + 4*vlo
)[1])))), (pldOut->rgba + 4*mid)[2] = ((unsigned char)((((
0.5f))*(((pldIn->rgba + 4*vhi)[2]) - ((pldIn->rgba + 4*
vlo)[2])) + ((pldIn->rgba + 4*vlo)[2])))), (pldOut->rgba
+ 4*mid)[3] = ((unsigned char)((((0.5f))*(((pldIn->rgba +
4*vhi)[3]) - ((pldIn->rgba + 4*vlo)[3])) + ((pldIn->rgba
+ 4*vlo)[3])))));
2296 }
2297 if ((1 << limnPolyDataInfoNorm) & bitflag) {
2298 float tmp;
2299 ELL_3V_LERP(pldOut->norm + 3*mid, 0.5f,((pldOut->norm + 3*mid)[0] = (((0.5f))*(((pldIn->norm +
3*vhi)[0]) - ((pldIn->norm + 3*vlo)[0])) + ((pldIn->norm
+ 3*vlo)[0])), (pldOut->norm + 3*mid)[1] = (((0.5f))*(((pldIn
->norm + 3*vhi)[1]) - ((pldIn->norm + 3*vlo)[1])) + ((pldIn
->norm + 3*vlo)[1])), (pldOut->norm + 3*mid)[2] = (((0.5f
))*(((pldIn->norm + 3*vhi)[2]) - ((pldIn->norm + 3*vlo)
[2])) + ((pldIn->norm + 3*vlo)[2])))
2300 pldIn->norm + 3*vlo,((pldOut->norm + 3*mid)[0] = (((0.5f))*(((pldIn->norm +
3*vhi)[0]) - ((pldIn->norm + 3*vlo)[0])) + ((pldIn->norm
+ 3*vlo)[0])), (pldOut->norm + 3*mid)[1] = (((0.5f))*(((pldIn
->norm + 3*vhi)[1]) - ((pldIn->norm + 3*vlo)[1])) + ((pldIn
->norm + 3*vlo)[1])), (pldOut->norm + 3*mid)[2] = (((0.5f
))*(((pldIn->norm + 3*vhi)[2]) - ((pldIn->norm + 3*vlo)
[2])) + ((pldIn->norm + 3*vlo)[2])))
2301 pldIn->norm + 3*vhi)((pldOut->norm + 3*mid)[0] = (((0.5f))*(((pldIn->norm +
3*vhi)[0]) - ((pldIn->norm + 3*vlo)[0])) + ((pldIn->norm
+ 3*vlo)[0])), (pldOut->norm + 3*mid)[1] = (((0.5f))*(((pldIn
->norm + 3*vhi)[1]) - ((pldIn->norm + 3*vlo)[1])) + ((pldIn
->norm + 3*vlo)[1])), (pldOut->norm + 3*mid)[2] = (((0.5f
))*(((pldIn->norm + 3*vhi)[2]) - ((pldIn->norm + 3*vlo)
[2])) + ((pldIn->norm + 3*vlo)[2])));
2302 ELL_3V_NORM_TT(pldOut->norm + 3*mid, float,(tmp = ((float)((sqrt((((pldOut->norm + 3*mid))[0]*((pldOut
->norm + 3*mid))[0] + ((pldOut->norm + 3*mid))[1]*((pldOut
->norm + 3*mid))[1] + ((pldOut->norm + 3*mid))[2]*((pldOut
->norm + 3*mid))[2]))))), ((pldOut->norm + 3*mid)[0] = (
(float)((1.0/tmp)*(pldOut->norm + 3*mid)[0])), (pldOut->
norm + 3*mid)[1] = ((float)((1.0/tmp)*(pldOut->norm + 3*mid
)[1])), (pldOut->norm + 3*mid)[2] = ((float)((1.0/tmp)*(pldOut
->norm + 3*mid)[2]))))
2303 pldOut->norm + 3*mid, tmp)(tmp = ((float)((sqrt((((pldOut->norm + 3*mid))[0]*((pldOut
->norm + 3*mid))[0] + ((pldOut->norm + 3*mid))[1]*((pldOut
->norm + 3*mid))[1] + ((pldOut->norm + 3*mid))[2]*((pldOut
->norm + 3*mid))[2]))))), ((pldOut->norm + 3*mid)[0] = (
(float)((1.0/tmp)*(pldOut->norm + 3*mid)[0])), (pldOut->
norm + 3*mid)[1] = ((float)((1.0/tmp)*(pldOut->norm + 3*mid
)[1])), (pldOut->norm + 3*mid)[2] = ((float)((1.0/tmp)*(pldOut
->norm + 3*mid)[2]))));
2304 }
2305 if ((1 << limnPolyDataInfoTex2) & bitflag) {
2306 ELL_2V_LERP(pldOut->tex2 + 2*mid, 0.5f,((pldOut->tex2 + 2*mid)[0] = (((0.5f))*(((pldIn->tex2 +
2*vhi)[0]) - ((pldIn->tex2 + 2*vlo)[0])) + ((pldIn->tex2
+ 2*vlo)[0])), (pldOut->tex2 + 2*mid)[1] = (((0.5f))*(((pldIn
->tex2 + 2*vhi)[1]) - ((pldIn->tex2 + 2*vlo)[1])) + ((pldIn
->tex2 + 2*vlo)[1])))
2307 pldIn->tex2 + 2*vlo,((pldOut->tex2 + 2*mid)[0] = (((0.5f))*(((pldIn->tex2 +
2*vhi)[0]) - ((pldIn->tex2 + 2*vlo)[0])) + ((pldIn->tex2
+ 2*vlo)[0])), (pldOut->tex2 + 2*mid)[1] = (((0.5f))*(((pldIn
->tex2 + 2*vhi)[1]) - ((pldIn->tex2 + 2*vlo)[1])) + ((pldIn
->tex2 + 2*vlo)[1])))
2308 pldIn->tex2 + 2*vhi)((pldOut->tex2 + 2*mid)[0] = (((0.5f))*(((pldIn->tex2 +
2*vhi)[0]) - ((pldIn->tex2 + 2*vlo)[0])) + ((pldIn->tex2
+ 2*vlo)[0])), (pldOut->tex2 + 2*mid)[1] = (((0.5f))*(((pldIn
->tex2 + 2*vhi)[1]) - ((pldIn->tex2 + 2*vlo)[1])) + ((pldIn
->tex2 + 2*vlo)[1])));
2309 }
2310 if ((1 << limnPolyDataInfoTang) & bitflag) {
2311 float tmp;
2312 ELL_3V_LERP(pldOut->tang + 3*mid, 0.5f,((pldOut->tang + 3*mid)[0] = (((0.5f))*(((pldIn->tang +
3*vhi)[0]) - ((pldIn->tang + 3*vlo)[0])) + ((pldIn->tang
+ 3*vlo)[0])), (pldOut->tang + 3*mid)[1] = (((0.5f))*(((pldIn
->tang + 3*vhi)[1]) - ((pldIn->tang + 3*vlo)[1])) + ((pldIn
->tang + 3*vlo)[1])), (pldOut->tang + 3*mid)[2] = (((0.5f
))*(((pldIn->tang + 3*vhi)[2]) - ((pldIn->tang + 3*vlo)
[2])) + ((pldIn->tang + 3*vlo)[2])))
2313 pldIn->tang + 3*vlo,((pldOut->tang + 3*mid)[0] = (((0.5f))*(((pldIn->tang +
3*vhi)[0]) - ((pldIn->tang + 3*vlo)[0])) + ((pldIn->tang
+ 3*vlo)[0])), (pldOut->tang + 3*mid)[1] = (((0.5f))*(((pldIn
->tang + 3*vhi)[1]) - ((pldIn->tang + 3*vlo)[1])) + ((pldIn
->tang + 3*vlo)[1])), (pldOut->tang + 3*mid)[2] = (((0.5f
))*(((pldIn->tang + 3*vhi)[2]) - ((pldIn->tang + 3*vlo)
[2])) + ((pldIn->tang + 3*vlo)[2])))
2314 pldIn->tang + 3*vhi)((pldOut->tang + 3*mid)[0] = (((0.5f))*(((pldIn->tang +
3*vhi)[0]) - ((pldIn->tang + 3*vlo)[0])) + ((pldIn->tang
+ 3*vlo)[0])), (pldOut->tang + 3*mid)[1] = (((0.5f))*(((pldIn
->tang + 3*vhi)[1]) - ((pldIn->tang + 3*vlo)[1])) + ((pldIn
->tang + 3*vlo)[1])), (pldOut->tang + 3*mid)[2] = (((0.5f
))*(((pldIn->tang + 3*vhi)[2]) - ((pldIn->tang + 3*vlo)
[2])) + ((pldIn->tang + 3*vlo)[2])));
2315 ELL_3V_NORM_TT(pldOut->tang + 3*mid, float,(tmp = ((float)((sqrt((((pldOut->tang + 3*mid))[0]*((pldOut
->tang + 3*mid))[0] + ((pldOut->tang + 3*mid))[1]*((pldOut
->tang + 3*mid))[1] + ((pldOut->tang + 3*mid))[2]*((pldOut
->tang + 3*mid))[2]))))), ((pldOut->tang + 3*mid)[0] = (
(float)((1.0/tmp)*(pldOut->tang + 3*mid)[0])), (pldOut->
tang + 3*mid)[1] = ((float)((1.0/tmp)*(pldOut->tang + 3*mid
)[1])), (pldOut->tang + 3*mid)[2] = ((float)((1.0/tmp)*(pldOut
->tang + 3*mid)[2]))))
2316 pldOut->tang + 3*mid, tmp)(tmp = ((float)((sqrt((((pldOut->tang + 3*mid))[0]*((pldOut
->tang + 3*mid))[0] + ((pldOut->tang + 3*mid))[1]*((pldOut
->tang + 3*mid))[1] + ((pldOut->tang + 3*mid))[2]*((pldOut
->tang + 3*mid))[2]))))), ((pldOut->tang + 3*mid)[0] = (
(float)((1.0/tmp)*(pldOut->tang + 3*mid)[0])), (pldOut->
tang + 3*mid)[1] = ((float)((1.0/tmp)*(pldOut->tang + 3*mid
)[1])), (pldOut->tang + 3*mid)[2] = ((float)((1.0/tmp)*(pldOut
->tang + 3*mid)[2]))));
2317 }
2318 }
2319 }
2320
2321 airMopOkay(mop);
2322 return 0;
2323}
2324
2325/* helper function for the limnPolyDataNeighborList below */
2326static void
2327registerNeighbor(unsigned int *nblist, size_t *len, unsigned int *maxnb,
2328 unsigned int u, unsigned int v) {
2329 unsigned int idx=nblist[u], pointer=u, depth=1;
2330 while (idx!=0) {
2331 if (nblist[idx]==v)
2332 return; /* has already been registered */
2333 pointer=idx+1;
2334 idx=nblist[pointer];
2335 depth++;
2336 }
2337 if (depth>*maxnb)
2338 *maxnb=depth;
2339 /* do the registration */
2340 nblist[pointer]=*len;
2341 nblist[*len]=v;
2342 nblist[*len+1]=0;
2343 (*len)+=2;
2344 /* now the other way around */
2345 idx=nblist[v]; pointer=v;
2346 while (idx!=0) {
2347 /* do not have to check nblist[idx]==u due to symmetry */
2348 pointer=idx+1;
2349 idx=nblist[pointer];
2350 }
2351 nblist[pointer]=*len;
2352 nblist[*len]=u;
2353 nblist[*len+1]=0;
2354 (*len)+=2;
2355}
2356
2357/* Here's the thing with all these limnPolyDataNeighbor* functions:
2358 *
2359 * *List is used for building up the information and called by all others
2360 * *Array is a great representation if all vertices have a similar number
2361 * of neighbors (in particular, no gross outliers)
2362 * The output of this is used by glyph coloring in elf.
2363 * *ArrayComp is a compressed representation that is best if vertices have
2364 * a more variable number of neighbors
2365 * The output of this is used by surface smoothing in limn.
2366 */
2367
2368/* Mallocs *nblist and fills it with a linked list that contains all neighbors
2369 * of all n vertices in the given limnPolyData. The format is as follows:
2370 * For v<n, (*nblist)[v] is an index i (i>=n) into *nblist, or 0 if vertex v
2371 * does not have any neighbors.
2372 * For an index i obtained this way, (*nblist)[i] is a neighbor of v,
2373 * (*nblist)[i+1] is the index of the next list element (or 0).
2374 * If non-NULL, *len is set to the required size of *nblist - the initial malloc
2375 * makes a conservative guess and you may want to realloc the result to *len
2376 * bytes in order to free memory that ended up unused
2377 * If non-NULL, *m is set to the maximum number of neighbors (over all vertices)
2378 * Return value is 0 upon success, -1 upon error. Biff is used for errors.
2379 */
2380int
2381limnPolyDataNeighborList(unsigned int **nblist, size_t *len,
2382 unsigned int *maxnb, const limnPolyData *pld) {
2383 static const char me[]="limnPolyDataNeighborList";
2384 unsigned int i, j, m, estimate=0, *indx;
2385 size_t last;
2386 /* estimate the maximum number of neighborhood relations */
2387 for (i=0; i<pld->primNum; i++) {
2388 switch (pld->type[i]) {
2389 case limnPrimitiveTriangles:
2390 case limnPrimitiveQuads:
2391 estimate+=pld->icnt[i]*2;
2392 break;
2393 case limnPrimitiveTriangleStrip:
2394 case limnPrimitiveTriangleFan:
2395 estimate+=4*(pld->icnt[i]-2)+2;
2396 break;
2397 case limnPrimitiveLineStrip:
2398 estimate+=2*(pld->icnt[i]-1);
2399 break;
2400 case limnPrimitiveLines:
2401 estimate+=pld->icnt[i];
2402 break;
2403 default: /* should be a noop; silently ignore it */
2404 break;
2405 }
2406 }
2407 /* allocate *nblist */
2408 *nblist = (unsigned int*) malloc(sizeof(unsigned int)*
2409 (pld->xyzwNum+2*estimate));
2410 if (*nblist==NULL((void*)0)) {
2411 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate nblist buffer", me);
2412 return -1;
2413 }
2414 /* populate the list */
2415 memset(*nblist, 0, sizeof(unsigned int)*pld->xyzwNum)__builtin___memset_chk (*nblist, 0, sizeof(unsigned int)*pld->
xyzwNum, __builtin_object_size (*nblist, 0));
2416 last=pld->xyzwNum; m=0; indx=pld->indx;
2417 for (i=0; i<pld->primNum; i++) {
2418 switch (pld->type[i]) {
2419 case limnPrimitiveTriangles:
2420 for (j=0; j<pld->icnt[i]; j+=3) { /* go through all triangles */
2421 registerNeighbor(*nblist, &last, &m, indx[j], indx[j+1]);
2422 registerNeighbor(*nblist, &last, &m, indx[j+1], indx[j+2]);
2423 registerNeighbor(*nblist, &last, &m, indx[j+2], indx[j]);
2424 }
2425 break;
2426 case limnPrimitiveTriangleStrip:
2427 if (pld->icnt[i]>0)
2428 registerNeighbor(*nblist, &last, &m, indx[0], indx[1]);
2429 for (j=0; j<pld->icnt[i]-2; j++) {
2430 registerNeighbor(*nblist, &last, &m, indx[j], indx[j+2]);
2431 registerNeighbor(*nblist, &last, &m, indx[j+1], indx[j+2]);
2432 }
2433 break;
2434 case limnPrimitiveTriangleFan:
2435 if (pld->icnt[i]>0)
2436 registerNeighbor(*nblist, &last, &m, indx[0], indx[1]);
2437 for (j=0; j<pld->icnt[i]-2; j++) {
2438 registerNeighbor(*nblist, &last, &m, indx[0], indx[j+2]);
2439 registerNeighbor(*nblist, &last, &m, indx[j+1], indx[j+2]);
2440 }
2441 break;
2442 case limnPrimitiveQuads:
2443 for (j=0; j<pld->icnt[i]; j+=4) { /* go through all quads */
2444 registerNeighbor(*nblist, &last, &m, indx[j], indx[j+1]);
2445 registerNeighbor(*nblist, &last, &m, indx[j+1], indx[j+2]);
2446 registerNeighbor(*nblist, &last, &m, indx[j+2], indx[j+3]);
2447 registerNeighbor(*nblist, &last, &m, indx[j+3], indx[j]);
2448 }
2449 break;
2450 case limnPrimitiveLineStrip:
2451 for (j=0; j<pld->icnt[i]-1; j++) {
2452 registerNeighbor(*nblist, &last, &m, indx[j], indx[j+1]);
2453 }
2454 break;
2455 case limnPrimitiveLines:
2456 for (j=0; j<pld->icnt[i]; j+=2) {
2457 registerNeighbor(*nblist, &last, &m, indx[j], indx[j+1]);
2458 }
2459 break;
2460 default: /* should be a noop; silently ignore it */
2461 break;
2462 }
2463 indx+=pld->icnt[i];
2464 }
2465 if (len!=NULL((void*)0)) *len=last*sizeof(unsigned int);
2466 if (maxnb!=NULL((void*)0)) *maxnb=m;
2467 return 0;
2468}
2469
2470/* Over the set of all n vertices in a given limnPolyData, finds the
2471 * maximum number m of neighbors. Sets *neighbors to a malloc'ed block
2472 * of (n*m) indices and lists the neighbors of vertex v at position
2473 * (v*m) of the list, padded with -1s.
2474 * *maxnb will be set to m (and is assumed to be non-NULL!)
2475 * Returns -1 and adds a biff error if there was a problem allocating memory.
2476 */
2477int
2478limnPolyDataNeighborArray(int **neighbors, unsigned int *maxnb,
2479 const limnPolyData *pld) {
2480 static const char me[]="limnPolyDataNeighborArray";
2481 unsigned int i, *nblist, m;
2482 /* get the neighbors as a linked list */
2483 if (-1==limnPolyDataNeighborList(&nblist, NULL((void*)0), &m, pld)) {
2484 return -1;
2485 }
2486 /* convert the result into an array */
2487 *neighbors = (int *) malloc(sizeof(int)*m*pld->xyzwNum);
2488 if (NULL((void*)0)==*neighbors) {
2489 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate neighbors buffer", me);
2490 free(nblist); return -1;
2491 }
2492 for (i=0; i<pld->xyzwNum; i++) {
2493 unsigned int aidx=0, lidx=nblist[i];
2494 while (lidx!=0) {
2495 (*neighbors)[m*i+aidx]=nblist[lidx];
2496 lidx=nblist[lidx+1];
2497 aidx++;
2498 }
2499 while (aidx<m) {
2500 (*neighbors)[m*i+aidx]=-1;
2501 aidx++;
2502 }
2503 }
2504 *maxnb=m;
2505 free(nblist);
2506 return 0;
2507}
2508
2509/* Returns a compressed form of the above, rather than padding with -1s.
2510 * *neighbors is malloc'ed to an array that holds the indices of all neighbors
2511 * *idx is malloc'ed to an array of length pld->xyzwNum+1;
2512 * (*idx)[i] will give you the position in *neighbors at which the neighbors of
2513 * vertex i start
2514 * Returns -1 and adds a biff error if there was a problem allocating memory.
2515 */
2516int
2517limnPolyDataNeighborArrayComp(int **neighbors, int **idx,
2518 const limnPolyData *pld) {
2519 static const char me[]="limnPolyDataNeighborArrayComp";
2520 unsigned int i, ct, *nblist;
2521 size_t len;
2522 airArray *mop;
2523 mop = airMopNew();
2524 /* get the neighbors as a linked list */
2525 if (-1==limnPolyDataNeighborList(&nblist, &len, NULL((void*)0), pld)) {
2526 return -1;
2527 }
2528 airMopAdd(mop, nblist, airFree, airMopAlways);
2529 /* convert the result into compressed form */
2530 *neighbors = (int *) malloc(sizeof(int)*(len-pld->xyzwNum)/2);
2531 if (NULL((void*)0)==*neighbors) {
2532 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate neighbors buffer", me);
2533 airMopError(mop); return -1;
2534 }
2535 airMopAdd(mop, neighbors, airFree, airMopOnError);
2536 *idx = (int *) malloc(sizeof(int)*(pld->xyzwNum+1));
2537 if (NULL((void*)0)==*idx) {
2538 biffAddf(LIMNlimnBiffKey, "%s: couldn't allocate index buffer", me);
2539 airMopError(mop); return -1;
2540 }
2541 airMopAdd(mop, idx, airFree, airMopOnError);
2542 for (ct=i=0; i<pld->xyzwNum; i++) {
2543 unsigned int lidx=nblist[i];
2544 (*idx)[i]=ct;
2545 while (lidx!=0) {
2546 (*neighbors)[ct]=nblist[lidx];
2547 lidx=nblist[lidx+1];
2548 ct++;
2549 }
2550 }
2551 (*idx)[pld->xyzwNum]=ct;
2552 airMopOkay(mop);
2553 return 0;
2554}