/Users/scm/git/github/teem/src/limn/polymod.c

Bug Summary

File:	src/limn/polymod.c
Location:	line 1791, column 21
Description:	Result of 'calloc' is converted to a pointer of type 'unsigned char', which is incompatible with sizeof operand type 'char'

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	*/
34	static unsigned int
35	flipListIntx(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	*/
65	static void
66	flipNeighborsGet(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	*/
111	static int
112	flipNeed(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 + 3triIdx)((vert)[0] = (vertWithTri + 3triIdx)[0], (vert)[1] = (vertWithTri + 3triIdx)[1], (vert)[2] = (vertWithTri + 3triIdx)[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	*/
142	static unsigned int
143	neighborsCheckPush(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	*/
236	static unsigned int
237	maxTrianglePerPrimitive(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	*/
250	static int
251	triangleWithVertex(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	*/
328	static int
329	vertexWithTriangle(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
369	static int
370	splitListExtract(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 + 5listLen)((edgeTmp)[0]=(edgeData + 5listLen)[0], (edgeTmp)[1]=(edgeData + 5listLen)[1], (edgeTmp)[2]=(edgeData + 5listLen)[2], (edgeTmp )[3]=(edgeData + 5listLen)[3], (edgeTmp)[4]=(edgeData + 5listLen )[4]);
425	ELL_5V_COPY(edgeData + 5listLen, edgeLine)((edgeData + 5listLen)[0]=(edgeLine)[0], (edgeData + 5listLen )[1]=(edgeLine)[1], (edgeData + 5listLen)[2]=(edgeLine)[2], ( edgeData + 5listLen)[3]=(edgeLine)[3], (edgeData + 5listLen )[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	*/
452	static unsigned int
453	sweepVertNext(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	*/
469	static int
470	sweepHave2(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	*/
481	static unsigned int
482	sweepTriNext(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	*/
515	static unsigned int
516	splitTriSweep(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	*/
591	static int
592	splitTriTrack(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
766	static int
767	splitVertDup(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, oldVertNum4sizeof(float))__builtin___memcpy_chk (pld->xyzw, pldTmp.xyzw, oldVertNum 4sizeof(float), __builtin_object_size (pld->xyzw, 0));
835	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
836	memcpy(pld->rgba, pldTmp.rgba, oldVertNum4sizeof(unsigned char))__builtin___memcpy_chk (pld->rgba, pldTmp.rgba, oldVertNum 4sizeof(unsigned char), __builtin_object_size (pld->rgba , 0));
837	}
838	if ((1 << limnPolyDataInfoNorm) & bitflag) {
839	memcpy(pld->norm, pldTmp.norm, oldVertNum3sizeof(float))__builtin___memcpy_chk (pld->norm, pldTmp.norm, oldVertNum 3sizeof(float), __builtin_object_size (pld->norm, 0));
840	}
841	if ((1 << limnPolyDataInfoTex2) & bitflag) {
842	memcpy(pld->tex2, pldTmp.tex2, oldVertNum2sizeof(float))__builtin___memcpy_chk (pld->tex2, pldTmp.tex2, oldVertNum 2sizeof(float), __builtin_object_size (pld->tex2, 0));
843	}
844	if ((1 << limnPolyDataInfoTang) & bitflag) {
845	memcpy(pld->tang, pldTmp.tang, oldVertNum3sizeof(float))__builtin___memcpy_chk (pld->tang, pldTmp.tang, oldVertNum 3sizeof(float), __builtin_object_size (pld->tang, 0));
846	}
847
848	vixLut = AIR_CAST(unsigned int , calloc(2vixLutLen,((unsigned int )(calloc(2vixLutLen, sizeof(unsigned int))))
849	sizeof(unsigned int)))((unsigned int )(calloc(2vixLutLen, sizeof(unsigned int))));
850	airMopAdd(mop, vixLut, airFree, airMopAlways);
851	if (looping) {
852	for (ii=0; ii<vixLutLen; ii++) {
853	vixLut[0 + 2ii] = edgeData[2 + 5ii];
854	vixLut[1 + 2*ii] = oldVertNum + ii;
855	}
856	} else {
857	for (ii=0; ii<vixLutLen-1; ii++) {
858	vixLut[0 + 2ii] = edgeData[2 + 5ii];
859	vixLut[1 + 2*ii] = oldVertNum + ii;
860	}
861	/* now ii == vixLutLen-1 == listLen */
862	vixLut[0 + 2ii] = 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 + 4vixLut[1 + 2ii],((pld->xyzw + 4vixLut[1 + 2ii])[0] = (pld->xyzw + 4vixLut [0 + 2ii])[0], (pld->xyzw + 4vixLut[1 + 2ii])[1] = (pld ->xyzw + 4vixLut[0 + 2ii])[1], (pld->xyzw + 4vixLut[ 1 + 2ii])[2] = (pld->xyzw + 4vixLut[0 + 2ii])[2], (pld-> xyzw + 4vixLut[1 + 2ii])[3] = (pld->xyzw + 4vixLut[0 + 2 ii])[3])
869	pld->xyzw + 4vixLut[0 + 2ii])((pld->xyzw + 4vixLut[1 + 2ii])[0] = (pld->xyzw + 4vixLut [0 + 2ii])[0], (pld->xyzw + 4vixLut[1 + 2ii])[1] = (pld ->xyzw + 4vixLut[0 + 2ii])[1], (pld->xyzw + 4vixLut[ 1 + 2ii])[2] = (pld->xyzw + 4vixLut[0 + 2ii])[2], (pld-> xyzw + 4vixLut[1 + 2ii])[3] = (pld->xyzw + 4vixLut[0 + 2 ii])[3]);
870	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
871	ELL_4V_COPY(pld->rgba + 4vixLut[1 + 2ii],((pld->rgba + 4vixLut[1 + 2ii])[0] = (pld->rgba + 4vixLut [0 + 2ii])[0], (pld->rgba + 4vixLut[1 + 2ii])[1] = (pld ->rgba + 4vixLut[0 + 2ii])[1], (pld->rgba + 4vixLut[ 1 + 2ii])[2] = (pld->rgba + 4vixLut[0 + 2ii])[2], (pld-> rgba + 4vixLut[1 + 2ii])[3] = (pld->rgba + 4vixLut[0 + 2 ii])[3])
872	pld->rgba + 4vixLut[0 + 2ii])((pld->rgba + 4vixLut[1 + 2ii])[0] = (pld->rgba + 4vixLut [0 + 2ii])[0], (pld->rgba + 4vixLut[1 + 2ii])[1] = (pld ->rgba + 4vixLut[0 + 2ii])[1], (pld->rgba + 4vixLut[ 1 + 2ii])[2] = (pld->rgba + 4vixLut[0 + 2ii])[2], (pld-> rgba + 4vixLut[1 + 2ii])[3] = (pld->rgba + 4vixLut[0 + 2 ii])[3]);
873	}
874	if ((1 << limnPolyDataInfoNorm) & bitflag) {
875	ELL_3V_COPY(pld->norm + 3vixLut[1 + 2ii],((pld->norm + 3vixLut[1 + 2ii])[0] = (pld->norm + 3vixLut [0 + 2ii])[0], (pld->norm + 3vixLut[1 + 2ii])[1] = (pld ->norm + 3vixLut[0 + 2ii])[1], (pld->norm + 3vixLut[ 1 + 2ii])[2] = (pld->norm + 3vixLut[0 + 2ii])[2])
876	pld->norm + 3vixLut[0 + 2ii])((pld->norm + 3vixLut[1 + 2ii])[0] = (pld->norm + 3vixLut [0 + 2ii])[0], (pld->norm + 3vixLut[1 + 2ii])[1] = (pld ->norm + 3vixLut[0 + 2ii])[1], (pld->norm + 3vixLut[ 1 + 2ii])[2] = (pld->norm + 3vixLut[0 + 2ii])[2]);
877	}
878	if ((1 << limnPolyDataInfoTex2) & bitflag) {
879	ELL_2V_COPY(pld->tex2 + 2vixLut[1 + 2ii],((pld->tex2 + 2vixLut[1 + 2ii])[0] = (pld->tex2 + 2vixLut [0 + 2ii])[0], (pld->tex2 + 2vixLut[1 + 2ii])[1] = (pld ->tex2 + 2vixLut[0 + 2ii])[1])
880	pld->tex2 + 2vixLut[0 + 2ii])((pld->tex2 + 2vixLut[1 + 2ii])[0] = (pld->tex2 + 2vixLut [0 + 2ii])[0], (pld->tex2 + 2vixLut[1 + 2ii])[1] = (pld ->tex2 + 2vixLut[0 + 2ii])[1]);
881	}
882	if ((1 << limnPolyDataInfoTang) & bitflag) {
883	ELL_3V_COPY(pld->tang + 3vixLut[1 + 2ii],((pld->tang + 3vixLut[1 + 2ii])[0] = (pld->tang + 3vixLut [0 + 2ii])[0], (pld->tang + 3vixLut[1 + 2ii])[1] = (pld ->tang + 3vixLut[0 + 2ii])[1], (pld->tang + 3vixLut[ 1 + 2ii])[2] = (pld->tang + 3vixLut[0 + 2ii])[2])
884	pld->tang + 3vixLut[0 + 2ii])((pld->tang + 3vixLut[1 + 2ii])[0] = (pld->tang + 3vixLut [0 + 2ii])[0], (pld->tang + 3vixLut[1 + 2ii])[1] = (pld ->tang + 3vixLut[0 + 2ii])[1], (pld->tang + 3vixLut[ 1 + 2ii])[2] = (pld->tang + 3vixLut[0 + 2ii])[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	*/
922	static int
923	doSplitting(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(maxTriPerVertedgeArr->len,((unsigned int )(calloc(maxTriPerVertedgeArr->len, sizeof (unsigned int))))
948	sizeof(unsigned int)))((unsigned int )(calloc(maxTriPerVertedgeArr->len, sizeof (unsigned int))));
949	track1 = AIR_CAST(unsigned int , calloc(maxTriPerVertedgeArr->len,((unsigned int )(calloc(maxTriPerVertedgeArr->len, sizeof (unsigned int))))
950	sizeof(unsigned int)))((unsigned int )(calloc(maxTriPerVertedgeArr->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
1082	int
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), 5sizeof(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	*/
1264	int
1265	limnPolyDataVertexWindingFix(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
1283	int
1284	limnPolyDataCCFind(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 + 3realTriIdx)((baseIndx)[0] = (indxOld + 3realTriIdx)[0], (baseIndx)[1] = (indxOld + 3realTriIdx)[1], (baseIndx)[2] = (indxOld + 3realTriIdx )[2]);
1402	baseIndx += 3;
1403	pld->icnt[primIdxNew] += 3;
1404	}
1405	}
1406	}
1407
1408	airMopOkay(mop);
1409	return 0;
1410	}
1411
1412	int
1413	limnPolyDataPrimitiveSort(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)) {
1424	biffAddf(LIMNlimnBiffKey, "%s: got NULL pointer", me);
1425	return 1;
1426	}
1427	if (!(1 == _nval->dim
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);
1445	ntwo[0] = nval;
1446	ntwo[1] = nval;
1447	if (!E) E \|= nrrdJoin(nrec, ntwo, 2, 0, AIR_TRUE1);
1448	if (E) {
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++) {
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) )))
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)) {
1468	biffAddf(LIMNlimnBiffKey, "%s: couldn't allocated temp buffers", me);
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 + 2primIdx])((unsigned int)(rec[1 + 2primIdx]));
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
1500	int
1501	limnPolyDataVertexWindingFlip(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
1529	int
1530	limnPolyDataPrimitiveSelect(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 + 4newVertIdx, pldIn->xyzw + 4oldVertIdx)((pldOut->xyzw + 4newVertIdx)[0] = (pldIn->xyzw + 4oldVertIdx )[0], (pldOut->xyzw + 4newVertIdx)[1] = (pldIn->xyzw + 4oldVertIdx)[1], (pldOut->xyzw + 4newVertIdx)[2] = (pldIn ->xyzw + 4oldVertIdx)[2], (pldOut->xyzw + 4newVertIdx )[3] = (pldIn->xyzw + 4oldVertIdx)[3]);
1639	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1640	ELL_4V_COPY(pldOut->rgba + 4newVertIdx, pldIn->rgba + 4oldVertIdx)((pldOut->rgba + 4newVertIdx)[0] = (pldIn->rgba + 4oldVertIdx )[0], (pldOut->rgba + 4newVertIdx)[1] = (pldIn->rgba + 4oldVertIdx)[1], (pldOut->rgba + 4newVertIdx)[2] = (pldIn ->rgba + 4oldVertIdx)[2], (pldOut->rgba + 4newVertIdx )[3] = (pldIn->rgba + 4oldVertIdx)[3]);
1641	}
1642	if ((1 << limnPolyDataInfoNorm) & bitflag) {
1643	ELL_3V_COPY(pldOut->norm + 3newVertIdx, pldIn->norm + 3oldVertIdx)((pldOut->norm + 3newVertIdx)[0] = (pldIn->norm + 3oldVertIdx )[0], (pldOut->norm + 3newVertIdx)[1] = (pldIn->norm + 3oldVertIdx)[1], (pldOut->norm + 3newVertIdx)[2] = (pldIn ->norm + 3oldVertIdx)[2]);
1644	}
1645	if ((1 << limnPolyDataInfoTex2) & bitflag) {
1646	ELL_3V_COPY(pldOut->tex2 + 2newVertIdx, pldIn->tex2 + 2oldVertIdx)((pldOut->tex2 + 2newVertIdx)[0] = (pldIn->tex2 + 2oldVertIdx )[0], (pldOut->tex2 + 2newVertIdx)[1] = (pldIn->tex2 + 2oldVertIdx)[1], (pldOut->tex2 + 2newVertIdx)[2] = (pldIn ->tex2 + 2oldVertIdx)[2]);
1647	}
1648	if ((1 << limnPolyDataInfoTang) & bitflag) {
1649	ELL_3V_COPY(pldOut->tang + 3newVertIdx, pldIn->tang + 3oldVertIdx)((pldOut->tang + 3newVertIdx)[0] = (pldIn->tang + 3oldVertIdx )[0], (pldOut->tang + 3newVertIdx)[1] = (pldIn->tang + 3oldVertIdx)[1], (pldOut->tang + 3newVertIdx)[2] = (pldIn ->tang + 3oldVertIdx)[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	*/
1662	static int
1663	clipEdge(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+4q, float, alpha, pld->xyzw+4disc, pld->xyzw+4kept)((newpld->xyzw+4q)[0] = ((float)((((alpha))(((pld->xyzw +4kept)[0]) - ((pld->xyzw+4disc)[0])) + ((pld->xyzw+4 disc)[0])))), (newpld->xyzw+4q)[1] = ((float)((((alpha)) (((pld->xyzw+4kept)[1]) - ((pld->xyzw+4disc)[1])) + ( (pld->xyzw+4disc)[1])))), (newpld->xyzw+4q)[2] = ((float )((((alpha))(((pld->xyzw+4kept)[2]) - ((pld->xyzw+4disc )[2])) + ((pld->xyzw+4disc)[2])))), (newpld->xyzw+4q) [3] = ((float)((((alpha))(((pld->xyzw+4kept)[3]) - ((pld ->xyzw+4disc)[3])) + ((pld->xyzw+4*disc)[3])))));
1691	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1692	airArrayLenIncr(rgbaArr, 1);
1693	ELL_4V_LERP_TT(newpld->rgba+4q, unsigned char, alpha, pld->rgba+4disc, pld->rgba+4kept)((newpld->rgba+4q)[0] = ((unsigned char)((((alpha))(((pld ->rgba+4kept)[0]) - ((pld->rgba+4disc)[0])) + ((pld-> rgba+4disc)[0])))), (newpld->rgba+4q)[1] = ((unsigned char )((((alpha))(((pld->rgba+4kept)[1]) - ((pld->rgba+4disc )[1])) + ((pld->rgba+4disc)[1])))), (newpld->rgba+4q) [2] = ((unsigned char)((((alpha))(((pld->rgba+4kept)[2]) - ((pld->rgba+4disc)[2])) + ((pld->rgba+4disc)[2]))) ), (newpld->rgba+4q)[3] = ((unsigned char)((((alpha))((( pld->rgba+4kept)[3]) - ((pld->rgba+4disc)[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+3disc, pld->norm+3kept)((pld->norm+3disc)[0](pld->norm+3kept)[0] + (pld-> norm+3disc)[1](pld->norm+3kept)[1] + (pld->norm+3disc )[2](pld->norm+3*kept)[2])<0) {
1700	ELL_3V_SCALE_TT(fnorm, float, -1.0, pld->norm+3kept)((fnorm)[0] = ((float)((-1.0)(pld->norm+3kept)[0])), (fnorm )[1] = ((float)((-1.0)(pld->norm+3kept)[1])), (fnorm)[2] = ((float)((-1.0)(pld->norm+3*kept)[2])));
1701	} else {
1702	ELL_3V_COPY(fnorm, pld->norm+3kept)((fnorm)[0] = (pld->norm+3kept)[0], (fnorm)[1] = (pld-> norm+3kept)[1], (fnorm)[2] = (pld->norm+3kept)[2]);
1703	}
1704	airArrayLenIncr(normArr, 1);
1705	ELL_3V_LERP_TT(newpld->norm+3q, float, alpha, pld->norm+3disc, fnorm)((newpld->norm+3q)[0] = ((float)((((alpha))(((fnorm)[0]) - ((pld->norm+3disc)[0])) + ((pld->norm+3disc)[0]))) ), (newpld->norm+3q)[1] = ((float)((((alpha))(((fnorm)[1 ]) - ((pld->norm+3disc)[1])) + ((pld->norm+3disc)[1]) ))), (newpld->norm+3q)[2] = ((float)((((alpha))(((fnorm) [2]) - ((pld->norm+3disc)[2])) + ((pld->norm+3disc)[2 ])))));
1706	/* re-normalize */
1707	len=ELL_3V_LEN(newpld->norm+3q)(sqrt((((newpld->norm+3q))[0]((newpld->norm+3q))[0] + ((newpld->norm+3q))[1]((newpld->norm+3q))[1] + ((newpld ->norm+3q))[2]((newpld->norm+3q))[2])));
1708	if (len>1e-20) {
1709	ELL_3V_SCALE_TT(newpld->norm+3q, float, 1.0/len, newpld->norm+3q)((newpld->norm+3q)[0] = ((float)((1.0/len)(newpld->norm +3q)[0])), (newpld->norm+3q)[1] = ((float)((1.0/len)(newpld ->norm+3q)[1])), (newpld->norm+3q)[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+2q, float, alpha, pld->tex2+2disc, pld->tex2+2kept)((newpld->tex2+2q)[0] = ((float)((((alpha))(((pld->tex2 +2kept)[0]) - ((pld->tex2+2disc)[0])) + ((pld->tex2+2 disc)[0])))), (newpld->tex2+2q)[1] = ((float)((((alpha)) (((pld->tex2+2kept)[1]) - ((pld->tex2+2disc)[1])) + ( (pld->tex2+2*disc)[1])))));
1715	}
1716	if ((1 << limnPolyDataInfoTang) & bitflag) {
1717	airArrayLenIncr(tangArr, 1);
1718	ELL_3V_LERP_TT(newpld->tang+3q, float, alpha, pld->tang+3disc, pld->tang+3kept)((newpld->tang+3q)[0] = ((float)((((alpha))(((pld->tang +3kept)[0]) - ((pld->tang+3disc)[0])) + ((pld->tang+3 disc)[0])))), (newpld->tang+3q)[1] = ((float)((((alpha)) (((pld->tang+3kept)[1]) - ((pld->tang+3disc)[1])) + ( (pld->tang+3disc)[1])))), (newpld->tang+3q)[2] = ((float )((((alpha))(((pld->tang+3kept)[2]) - ((pld->tang+3disc )[2])) + ((pld->tang+3disc)[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	*/
1743	int
1744	limnPolyDataClipMulti(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))))
	Result of 'calloc' is converted to a pointer of type 'unsigned char', which is incompatible with sizeof operand type '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->xyzwNumsizeof(int)))((int )(malloc(pld->xyzwNumsizeof(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), 3sizeof(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+4q, pld->xyzw+4oldidx)((newpld->xyzw+4q)[0] = (pld->xyzw+4oldidx)[0], (newpld ->xyzw+4q)[1] = (pld->xyzw+4oldidx)[1], (newpld->xyzw +4q)[2] = (pld->xyzw+4oldidx)[2], (newpld->xyzw+4q)[ 3] = (pld->xyzw+4oldidx)[3]);
1934	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
1935	airArrayLenIncr(rgbaArr, 1);
1936	ELL_4V_COPY(newpld->rgba+4q, pld->rgba+4oldidx)((newpld->rgba+4q)[0] = (pld->rgba+4oldidx)[0], (newpld ->rgba+4q)[1] = (pld->rgba+4oldidx)[1], (newpld->rgba +4q)[2] = (pld->rgba+4oldidx)[2], (newpld->rgba+4q)[ 3] = (pld->rgba+4oldidx)[3]);
1937	}
1938	if ((1 << limnPolyDataInfoNorm) & bitflag) {
1939	airArrayLenIncr(normArr, 1);
1940	ELL_3V_COPY(newpld->norm+3q, pld->norm+3oldidx)((newpld->norm+3q)[0] = (pld->norm+3oldidx)[0], (newpld ->norm+3q)[1] = (pld->norm+3oldidx)[1], (newpld->norm +3q)[2] = (pld->norm+3oldidx)[2]);
1941	}
1942	if ((1 << limnPolyDataInfoTex2) & bitflag) {
1943	airArrayLenIncr(tex2Arr, 1);
1944	ELL_2V_COPY(newpld->tex2+2q, pld->tex2+2oldidx)((newpld->tex2+2q)[0] = (pld->tex2+2oldidx)[0], (newpld ->tex2+2q)[1] = (pld->tex2+2oldidx)[1]);
1945	}
1946	if ((1 << limnPolyDataInfoTang) & bitflag) {
1947	airArrayLenIncr(tangArr, 1);
1948	ELL_3V_COPY(newpld->tang+3q, pld->tang+3oldidx)((newpld->tang+3q)[0] = (pld->tang+3oldidx)[0], (newpld ->tang+3q)[1] = (pld->tang+3oldidx)[1], (newpld->tang +3q)[2] = (pld->tang+3oldidx)[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	*/
2026	int
2027	limnPolyDataClip(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	*/
2036	limnPolyData *
2037	limnPolyDataCompress(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+4vertMap[i], pld->xyzw+4i)((ret->xyzw+4vertMap[i])[0] = (pld->xyzw+4i)[0], (ret ->xyzw+4vertMap[i])[1] = (pld->xyzw+4i)[1], (ret-> xyzw+4vertMap[i])[2] = (pld->xyzw+4i)[2], (ret->xyzw+ 4vertMap[i])[3] = (pld->xyzw+4i)[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+4vertMap[i], pld->rgba+4i)((ret->rgba+4vertMap[i])[0] = (pld->rgba+4i)[0], (ret ->rgba+4vertMap[i])[1] = (pld->rgba+4i)[1], (ret-> rgba+4vertMap[i])[2] = (pld->rgba+4i)[2], (ret->rgba+ 4vertMap[i])[3] = (pld->rgba+4i)[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+3vertMap[i], pld->norm+3i)((ret->norm+3vertMap[i])[0] = (pld->norm+3i)[0], (ret ->norm+3vertMap[i])[1] = (pld->norm+3i)[1], (ret-> norm+3vertMap[i])[2] = (pld->norm+3i)[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+2vertMap[i], pld->tex2+2i)((ret->tex2+2vertMap[i])[0] = (pld->tex2+2i)[0], (ret ->tex2+2vertMap[i])[1] = (pld->tex2+2i)[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+3vertMap[i], pld->tang+3i)((ret->tang+3vertMap[i])[0] = (pld->tang+3i)[0], (ret ->tang+3vertMap[i])[1] = (pld->tang+3i)[1], (ret-> tang+3vertMap[i])[2] = (pld->tang+3i)[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	*/
2127	limnPolyData 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+4vertNum, plds[i]->xyzw,__builtin___memcpy_chk (ret->xyzw+4vertNum, plds[i]->xyzw , sizeof(float)4plds[i]->xyzwNum, __builtin_object_size ( ret->xyzw+4*vertNum, 0))
2163	sizeof(float)4plds[i]->xyzwNum)__builtin___memcpy_chk (ret->xyzw+4vertNum, plds[i]->xyzw , sizeof(float)4plds[i]->xyzwNum, __builtin_object_size ( ret->xyzw+4vertNum, 0));
2164	if (ret->rgba!=NULL((void*)0)) {
2165	memcpy(ret->rgba+4vertNum, plds[i]->rgba,__builtin___memcpy_chk (ret->rgba+4vertNum, plds[i]->rgba , sizeof(unsigned char)4plds[i]->xyzwNum, __builtin_object_size (ret->rgba+4*vertNum, 0))
2166	sizeof(unsigned char)4plds[i]->xyzwNum)__builtin___memcpy_chk (ret->rgba+4vertNum, plds[i]->rgba , sizeof(unsigned char)4plds[i]->xyzwNum, __builtin_object_size (ret->rgba+4vertNum, 0));
2167	}
2168	if (ret->norm!=NULL((void*)0)) {
2169	memcpy(ret->norm+3vertNum, plds[i]->norm,__builtin___memcpy_chk (ret->norm+3vertNum, plds[i]->norm , sizeof(float)3plds[i]->xyzwNum, __builtin_object_size ( ret->norm+3*vertNum, 0))
2170	sizeof(float)3plds[i]->xyzwNum)__builtin___memcpy_chk (ret->norm+3vertNum, plds[i]->norm , sizeof(float)3plds[i]->xyzwNum, __builtin_object_size ( ret->norm+3vertNum, 0));
2171	}
2172	if (ret->tex2!=NULL((void*)0)) {
2173	memcpy(ret->tex2+2vertNum, plds[i]->tex2,__builtin___memcpy_chk (ret->tex2+2vertNum, plds[i]->tex2 , sizeof(float)2plds[i]->xyzwNum, __builtin_object_size ( ret->tex2+2*vertNum, 0))
2174	sizeof(float)2plds[i]->xyzwNum)__builtin___memcpy_chk (ret->tex2+2vertNum, plds[i]->tex2 , sizeof(float)2plds[i]->xyzwNum, __builtin_object_size ( ret->tex2+2vertNum, 0));
2175	}
2176	if (ret->tang!=NULL((void*)0)) {
2177	memcpy(ret->tang+3vertNum, plds[i]->tang,__builtin___memcpy_chk (ret->tang+3vertNum, plds[i]->tang , sizeof(float)3plds[i]->xyzwNum, __builtin_object_size ( ret->tang+3*vertNum, 0))
2178	sizeof(float)3plds[i]->xyzwNum)__builtin___memcpy_chk (ret->tang+3vertNum, plds[i]->tang , sizeof(float)3plds[i]->xyzwNum, __builtin_object_size ( ret->tang+3vertNum, 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
2196	int
2197	limnPolyDataEdgeHalve(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 + nvoldvhi] = newvert[vhi + nvoldvlo] = 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 + nvoldvhi] = newvert[vhi + nvoldvlo] = 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 + nvoldvhi] = newvert[vhi + nvoldvlo] = nvidx++;
2241	}
2242	}
2243
2244	/* allocate output */
2245	bitflag = limnPolyDataInfoBitFlag(pldIn);
2246	if (limnPolyDataAlloc(pldOut, bitflag, nvidx, 34trinum, 1)) {
2247	biffAddf(LIMNlimnBiffKey, "%s: trouble allocating output", me);
2248	airMopError(mop); return 1;
2249	}
2250	pldOut->type[0] = limnPrimitiveTriangles;
2251	pldOut->icnt[0] = 34trinum;
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 + 4triidx), aa, ab, ac)((pldOut->indx + 3(0 + 4triidx))[0] = (aa), (pldOut-> indx + 3(0 + 4triidx))[1] = (ab), (pldOut->indx + 3(0 + 4triidx))[2] = (ac));
2263	ELL_3V_SET(pldOut->indx + 3(1 + 4triidx), ab, bc, ac)((pldOut->indx + 3(1 + 4triidx))[0] = (ab), (pldOut-> indx + 3(1 + 4triidx))[1] = (bc), (pldOut->indx + 3(1 + 4triidx))[2] = (ac));
2264	ELL_3V_SET(pldOut->indx + 3(2 + 4triidx), ab, bb, bc)((pldOut->indx + 3(2 + 4triidx))[0] = (ab), (pldOut-> indx + 3(2 + 4triidx))[1] = (bb), (pldOut->indx + 3(2 + 4triidx))[2] = (bc));
2265	ELL_3V_SET(pldOut->indx + 3(3 + 4triidx), ac, bc, cc)((pldOut->indx + 3(3 + 4triidx))[0] = (ac), (pldOut-> indx + 3(3 + 4triidx))[1] = (bc), (pldOut->indx + 3(3 + 4triidx))[2] = (cc));
2266	}
2267
2268	/* set output vertex info */
2269	for (vlo=0; vlo<nvold; vlo++) {
2270	ELL_4V_COPY(pldOut->xyzw + 4vlo, pldIn->xyzw + 4vlo)((pldOut->xyzw + 4vlo)[0] = (pldIn->xyzw + 4vlo)[0], ( pldOut->xyzw + 4vlo)[1] = (pldIn->xyzw + 4vlo)[1], (pldOut ->xyzw + 4vlo)[2] = (pldIn->xyzw + 4vlo)[2], (pldOut-> xyzw + 4vlo)[3] = (pldIn->xyzw + 4vlo)[3]);
2271	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
2272	ELL_4V_COPY(pldOut->rgba + 4vlo, pldIn->rgba + 4vlo)((pldOut->rgba + 4vlo)[0] = (pldIn->rgba + 4vlo)[0], ( pldOut->rgba + 4vlo)[1] = (pldIn->rgba + 4vlo)[1], (pldOut ->rgba + 4vlo)[2] = (pldIn->rgba + 4vlo)[2], (pldOut-> rgba + 4vlo)[3] = (pldIn->rgba + 4vlo)[3]);
2273	}
2274	if ((1 << limnPolyDataInfoNorm) & bitflag) {
2275	ELL_3V_COPY(pldOut->norm + 3vlo, pldIn->norm + 3vlo)((pldOut->norm + 3vlo)[0] = (pldIn->norm + 3vlo)[0], ( pldOut->norm + 3vlo)[1] = (pldIn->norm + 3vlo)[1], (pldOut ->norm + 3vlo)[2] = (pldIn->norm + 3vlo)[2]);
2276	}
2277	if ((1 << limnPolyDataInfoTex2) & bitflag) {
2278	ELL_2V_COPY(pldOut->tex2 + 2vlo, pldIn->tex2 + 2vlo)((pldOut->tex2 + 2vlo)[0] = (pldIn->tex2 + 2vlo)[0], ( pldOut->tex2 + 2vlo)[1] = (pldIn->tex2 + 2vlo)[1]);
2279	}
2280	if ((1 << limnPolyDataInfoTang) & bitflag) {
2281	ELL_3V_COPY(pldOut->tang + 3vlo, pldIn->tang + 3vlo)((pldOut->tang + 3vlo)[0] = (pldIn->tang + 3vlo)[0], ( pldOut->tang + 3vlo)[1] = (pldIn->tang + 3vlo)[1], (pldOut ->tang + 3vlo)[2] = (pldIn->tang + 3vlo)[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 + 4mid, 0.5f,((pldOut->xyzw + 4mid)[0] = (((0.5f))(((pldIn->xyzw + 4vhi)[0]) - ((pldIn->xyzw + 4vlo)[0])) + ((pldIn->xyzw + 4vlo)[0])), (pldOut->xyzw + 4mid)[1] = (((0.5f))(((pldIn ->xyzw + 4vhi)[1]) - ((pldIn->xyzw + 4vlo)[1])) + ((pldIn ->xyzw + 4vlo)[1])), (pldOut->xyzw + 4mid)[2] = (((0.5f ))(((pldIn->xyzw + 4vhi)[2]) - ((pldIn->xyzw + 4vlo) [2])) + ((pldIn->xyzw + 4vlo)[2])), (pldOut->xyzw + 4* mid)[3] = (((0.5f))(((pldIn->xyzw + 4vhi)[3]) - ((pldIn-> xyzw + 4vlo)[3])) + ((pldIn->xyzw + 4vlo)[3])))
2290	pldIn->xyzw + 4vlo,((pldOut->xyzw + 4mid)[0] = (((0.5f))(((pldIn->xyzw + 4vhi)[0]) - ((pldIn->xyzw + 4vlo)[0])) + ((pldIn->xyzw + 4vlo)[0])), (pldOut->xyzw + 4mid)[1] = (((0.5f))(((pldIn ->xyzw + 4vhi)[1]) - ((pldIn->xyzw + 4vlo)[1])) + ((pldIn ->xyzw + 4vlo)[1])), (pldOut->xyzw + 4mid)[2] = (((0.5f ))(((pldIn->xyzw + 4vhi)[2]) - ((pldIn->xyzw + 4vlo) [2])) + ((pldIn->xyzw + 4vlo)[2])), (pldOut->xyzw + 4* mid)[3] = (((0.5f))(((pldIn->xyzw + 4vhi)[3]) - ((pldIn-> xyzw + 4vlo)[3])) + ((pldIn->xyzw + 4vlo)[3])))
2291	pldIn->xyzw + 4vhi)((pldOut->xyzw + 4mid)[0] = (((0.5f))(((pldIn->xyzw + 4vhi)[0]) - ((pldIn->xyzw + 4vlo)[0])) + ((pldIn->xyzw + 4vlo)[0])), (pldOut->xyzw + 4mid)[1] = (((0.5f))(((pldIn ->xyzw + 4vhi)[1]) - ((pldIn->xyzw + 4vlo)[1])) + ((pldIn ->xyzw + 4vlo)[1])), (pldOut->xyzw + 4mid)[2] = (((0.5f ))(((pldIn->xyzw + 4vhi)[2]) - ((pldIn->xyzw + 4vlo) [2])) + ((pldIn->xyzw + 4vlo)[2])), (pldOut->xyzw + 4* mid)[3] = (((0.5f))(((pldIn->xyzw + 4vhi)[3]) - ((pldIn-> xyzw + 4vlo)[3])) + ((pldIn->xyzw + 4vlo)[3])));
2292	if ((1 << limnPolyDataInfoRGBA) & bitflag) {
2293	ELL_4V_LERP_TT(pldOut->rgba + 4mid, unsigned char, 0.5f,((pldOut->rgba + 4mid)[0] = ((unsigned char)((((0.5f))(( (pldIn->rgba + 4vhi)[0]) - ((pldIn->rgba + 4vlo)[0])) + ((pldIn->rgba + 4vlo)[0])))), (pldOut->rgba + 4mid )[1] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[ 1]) - ((pldIn->rgba + 4vlo)[1])) + ((pldIn->rgba + 4vlo )[1])))), (pldOut->rgba + 4mid)[2] = ((unsigned char)(((( 0.5f))(((pldIn->rgba + 4vhi)[2]) - ((pldIn->rgba + 4* vlo)[2])) + ((pldIn->rgba + 4vlo)[2])))), (pldOut->rgba + 4mid)[3] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[3]) - ((pldIn->rgba + 4vlo)[3])) + ((pldIn->rgba + 4vlo)[3])))))
2294	pldIn->rgba + 4vlo,((pldOut->rgba + 4mid)[0] = ((unsigned char)((((0.5f))(( (pldIn->rgba + 4vhi)[0]) - ((pldIn->rgba + 4vlo)[0])) + ((pldIn->rgba + 4vlo)[0])))), (pldOut->rgba + 4mid )[1] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[ 1]) - ((pldIn->rgba + 4vlo)[1])) + ((pldIn->rgba + 4vlo )[1])))), (pldOut->rgba + 4mid)[2] = ((unsigned char)(((( 0.5f))(((pldIn->rgba + 4vhi)[2]) - ((pldIn->rgba + 4* vlo)[2])) + ((pldIn->rgba + 4vlo)[2])))), (pldOut->rgba + 4mid)[3] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[3]) - ((pldIn->rgba + 4vlo)[3])) + ((pldIn->rgba + 4vlo)[3])))))
2295	pldIn->rgba + 4vhi)((pldOut->rgba + 4mid)[0] = ((unsigned char)((((0.5f))(( (pldIn->rgba + 4vhi)[0]) - ((pldIn->rgba + 4vlo)[0])) + ((pldIn->rgba + 4vlo)[0])))), (pldOut->rgba + 4mid )[1] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[ 1]) - ((pldIn->rgba + 4vlo)[1])) + ((pldIn->rgba + 4vlo )[1])))), (pldOut->rgba + 4mid)[2] = ((unsigned char)(((( 0.5f))(((pldIn->rgba + 4vhi)[2]) - ((pldIn->rgba + 4* vlo)[2])) + ((pldIn->rgba + 4vlo)[2])))), (pldOut->rgba + 4mid)[3] = ((unsigned char)((((0.5f))(((pldIn->rgba + 4vhi)[3]) - ((pldIn->rgba + 4vlo)[3])) + ((pldIn->rgba + 4vlo)[3])))));
2296	}
2297	if ((1 << limnPolyDataInfoNorm) & bitflag) {
2298	float tmp;
2299	ELL_3V_LERP(pldOut->norm + 3mid, 0.5f,((pldOut->norm + 3mid)[0] = (((0.5f))(((pldIn->norm + 3vhi)[0]) - ((pldIn->norm + 3vlo)[0])) + ((pldIn->norm + 3vlo)[0])), (pldOut->norm + 3mid)[1] = (((0.5f))(((pldIn ->norm + 3vhi)[1]) - ((pldIn->norm + 3vlo)[1])) + ((pldIn ->norm + 3vlo)[1])), (pldOut->norm + 3mid)[2] = (((0.5f ))(((pldIn->norm + 3vhi)[2]) - ((pldIn->norm + 3vlo) [2])) + ((pldIn->norm + 3vlo)[2])))
2300	pldIn->norm + 3vlo,((pldOut->norm + 3mid)[0] = (((0.5f))(((pldIn->norm + 3vhi)[0]) - ((pldIn->norm + 3vlo)[0])) + ((pldIn->norm + 3vlo)[0])), (pldOut->norm + 3mid)[1] = (((0.5f))(((pldIn ->norm + 3vhi)[1]) - ((pldIn->norm + 3vlo)[1])) + ((pldIn ->norm + 3vlo)[1])), (pldOut->norm + 3mid)[2] = (((0.5f ))(((pldIn->norm + 3vhi)[2]) - ((pldIn->norm + 3vlo) [2])) + ((pldIn->norm + 3vlo)[2])))
2301	pldIn->norm + 3vhi)((pldOut->norm + 3mid)[0] = (((0.5f))(((pldIn->norm + 3vhi)[0]) - ((pldIn->norm + 3vlo)[0])) + ((pldIn->norm + 3vlo)[0])), (pldOut->norm + 3mid)[1] = (((0.5f))(((pldIn ->norm + 3vhi)[1]) - ((pldIn->norm + 3vlo)[1])) + ((pldIn ->norm + 3vlo)[1])), (pldOut->norm + 3mid)[2] = (((0.5f ))(((pldIn->norm + 3vhi)[2]) - ((pldIn->norm + 3vlo) [2])) + ((pldIn->norm + 3vlo)[2])));
2302	ELL_3V_NORM_TT(pldOut->norm + 3mid, float,(tmp = ((float)((sqrt((((pldOut->norm + 3mid))[0]((pldOut ->norm + 3mid))[0] + ((pldOut->norm + 3mid))[1]((pldOut ->norm + 3mid))[1] + ((pldOut->norm + 3mid))[2]((pldOut ->norm + 3mid))[2]))))), ((pldOut->norm + 3mid)[0] = ( (float)((1.0/tmp)(pldOut->norm + 3mid)[0])), (pldOut-> norm + 3mid)[1] = ((float)((1.0/tmp)(pldOut->norm + 3mid )[1])), (pldOut->norm + 3mid)[2] = ((float)((1.0/tmp)(pldOut ->norm + 3*mid)[2]))))
2303	pldOut->norm + 3mid, tmp)(tmp = ((float)((sqrt((((pldOut->norm + 3mid))[0]((pldOut ->norm + 3mid))[0] + ((pldOut->norm + 3mid))[1]((pldOut ->norm + 3mid))[1] + ((pldOut->norm + 3mid))[2]((pldOut ->norm + 3mid))[2]))))), ((pldOut->norm + 3mid)[0] = ( (float)((1.0/tmp)(pldOut->norm + 3mid)[0])), (pldOut-> norm + 3mid)[1] = ((float)((1.0/tmp)(pldOut->norm + 3mid )[1])), (pldOut->norm + 3mid)[2] = ((float)((1.0/tmp)(pldOut ->norm + 3*mid)[2]))));
2304	}
2305	if ((1 << limnPolyDataInfoTex2) & bitflag) {
2306	ELL_2V_LERP(pldOut->tex2 + 2mid, 0.5f,((pldOut->tex2 + 2mid)[0] = (((0.5f))(((pldIn->tex2 + 2vhi)[0]) - ((pldIn->tex2 + 2vlo)[0])) + ((pldIn->tex2 + 2vlo)[0])), (pldOut->tex2 + 2mid)[1] = (((0.5f))(((pldIn ->tex2 + 2vhi)[1]) - ((pldIn->tex2 + 2vlo)[1])) + ((pldIn ->tex2 + 2*vlo)[1])))
2307	pldIn->tex2 + 2vlo,((pldOut->tex2 + 2mid)[0] = (((0.5f))(((pldIn->tex2 + 2vhi)[0]) - ((pldIn->tex2 + 2vlo)[0])) + ((pldIn->tex2 + 2vlo)[0])), (pldOut->tex2 + 2mid)[1] = (((0.5f))(((pldIn ->tex2 + 2vhi)[1]) - ((pldIn->tex2 + 2vlo)[1])) + ((pldIn ->tex2 + 2*vlo)[1])))
2308	pldIn->tex2 + 2vhi)((pldOut->tex2 + 2mid)[0] = (((0.5f))(((pldIn->tex2 + 2vhi)[0]) - ((pldIn->tex2 + 2vlo)[0])) + ((pldIn->tex2 + 2vlo)[0])), (pldOut->tex2 + 2mid)[1] = (((0.5f))(((pldIn ->tex2 + 2vhi)[1]) - ((pldIn->tex2 + 2vlo)[1])) + ((pldIn ->tex2 + 2*vlo)[1])));
2309	}
2310	if ((1 << limnPolyDataInfoTang) & bitflag) {
2311	float tmp;
2312	ELL_3V_LERP(pldOut->tang + 3mid, 0.5f,((pldOut->tang + 3mid)[0] = (((0.5f))(((pldIn->tang + 3vhi)[0]) - ((pldIn->tang + 3vlo)[0])) + ((pldIn->tang + 3vlo)[0])), (pldOut->tang + 3mid)[1] = (((0.5f))(((pldIn ->tang + 3vhi)[1]) - ((pldIn->tang + 3vlo)[1])) + ((pldIn ->tang + 3vlo)[1])), (pldOut->tang + 3mid)[2] = (((0.5f ))(((pldIn->tang + 3vhi)[2]) - ((pldIn->tang + 3vlo) [2])) + ((pldIn->tang + 3vlo)[2])))
2313	pldIn->tang + 3vlo,((pldOut->tang + 3mid)[0] = (((0.5f))(((pldIn->tang + 3vhi)[0]) - ((pldIn->tang + 3vlo)[0])) + ((pldIn->tang + 3vlo)[0])), (pldOut->tang + 3mid)[1] = (((0.5f))(((pldIn ->tang + 3vhi)[1]) - ((pldIn->tang + 3vlo)[1])) + ((pldIn ->tang + 3vlo)[1])), (pldOut->tang + 3mid)[2] = (((0.5f ))(((pldIn->tang + 3vhi)[2]) - ((pldIn->tang + 3vlo) [2])) + ((pldIn->tang + 3vlo)[2])))
2314	pldIn->tang + 3vhi)((pldOut->tang + 3mid)[0] = (((0.5f))(((pldIn->tang + 3vhi)[0]) - ((pldIn->tang + 3vlo)[0])) + ((pldIn->tang + 3vlo)[0])), (pldOut->tang + 3mid)[1] = (((0.5f))(((pldIn ->tang + 3vhi)[1]) - ((pldIn->tang + 3vlo)[1])) + ((pldIn ->tang + 3vlo)[1])), (pldOut->tang + 3mid)[2] = (((0.5f ))(((pldIn->tang + 3vhi)[2]) - ((pldIn->tang + 3vlo) [2])) + ((pldIn->tang + 3vlo)[2])));
2315	ELL_3V_NORM_TT(pldOut->tang + 3mid, float,(tmp = ((float)((sqrt((((pldOut->tang + 3mid))[0]((pldOut ->tang + 3mid))[0] + ((pldOut->tang + 3mid))[1]((pldOut ->tang + 3mid))[1] + ((pldOut->tang + 3mid))[2]((pldOut ->tang + 3mid))[2]))))), ((pldOut->tang + 3mid)[0] = ( (float)((1.0/tmp)(pldOut->tang + 3mid)[0])), (pldOut-> tang + 3mid)[1] = ((float)((1.0/tmp)(pldOut->tang + 3mid )[1])), (pldOut->tang + 3mid)[2] = ((float)((1.0/tmp)(pldOut ->tang + 3*mid)[2]))))
2316	pldOut->tang + 3mid, tmp)(tmp = ((float)((sqrt((((pldOut->tang + 3mid))[0]((pldOut ->tang + 3mid))[0] + ((pldOut->tang + 3mid))[1]((pldOut ->tang + 3mid))[1] + ((pldOut->tang + 3mid))[2]((pldOut ->tang + 3mid))[2]))))), ((pldOut->tang + 3mid)[0] = ( (float)((1.0/tmp)(pldOut->tang + 3mid)[0])), (pldOut-> tang + 3mid)[1] = ((float)((1.0/tmp)(pldOut->tang + 3mid )[1])), (pldOut->tang + 3mid)[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 */
2326	static void
2327	registerNeighbor(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	*/
2380	int
2381	limnPolyDataNeighborList(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	*/
2477	int
2478	limnPolyDataNeighborArray(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)mpld->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)[mi+aidx]=nblist[lidx];
2496	lidx=nblist[lidx+1];
2497	aidx++;
2498	}
2499	while (aidx<m) {
2500	(neighbors)[mi+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	*/
2516	int
2517	limnPolyDataNeighborArrayComp(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	}