Bug Summary

File:src/seek/extract.c
Location:line 903, column 5
Description:Potential leak of memory pointed to by 'bag'

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
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22*/
23
24#include "seek.h"
25#include "privateSeek.h"
26
27static baggage *
28baggageNew(seekContext *sctx) {
29 baggage *bag;
30 unsigned int sx;
31
32 bag = AIR_CALLOC(1, baggage)(baggage*)(calloc((1), sizeof(baggage)));
7
←
Within the expansion of the macro 'AIR_CALLOC':
→
a
Memory is allocated
33
34 /* this is basically the mapping from the 12 edges on each voxel to
35 the 5 unique edges for each sample on the slab, based on the lay-out
36 defined in the beginning of tables.c */
37 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
38 /* X Y */
39 bag->evti[ 0] = 0 + 5*(0 + sx*0);
40 bag->evti[ 1] = 1 + 5*(0 + sx*0);
41 bag->evti[ 2] = 1 + 5*(1 + sx*0);
42 bag->evti[ 3] = 0 + 5*(0 + sx*1);
43 bag->evti[ 4] = 2 + 5*(0 + sx*0);
44 bag->evti[ 5] = 2 + 5*(1 + sx*0);
45 bag->evti[ 6] = 2 + 5*(0 + sx*1);
46 bag->evti[ 7] = 2 + 5*(1 + sx*1);
47 bag->evti[ 8] = 3 + 5*(0 + sx*0);
48 bag->evti[ 9] = 4 + 5*(0 + sx*0);
49 bag->evti[10] = 4 + 5*(1 + sx*0);
50 bag->evti[11] = 3 + 5*(0 + sx*1);
51
52 switch (sctx->type) {
8
←
Control jumps to 'case seekTypeRidgeSurface:' at line 53
→
53 case seekTypeRidgeSurface:
54 case seekTypeRidgeSurfaceOP:
55 case seekTypeRidgeSurfaceT:
56 bag->esIdx = 2;
57 bag->modeSign = -1;
58 break;
9
←
 Execution continues on line 84
→
59 case seekTypeValleySurface:
60 case seekTypeValleySurfaceOP:
61 case seekTypeValleySurfaceT:
62 bag->esIdx = 0;
63 bag->modeSign = +1;
64 break;
65 case seekTypeMaximalSurface:
66 bag->esIdx = 0;
67 bag->modeSign = -1;
68 break;
69 case seekTypeMinimalSurface:
70 bag->esIdx = 2;
71 bag->modeSign = +1;
72 break;
73 default:
74 /* biffAddf(SEEK, "%s: feature type %s not handled", me,
75 airEnumStr(seekType, sctx->type));
76 return 1;
77 */
78 /* without biff, we get as nasty as possible */
79 bag->esIdx = UINT_MAX(2147483647 *2U +1U);
80 bag->modeSign = 0;
81 break;
82 }
83
84 if (seekTypeIsocontour == sctx->type) {
10
←
Taking false branch
→
85 if (sctx->ninscl) {
86 bag->scllup = nrrdDLookup[sctx->ninscl->type];
87 bag->scldata = sctx->ninscl->data;
88 } else {
89 bag->scllup = nrrdDLookup[sctx->nsclDerived->type];
90 bag->scldata = sctx->nsclDerived->data;
91 }
92 } else {
93 bag->scllup = NULL((void*)0);
94 bag->scldata = NULL((void*)0);
95 }
96
97 bag->xyzwArr = NULL((void*)0);
98 bag->normArr = NULL((void*)0);
99 bag->indxArr = NULL((void*)0);
100 return bag;
101}
102
103static baggage *
104baggageNix(baggage *bag) {
105
106 if (bag) {
107 airArrayNix(bag->normArr);
108 airArrayNix(bag->xyzwArr);
109 airArrayNix(bag->indxArr);
110
111 airFree(bag);
112 }
113 return NULL((void*)0);
114}
115
116static int
117outputInit(seekContext *sctx, baggage *bag, limnPolyData *lpld) {
118 static const char me[]="outputInit";
119 unsigned int estVertNum, estFaceNum, minI, maxI, valI, *spanHist;
120 airPtrPtrUnion appu;
121 int E;
122
123 if (seekTypeIsocontour == sctx->type
124 && AIR_IN_OP(sctx->range->min, sctx->isovalue, sctx->range->max)((sctx->range->min) < (sctx->isovalue) &&
(sctx->isovalue) < (sctx->range->max))) {
125 unsigned int estVoxNum=0;
126 /* estimate number of voxels, faces, and vertices involved */
127 spanHist = AIR_CAST(unsigned int*, sctx->nspanHist->data)((unsigned int*)(sctx->nspanHist->data));
128 valI = airIndex(sctx->range->min, sctx->isovalue, sctx->range->max,
129 sctx->spanSize);
130 for (minI=0; minI<=valI; minI++) {
131 for (maxI=valI; maxI<sctx->spanSize; maxI++) {
132 estVoxNum += spanHist[minI + sctx->spanSize*maxI];
133 }
134 }
135 estVertNum = AIR_CAST(unsigned int, estVoxNum*(sctx->vertsPerVoxel))((unsigned int)(estVoxNum*(sctx->vertsPerVoxel)));
136 estFaceNum = AIR_CAST(unsigned int, estVoxNum*(sctx->facesPerVoxel))((unsigned int)(estVoxNum*(sctx->facesPerVoxel)));
137 if (sctx->verbose) {
138 fprintf(stderr__stderrp, "%s: estimated vox --> vert, face: %u --> %u, %u\n", me,
139 estVoxNum, estVertNum, estFaceNum);
140 }
141 } else {
142 estVertNum = 0;
143 estFaceNum = 0;
144 }
145 /* need something non-zero so that pre-allocations below aren't no-ops */
146 estVertNum = AIR_MAX(1, estVertNum)((1) > (estVertNum) ? (1) : (estVertNum));
147 estFaceNum = AIR_MAX(1, estFaceNum)((1) > (estFaceNum) ? (1) : (estFaceNum));
148
149 /* initialize limnPolyData with estimated # faces and vertices */
150 /* we will manage the innards of the limnPolyData entirely ourselves */
151 if (limnPolyDataAlloc(lpld, 0, 0, 0, 0)) {
152 biffAddf(SEEKseekBiffKey, "%s: trouble emptying given polydata", me);
153 return 1;
154 }
155 bag->xyzwArr = airArrayNew((appu.f = &(lpld->xyzw), appu.v),
156 &(lpld->xyzwNum),
157 4*sizeof(float), sctx->pldArrIncr);
158 if (sctx->normalsFind) {
159 bag->normArr = airArrayNew((appu.f = &(lpld->norm), appu.v),
160 &(lpld->normNum),
161 3*sizeof(float), sctx->pldArrIncr);
162 } else {
163 bag->normArr = NULL((void*)0);
164 }
165 bag->indxArr = airArrayNew((appu.ui = &(lpld->indx), appu.v),
166 &(lpld->indxNum),
167 sizeof(unsigned int), sctx->pldArrIncr);
168 lpld->primNum = 1; /* for now, its just triangle soup */
169 lpld->type = AIR_CALLOC(lpld->primNum, unsigned char)(unsigned char*)(calloc((lpld->primNum), sizeof(unsigned char
)));
170 lpld->icnt = AIR_CALLOC(lpld->primNum, unsigned int)(unsigned int*)(calloc((lpld->primNum), sizeof(unsigned int
)));
171 lpld->type[0] = limnPrimitiveTriangles;
172 lpld->icnt[0] = 0; /* incremented below */
173
174 E = 0;
175 airArrayLenPreSet(bag->xyzwArr, estVertNum);
176 E |= !(bag->xyzwArr->data);
177 if (sctx->normalsFind) {
178 airArrayLenPreSet(bag->normArr, estVertNum);
179 E |= !(bag->normArr->data);
180 }
181 airArrayLenPreSet(bag->indxArr, 3*estFaceNum);
182 E |= !(bag->indxArr->data);
183 if (E) {
184 biffAddf(SEEKseekBiffKey, "%s: couldn't pre-allocate contour geometry (%p %p %p)", me,
185 bag->xyzwArr->data,
186 (sctx->normalsFind ? bag->normArr->data : NULL((void*)0)),
187 bag->indxArr->data);
188 return 1;
189 }
190
191 /* initialize output summary info */
192 sctx->voxNum = 0;
193 sctx->vertNum = 0;
194 sctx->faceNum = 0;
195
196 return 0;
197}
198
199static double
200sclGet(seekContext *sctx, baggage *bag,
201 unsigned int xi, unsigned int yi, unsigned int zi) {
202
203 zi = AIR_MIN(sctx->sz-1, zi)((sctx->sz-1) < (zi) ? (sctx->sz-1) : (zi));
204 return bag->scllup(bag->scldata, xi + sctx->sx*(yi + sctx->sy*zi));
205}
206
207void
208_seekIdxProbe(seekContext *sctx, baggage *bag,
209 double xi, double yi, double zi) {
210 double idxOut[4], idxIn[4];
211 AIR_UNUSED(bag)(void)(bag);
212
213 ELL_4V_SET(idxOut, xi, yi, zi, 1)((idxOut)[0] = (xi), (idxOut)[1] = (yi), (idxOut)[2] = (zi), (
idxOut)[3] = (1));
214 ELL_4MV_MUL(idxIn, sctx->txfIdx, idxOut)((idxIn)[0]=(sctx->txfIdx)[ 0]*(idxOut)[0]+(sctx->txfIdx
)[ 1]*(idxOut)[1]+(sctx->txfIdx)[ 2]*(idxOut)[2]+(sctx->
txfIdx)[ 3]*(idxOut)[3], (idxIn)[1]=(sctx->txfIdx)[ 4]*(idxOut
)[0]+(sctx->txfIdx)[ 5]*(idxOut)[1]+(sctx->txfIdx)[ 6]*
(idxOut)[2]+(sctx->txfIdx)[ 7]*(idxOut)[3], (idxIn)[2]=(sctx
->txfIdx)[ 8]*(idxOut)[0]+(sctx->txfIdx)[ 9]*(idxOut)[1
]+(sctx->txfIdx)[10]*(idxOut)[2]+(sctx->txfIdx)[11]*(idxOut
)[3], (idxIn)[3]=(sctx->txfIdx)[12]*(idxOut)[0]+(sctx->
txfIdx)[13]*(idxOut)[1]+(sctx->txfIdx)[14]*(idxOut)[2]+(sctx
->txfIdx)[15]*(idxOut)[3]);
215 ELL_4V_HOMOG(idxIn, idxIn)((idxIn)[0] = (idxIn)[0]/(idxIn)[3], (idxIn)[1] = (idxIn)[1]/
(idxIn)[3], (idxIn)[2] = (idxIn)[2]/(idxIn)[3], (idxIn)[3] = 1.0
);
216 gageProbe(sctx->gctx, idxIn[0], idxIn[1], idxIn[2]);
217 return;
218}
219
220/*
221** this is one of the few things that has to operate on more than one
222** zi plane at once, and it is honestly probably the primary motivation
223** for putting zi into the baggage.
224**
225** NOTE: this is doing some bounds (on the positive x, y, z edges of the
226** volume) that probably should be done closer to the caller
227*/
228static int
229evecFlipProbe(seekContext *sctx, baggage *bag,
230 signed char *flip, /* OUTPUT HERE */
231 unsigned int xi, unsigned int yi, unsigned int ziOff,
232 unsigned int dx, unsigned int dy, unsigned int dz) {
233 static const char me[]="evecFlipProbe";
234 unsigned int sx, sy, sz;
235 double u, du, dot, wantDot, minDu, mode;
236 double current[3], next[3], posNext[3], posA[3], posB[3], evecA[3], evecB[3];
237
238 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
239 sy = AIR_CAST(unsigned int, sctx->sy)((unsigned int)(sctx->sy));
240 sz = AIR_CAST(unsigned int, sctx->sz)((unsigned int)(sctx->sz));
241
242 if (!(xi + dx < sx
243 && yi + dy < sy
244 && bag->zi + ziOff + dz < sz)) {
245 /* the edge we're being asked about is outside the volume */
246 *flip = 0;
247 return 0;
248 }
249
250 /* Note: Strength checking is no longer performed here.
251 * TS 2009-08-18 */
252
253 /* this edge is in bounds */
254
255 ELL_3V_SET(posA, xi, yi, bag->zi+ziOff)((posA)[0] = (xi), (posA)[1] = (yi), (posA)[2] = (bag->zi+
ziOff));
256 ELL_3V_SET(posB, xi+dx, yi+dy, bag->zi+ziOff+dz)((posB)[0] = (xi+dx), (posB)[1] = (yi+dy), (posB)[2] = (bag->
zi+ziOff+dz));
257 ELL_3V_COPY(evecA, sctx->evec((evecA)[0] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff + 2
*(xi + sx*yi))))[0], (evecA)[1] = (sctx->evec + 3*(bag->
esIdx + 3*(ziOff + 2*(xi + sx*yi))))[1], (evecA)[2] = (sctx->
evec + 3*(bag->esIdx + 3*(ziOff + 2*(xi + sx*yi))))[2])
258 + 3*(bag->esIdx + 3*(ziOff + 2*(xi + sx*yi))))((evecA)[0] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff + 2
*(xi + sx*yi))))[0], (evecA)[1] = (sctx->evec + 3*(bag->
esIdx + 3*(ziOff + 2*(xi + sx*yi))))[1], (evecA)[2] = (sctx->
evec + 3*(bag->esIdx + 3*(ziOff + 2*(xi + sx*yi))))[2]);
259 ELL_3V_COPY(evecB, sctx->evec((evecB)[0] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff+dz
+ 2*(xi+dx + sx*(yi+dy)))))[0], (evecB)[1] = (sctx->evec +
3*(bag->esIdx + 3*(ziOff+dz + 2*(xi+dx + sx*(yi+dy)))))[1
], (evecB)[2] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff+
dz + 2*(xi+dx + sx*(yi+dy)))))[2])
260 + 3*(bag->esIdx + 3*(ziOff+dz + 2*(xi+dx + sx*(yi+dy)))))((evecB)[0] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff+dz
+ 2*(xi+dx + sx*(yi+dy)))))[0], (evecB)[1] = (sctx->evec +
3*(bag->esIdx + 3*(ziOff+dz + 2*(xi+dx + sx*(yi+dy)))))[1
], (evecB)[2] = (sctx->evec + 3*(bag->esIdx + 3*(ziOff+
dz + 2*(xi+dx + sx*(yi+dy)))))[2]);
261
262#define SETNEXT(uu) \
263 ELL_3V_SCALE_ADD2(posNext, 1.0-(uu), posA, (uu), posB)((posNext)[0] = (1.0-(uu))*(posA)[0] + ((uu))*(posB)[0], (posNext
)[1] = (1.0-(uu))*(posA)[1] + ((uu))*(posB)[1], (posNext)[2] =
(1.0-(uu))*(posA)[2] + ((uu))*(posB)[2]); \
264 _seekIdxProbe(sctx, bag, posNext[0], posNext[1], posNext[2]); \
265 ELL_3V_COPY(next, sctx->evecAns + 3*bag->esIdx)((next)[0] = (sctx->evecAns + 3*bag->esIdx)[0], (next)[
1] = (sctx->evecAns + 3*bag->esIdx)[1], (next)[2] = (sctx
->evecAns + 3*bag->esIdx)[2]); \
266 if (ELL_3V_DOT(current, next)((current)[0]*(next)[0] + (current)[1]*(next)[1] + (current)[
2]*(next)[2]) < 0) { \
267 ELL_3V_SCALE(next, -1, next)((next)[0] = (-1)*(next)[0], (next)[1] = (-1)*(next)[1], (next
)[2] = (-1)*(next)[2]); \
268 } \
269 dot = ELL_3V_DOT(current, next)((current)[0]*(next)[0] + (current)[1]*(next)[1] + (current)[
2]*(next)[2]); \
270 mode = bag->modeSign*airMode3_d(sctx->evalAns);
271
272 ELL_3V_COPY(current, evecA)((current)[0] = (evecA)[0], (current)[1] = (evecA)[1], (current
)[2] = (evecA)[2]);
273 u = 0;
274 du = 0.49999;
275 wantDot = 0.9; /* between cos(pi/6) and cos(pi/8) */
276 minDu = 0.0002;
277 while (u + du < 1.0) {
278 SETNEXT(u+du);
279 /* Note: This was set to -0.8 in the original code. Again, I found
280 * that increasing it could eliminate spurious holes in the
281 * mesh. TS 2009-08-18 */
282 if (mode < -0.99) {
283 /* sorry, eigenvalue mode got too close to 2nd order isotropy */
284 *flip = 0;
285 return 0;
286 }
287 while (dot < wantDot) {
288 /* angle between current and next is too big; reduce step */
289 du /= 2;
290 if (du < minDu) {
291 fprintf(stderr__stderrp, "%s: evector wild @ u=%g: du=%g < minDu=%g; "
292 "dot=%g; mode = %g; "
293 "(xi,yi,zi)=(%u,%u,%u+%u); (dx,dy,dz)=(%u,%u,%u) ",
294 me, u, du, minDu,
295 dot, mode,
296 xi, yi, bag->zi, ziOff, dx, dy, dz);
297 *flip = 0;
298 return 0;
299 }
300 SETNEXT(u+du);
301 if (mode < -0.99) {
302 /* sorry, eigenvalue mode got too close to 2nd order isotropy */
303 *flip = 0;
304 return 0;
305 }
306 }
307 /* current and next have a small angle between them */
308 ELL_3V_COPY(current, next)((current)[0] = (next)[0], (current)[1] = (next)[1], (current
)[2] = (next)[2]);
309 u += du;
310 }
311 /* last fake iteration, to check endpoint explicitly */
312 u = 1.0;
313 SETNEXT(u);
314 if (dot < wantDot) {
315 biffAddf(SEEKseekBiffKey, "%s: confused at end of edge", me);
316 return 1;
317 }
318 ELL_3V_COPY(current, next)((current)[0] = (next)[0], (current)[1] = (next)[1], (current
)[2] = (next)[2]);
319
320#undef SETNEXT
321
322 /* we have now tracked the eigenvector along the edge */
323 dot = ELL_3V_DOT(current, evecB)((current)[0]*(evecB)[0] + (current)[1]*(evecB)[1] + (current
)[2]*(evecB)[2]);
324 *flip = (dot > 0
325 ? +1
326 : -1);
327 return 0;
328}
329
330/*
331** !!! this has to be done as a separate second pass because of how
332** !!! the flip quantities correspond to the edges between voxels
333**
334** For efficiency, evecFlipProbe is only executed on request (by
335** setting sctx->treated: 0x01 requests all edges of this voxel, 0x02
336** states that unique edge 3 was treated, 0x04 for unique edge 4) TS
337*/
338static int
339evecFlipShuffleProbe(seekContext *sctx, baggage *bag) {
340 static const char me[]="evecFlipShuffleProbe";
341 unsigned int xi, yi, sx, sy, si;
342 signed char flipA, flipB, flipC;
343
344 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
345 sy = AIR_CAST(unsigned int, sctx->sy)((unsigned int)(sctx->sy));
346
347 /* NOTE: these have to go all the way to sy-1 and sx-1, instead of
348 sy-2 and sx-2 (like shuffleProbe() below) because of the need to
349 set the flips on the edges at the positive X and Y volume
350 boundary. The necessary bounds checking happens in
351 evecFlipProbe(): messy, I know */
352 for (yi=0; yi<sy; yi++) {
353 for (xi=0; xi<sx; xi++) {
354 si = xi + sx*yi;
355 /* ================================================= */
356 if (sctx->treated[si]&0x02) { /* has been treated, just copy result */
357 sctx->flip[0 + 5*si] = sctx->flip[3 + 5*si];
358 } else if (sctx->treated[si]&0x01 ||
359 (yi!=0 && sctx->treated[xi+sx*(yi-1)]&0x01)) {
360 /* need to treat this */
361 if (evecFlipProbe(sctx, bag, &flipA, xi, yi, 0, 1, 0, 0)) {
362 biffAddf(SEEKseekBiffKey, "%s: problem at (xi,yi) = (%u,%u), zi=0", me, xi, yi);
363 return 1;
364 }
365 sctx->flip[0 + 5*si] = flipA;
366 }
367 if (sctx->treated[si]&0x04) { /* has been treated, just copy */
368 sctx->flip[1 + 5*si] = sctx->flip[4 + 5*si];
369 } else if (sctx->treated[si]&0x01 ||
370 (xi!=0 && sctx->treated[xi-1+sx*yi]&0x01)) {
371 if (evecFlipProbe(sctx, bag, &flipB, xi, yi, 0, 0, 1, 0)) {
372 biffAddf(SEEKseekBiffKey, "%s: problem at (xi,yi) = (%u,%u), zi=0", me, xi, yi);
373 return 1;
374 }
375 sctx->flip[1 + 5*si] = flipB;
376 }
377 if (sctx->treated[si]&0x01 || (xi!=0 && sctx->treated[xi-1+sx*yi]&0x01) ||
378 (yi!=0 && sctx->treated[xi+sx*(yi-1)]&0x01) ||
379 (xi!=0 && yi!=0 && sctx->treated[xi-1+sx*(yi-1)]&0x01)) {
380 if (evecFlipProbe(sctx, bag, &flipA, xi, yi, 0, 0, 0, 1)) {
381 biffAddf(SEEKseekBiffKey, "%s: problem at (xi,yi,zi) = (%u,%u,%u)", me,
382 xi, yi, bag->zi);
383 return 1;
384 }
385 sctx->flip[2 + 5*si] = flipA;
386 }
387 if (sctx->treated[si]&0x01 ||
388 (yi!=0 && sctx->treated[xi+sx*(yi-1)]&0x01)) {
389 if (evecFlipProbe(sctx, bag, &flipB, xi, yi, 1, 1, 0, 0)) {
390 biffAddf(SEEKseekBiffKey, "%s: problem at (xi,yi,zi) = (%u,%u,%u)", me,
391 xi, yi, bag->zi);
392 return 1;
393 }
394 sctx->flip[3 + 5*si] = flipB;
395 sctx->treated[si]|=0x02;
396 } else
397 sctx->treated[si]&=0xFD;
398 if (sctx->treated[si]&0x01 ||
399 (xi!=0 && sctx->treated[xi-1+sx*yi]&0x01)) {
400 if (evecFlipProbe(sctx, bag, &flipC, xi, yi, 1, 0, 1, 0)) {
401 biffAddf(SEEKseekBiffKey, "%s: problem at (xi,yi,zi) = (%u,%u,%u)", me,
402 xi, yi, bag->zi);
403 return 1;
404 }
405 sctx->flip[4 + 5*si] = flipC;
406 sctx->treated[si]|=0x04;
407 } else
408 sctx->treated[si]&=0xFB;
409 /* ================================================= */
410 }
411 }
412
413 return 0;
414}
415
416static int
417shuffleProbe(seekContext *sctx, baggage *bag) {
418 static const char me[]="shuffleProbe";
419 unsigned int xi, yi, sx, sy, si, spi;
420
421 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
422 sy = AIR_CAST(unsigned int, sctx->sy)((unsigned int)(sctx->sy));
423
424 if (!sctx->strengthUse) { /* just request all edges */
425 memset(sctx->treated, 0x01, sizeof(char)*sctx->sx*sctx->sy)__builtin___memset_chk (sctx->treated, 0x01, sizeof(char)*
sctx->sx*sctx->sy, __builtin_object_size (sctx->treated
, 0));
426 } else {
427 if (!bag->zi) {
428 /* clear full treated array */
429 memset(sctx->treated, 0, sizeof(char)*sctx->sx*sctx->sy)__builtin___memset_chk (sctx->treated, 0, sizeof(char)*sctx
->sx*sctx->sy, __builtin_object_size (sctx->treated,
0));
430 } else {
431 /* only clear requests for edge orientation */
432 for (yi=0; yi<sy; yi++) {
433 for (xi=0; xi<sx; xi++) {
434 sctx->treated[xi+sx*yi] &= 0xFE;
435 }
436 }
437 }
438 }
439
440 for (yi=0; yi<sy; yi++) {
441 for (xi=0; xi<sx; xi++) {
442 si = xi + sx*yi;
443 spi = (xi+1) + (sx+2)*(yi+1);
444 /* ================================================= */
445 if (!bag->zi) {
446 /* ----------------- set/probe bottom of initial slab */
447 sctx->vidx[0 + 5*si] = -1;
448 sctx->vidx[1 + 5*si] = -1;
449 if (sctx->gctx) { /* HEY: need this check, what's the right way? */
450 _seekIdxProbe(sctx, bag, xi, yi, 0);
451 }
452 if (sctx->strengthUse) {
453 sctx->stng[0 + 2*si] = sctx->strengthSign*sctx->stngAns[0];
454 if (!si) {
455 sctx->strengthSeenMax = sctx->stng[0 + 2*si];
456 }
457 sctx->strengthSeenMax = AIR_MAX(sctx->strengthSeenMax,((sctx->strengthSeenMax) > (sctx->stng[0 + 2*si]) ? (
sctx->strengthSeenMax) : (sctx->stng[0 + 2*si]))
458 sctx->stng[0 + 2*si])((sctx->strengthSeenMax) > (sctx->stng[0 + 2*si]) ? (
sctx->strengthSeenMax) : (sctx->stng[0 + 2*si]));
459 }
460 switch (sctx->type) {
461 case seekTypeIsocontour:
462 sctx->sclv[0 + 4*spi] = (sclGet(sctx, bag, xi, yi, 0)
463 - sctx->isovalue);
464 sctx->sclv[1 + 4*spi] = (sclGet(sctx, bag, xi, yi, 0)
465 - sctx->isovalue);
466 sctx->sclv[2 + 4*spi] = (sclGet(sctx, bag, xi, yi, 1)
467 - sctx->isovalue);
468 break;
469 case seekTypeRidgeSurface:
470 case seekTypeValleySurface:
471 case seekTypeMaximalSurface:
472 case seekTypeMinimalSurface:
473 case seekTypeRidgeSurfaceOP:
474 case seekTypeValleySurfaceOP:
475 ELL_3V_COPY(sctx->grad + 3*(0 + 2*si), sctx->gradAns)((sctx->grad + 3*(0 + 2*si))[0] = (sctx->gradAns)[0], (
sctx->grad + 3*(0 + 2*si))[1] = (sctx->gradAns)[1], (sctx
->grad + 3*(0 + 2*si))[2] = (sctx->gradAns)[2]);
476 ELL_3V_COPY(sctx->eval + 3*(0 + 2*si), sctx->evalAns)((sctx->eval + 3*(0 + 2*si))[0] = (sctx->evalAns)[0], (
sctx->eval + 3*(0 + 2*si))[1] = (sctx->evalAns)[1], (sctx
->eval + 3*(0 + 2*si))[2] = (sctx->evalAns)[2]);
477 ELL_3M_COPY(sctx->evec + 9*(0 + 2*si), sctx->evecAns)((((sctx->evec + 9*(0 + 2*si))+0)[0] = ((sctx->evecAns)
+0)[0], ((sctx->evec + 9*(0 + 2*si))+0)[1] = ((sctx->evecAns
)+0)[1], ((sctx->evec + 9*(0 + 2*si))+0)[2] = ((sctx->evecAns
)+0)[2]), (((sctx->evec + 9*(0 + 2*si))+3)[0] = ((sctx->
evecAns)+3)[0], ((sctx->evec + 9*(0 + 2*si))+3)[1] = ((sctx
->evecAns)+3)[1], ((sctx->evec + 9*(0 + 2*si))+3)[2] = (
(sctx->evecAns)+3)[2]), (((sctx->evec + 9*(0 + 2*si))+6
)[0] = ((sctx->evecAns)+6)[0], ((sctx->evec + 9*(0 + 2*
si))+6)[1] = ((sctx->evecAns)+6)[1], ((sctx->evec + 9*(
0 + 2*si))+6)[2] = ((sctx->evecAns)+6)[2]));
478 break;
479 }
480 } else {
481 /* ------------------- shuffle to bottom from top of slab */
482 sctx->vidx[0 + 5*si] = sctx->vidx[3 + 5*si];
483 sctx->vidx[1 + 5*si] = sctx->vidx[4 + 5*si];
484 if (sctx->strengthUse) {
485 sctx->stng[0 + 2*si] = sctx->stng[1 + 2*si];
486 }
487 switch (sctx->type) {
488 case seekTypeIsocontour:
489 sctx->sclv[0 + 4*spi] = sctx->sclv[1 + 4*spi];
490 sctx->sclv[1 + 4*spi] = sctx->sclv[2 + 4*spi];
491 sctx->sclv[2 + 4*spi] = sctx->sclv[3 + 4*spi];
492 break;
493 case seekTypeRidgeSurface:
494 case seekTypeValleySurface:
495 case seekTypeMaximalSurface:
496 case seekTypeMinimalSurface:
497 case seekTypeRidgeSurfaceOP:
498 case seekTypeValleySurfaceOP:
499 ELL_3V_COPY(sctx->grad + 3*(0 + 2*si), sctx->grad + 3*(1 + 2*si))((sctx->grad + 3*(0 + 2*si))[0] = (sctx->grad + 3*(1 + 2
*si))[0], (sctx->grad + 3*(0 + 2*si))[1] = (sctx->grad +
3*(1 + 2*si))[1], (sctx->grad + 3*(0 + 2*si))[2] = (sctx->
grad + 3*(1 + 2*si))[2]);
500 ELL_3V_COPY(sctx->eval + 3*(0 + 2*si), sctx->eval + 3*(1 + 2*si))((sctx->eval + 3*(0 + 2*si))[0] = (sctx->eval + 3*(1 + 2
*si))[0], (sctx->eval + 3*(0 + 2*si))[1] = (sctx->eval +
3*(1 + 2*si))[1], (sctx->eval + 3*(0 + 2*si))[2] = (sctx->
eval + 3*(1 + 2*si))[2]);
501 ELL_3M_COPY(sctx->evec + 9*(0 + 2*si), sctx->evec + 9*(1 + 2*si))((((sctx->evec + 9*(0 + 2*si))+0)[0] = ((sctx->evec + 9
*(1 + 2*si))+0)[0], ((sctx->evec + 9*(0 + 2*si))+0)[1] = (
(sctx->evec + 9*(1 + 2*si))+0)[1], ((sctx->evec + 9*(0 +
2*si))+0)[2] = ((sctx->evec + 9*(1 + 2*si))+0)[2]), (((sctx
->evec + 9*(0 + 2*si))+3)[0] = ((sctx->evec + 9*(1 + 2*
si))+3)[0], ((sctx->evec + 9*(0 + 2*si))+3)[1] = ((sctx->
evec + 9*(1 + 2*si))+3)[1], ((sctx->evec + 9*(0 + 2*si))+3
)[2] = ((sctx->evec + 9*(1 + 2*si))+3)[2]), (((sctx->evec
+ 9*(0 + 2*si))+6)[0] = ((sctx->evec + 9*(1 + 2*si))+6)[0
], ((sctx->evec + 9*(0 + 2*si))+6)[1] = ((sctx->evec + 9
*(1 + 2*si))+6)[1], ((sctx->evec + 9*(0 + 2*si))+6)[2] = (
(sctx->evec + 9*(1 + 2*si))+6)[2]));
502 break;
503 }
504 }
505 /* ----------------------- set/probe top of slab */
506 sctx->vidx[2 + 5*si] = -1;
507 sctx->vidx[3 + 5*si] = -1;
508 sctx->vidx[4 + 5*si] = -1;
509 if (sctx->gctx) { /* HEY: need this check, what's the right way? */
510 _seekIdxProbe(sctx, bag, xi, yi, bag->zi+1);
511 }
512 if (sctx->strengthUse) {
513 sctx->stng[1 + 2*si] = sctx->strengthSign*sctx->stngAns[0];
514 sctx->strengthSeenMax = AIR_MAX(sctx->strengthSeenMax,((sctx->strengthSeenMax) > (sctx->stng[1 + 2*si]) ? (
sctx->strengthSeenMax) : (sctx->stng[1 + 2*si]))
515 sctx->stng[1 + 2*si])((sctx->strengthSeenMax) > (sctx->stng[1 + 2*si]) ? (
sctx->strengthSeenMax) : (sctx->stng[1 + 2*si]));
516 if (sctx->stng[0+2*si]>sctx->strength ||
517 sctx->stng[1+2*si]>sctx->strength) {
518 /* mark up to four voxels as needed */
519 sctx->treated[si] |= 0x01;
520 if (xi!=0) sctx->treated[xi-1+sx*yi] |= 0x01;
521 if (yi!=0) sctx->treated[xi+sx*(yi-1)] |= 0x01;
522 if (xi!=0 && yi!=0) sctx->treated[xi-1+sx*(yi-1)] |= 0x01;
523 }
524 }
525 switch (sctx->type) {
526 case seekTypeIsocontour:
527 sctx->sclv[3 + 4*spi] = (sclGet(sctx, bag, xi, yi, bag->zi+2)
528 - sctx->isovalue);
529 break;
530 case seekTypeRidgeSurface:
531 case seekTypeValleySurface:
532 case seekTypeMaximalSurface:
533 case seekTypeMinimalSurface:
534 case seekTypeRidgeSurfaceOP:
535 case seekTypeValleySurfaceOP:
536 ELL_3V_COPY(sctx->grad + 3*(1 + 2*si), sctx->gradAns)((sctx->grad + 3*(1 + 2*si))[0] = (sctx->gradAns)[0], (
sctx->grad + 3*(1 + 2*si))[1] = (sctx->gradAns)[1], (sctx
->grad + 3*(1 + 2*si))[2] = (sctx->gradAns)[2]);
537 ELL_3V_COPY(sctx->eval + 3*(1 + 2*si), sctx->evalAns)((sctx->eval + 3*(1 + 2*si))[0] = (sctx->evalAns)[0], (
sctx->eval + 3*(1 + 2*si))[1] = (sctx->evalAns)[1], (sctx
->eval + 3*(1 + 2*si))[2] = (sctx->evalAns)[2]);
538 ELL_3M_COPY(sctx->evec + 9*(1 + 2*si), sctx->evecAns)((((sctx->evec + 9*(1 + 2*si))+0)[0] = ((sctx->evecAns)
+0)[0], ((sctx->evec + 9*(1 + 2*si))+0)[1] = ((sctx->evecAns
)+0)[1], ((sctx->evec + 9*(1 + 2*si))+0)[2] = ((sctx->evecAns
)+0)[2]), (((sctx->evec + 9*(1 + 2*si))+3)[0] = ((sctx->
evecAns)+3)[0], ((sctx->evec + 9*(1 + 2*si))+3)[1] = ((sctx
->evecAns)+3)[1], ((sctx->evec + 9*(1 + 2*si))+3)[2] = (
(sctx->evecAns)+3)[2]), (((sctx->evec + 9*(1 + 2*si))+6
)[0] = ((sctx->evecAns)+6)[0], ((sctx->evec + 9*(1 + 2*
si))+6)[1] = ((sctx->evecAns)+6)[1], ((sctx->evec + 9*(
1 + 2*si))+6)[2] = ((sctx->evecAns)+6)[2]));
539 break;
540 }
541 /* ================================================= */
542 }
543 /* copy ends of this scanline left/right to margin */
544 if (seekTypeIsocontour == sctx->type) {
545 ELL_4V_COPY(sctx->sclv + 4*(0 + (sx+2)*(yi+1)),((sctx->sclv + 4*(0 + (sx+2)*(yi+1)))[0] = (sctx->sclv +
4*(1 + (sx+2)*(yi+1)))[0], (sctx->sclv + 4*(0 + (sx+2)*(yi
+1)))[1] = (sctx->sclv + 4*(1 + (sx+2)*(yi+1)))[1], (sctx->
sclv + 4*(0 + (sx+2)*(yi+1)))[2] = (sctx->sclv + 4*(1 + (sx
+2)*(yi+1)))[2], (sctx->sclv + 4*(0 + (sx+2)*(yi+1)))[3] =
(sctx->sclv + 4*(1 + (sx+2)*(yi+1)))[3])
546 sctx->sclv + 4*(1 + (sx+2)*(yi+1)))((sctx->sclv + 4*(0 + (sx+2)*(yi+1)))[0] = (sctx->sclv +
4*(1 + (sx+2)*(yi+1)))[0], (sctx->sclv + 4*(0 + (sx+2)*(yi
+1)))[1] = (sctx->sclv + 4*(1 + (sx+2)*(yi+1)))[1], (sctx->
sclv + 4*(0 + (sx+2)*(yi+1)))[2] = (sctx->sclv + 4*(1 + (sx
+2)*(yi+1)))[2], (sctx->sclv + 4*(0 + (sx+2)*(yi+1)))[3] =
(sctx->sclv + 4*(1 + (sx+2)*(yi+1)))[3]);
547 ELL_4V_COPY(sctx->sclv + 4*(sx+1 + (sx+2)*(yi+1)),((sctx->sclv + 4*(sx+1 + (sx+2)*(yi+1)))[0] = (sctx->sclv
+ 4*(sx + (sx+2)*(yi+1)))[0], (sctx->sclv + 4*(sx+1 + (sx
+2)*(yi+1)))[1] = (sctx->sclv + 4*(sx + (sx+2)*(yi+1)))[1]
, (sctx->sclv + 4*(sx+1 + (sx+2)*(yi+1)))[2] = (sctx->sclv
+ 4*(sx + (sx+2)*(yi+1)))[2], (sctx->sclv + 4*(sx+1 + (sx
+2)*(yi+1)))[3] = (sctx->sclv + 4*(sx + (sx+2)*(yi+1)))[3]
)
548 sctx->sclv + 4*(sx + (sx+2)*(yi+1)))((sctx->sclv + 4*(sx+1 + (sx+2)*(yi+1)))[0] = (sctx->sclv
+ 4*(sx + (sx+2)*(yi+1)))[0], (sctx->sclv + 4*(sx+1 + (sx
+2)*(yi+1)))[1] = (sctx->sclv + 4*(sx + (sx+2)*(yi+1)))[1]
, (sctx->sclv + 4*(sx+1 + (sx+2)*(yi+1)))[2] = (sctx->sclv
+ 4*(sx + (sx+2)*(yi+1)))[2], (sctx->sclv + 4*(sx+1 + (sx
+2)*(yi+1)))[3] = (sctx->sclv + 4*(sx + (sx+2)*(yi+1)))[3]
);
549 }
550 }
551 /* copy top and bottom scanline up/down to margin */
552 if (seekTypeIsocontour == sctx->type) {
553 for (xi=0; xi<sx+2; xi++) {
554 ELL_4V_COPY(sctx->sclv + 4*(xi + (sx+2)*0),((sctx->sclv + 4*(xi + (sx+2)*0))[0] = (sctx->sclv + 4*
(xi + (sx+2)*1))[0], (sctx->sclv + 4*(xi + (sx+2)*0))[1] =
(sctx->sclv + 4*(xi + (sx+2)*1))[1], (sctx->sclv + 4*(
xi + (sx+2)*0))[2] = (sctx->sclv + 4*(xi + (sx+2)*1))[2], (
sctx->sclv + 4*(xi + (sx+2)*0))[3] = (sctx->sclv + 4*(xi
+ (sx+2)*1))[3])
555 sctx->sclv + 4*(xi + (sx+2)*1))((sctx->sclv + 4*(xi + (sx+2)*0))[0] = (sctx->sclv + 4*
(xi + (sx+2)*1))[0], (sctx->sclv + 4*(xi + (sx+2)*0))[1] =
(sctx->sclv + 4*(xi + (sx+2)*1))[1], (sctx->sclv + 4*(
xi + (sx+2)*0))[2] = (sctx->sclv + 4*(xi + (sx+2)*1))[2], (
sctx->sclv + 4*(xi + (sx+2)*0))[3] = (sctx->sclv + 4*(xi
+ (sx+2)*1))[3]);
556 ELL_4V_COPY(sctx->sclv + 4*(xi + (sx+2)*(sy+1)),((sctx->sclv + 4*(xi + (sx+2)*(sy+1)))[0] = (sctx->sclv
+ 4*(xi + (sx+2)*sy))[0], (sctx->sclv + 4*(xi + (sx+2)*(sy
+1)))[1] = (sctx->sclv + 4*(xi + (sx+2)*sy))[1], (sctx->
sclv + 4*(xi + (sx+2)*(sy+1)))[2] = (sctx->sclv + 4*(xi + (
sx+2)*sy))[2], (sctx->sclv + 4*(xi + (sx+2)*(sy+1)))[3] = (
sctx->sclv + 4*(xi + (sx+2)*sy))[3])
557 sctx->sclv + 4*(xi + (sx+2)*sy))((sctx->sclv + 4*(xi + (sx+2)*(sy+1)))[0] = (sctx->sclv
+ 4*(xi + (sx+2)*sy))[0], (sctx->sclv + 4*(xi + (sx+2)*(sy
+1)))[1] = (sctx->sclv + 4*(xi + (sx+2)*sy))[1], (sctx->
sclv + 4*(xi + (sx+2)*(sy+1)))[2] = (sctx->sclv + 4*(xi + (
sx+2)*sy))[2], (sctx->sclv + 4*(xi + (sx+2)*(sy+1)))[3] = (
sctx->sclv + 4*(xi + (sx+2)*sy))[3]);
558 }
559 }
560
561 /* this is done as a separate pass because it looks at values between
562 voxels (so its indexing is not trivial to fold into loops above) */
563 if (seekTypeRidgeSurface == sctx->type
564 || seekTypeValleySurface == sctx->type
565 || seekTypeMaximalSurface == sctx->type
566 || seekTypeMinimalSurface == sctx->type) {
567 if (evecFlipShuffleProbe(sctx, bag)) {
568 biffAddf(SEEKseekBiffKey, "%s: trouble at zi=%u\n", me, bag->zi);
569 return 1;
570 }
571 }
572
573 return 0;
574}
575
576#define VAL(xx, yy, zz) (val[4*( (xx) + (yy)*(sx+2) + spi) + (zz+1)])
577static void
578voxelGrads(double vgrad[8][3], double *val, int sx, int spi) {
579 ELL_3V_SET(vgrad[0],((vgrad[0])[0] = (VAL( 1, 0, 0) - VAL(-1, 0, 0)), (vgrad[0])[
1] = (VAL( 0, 1, 0) - VAL( 0, -1, 0)), (vgrad[0])[2] = (VAL( 0
, 0, 1) - VAL( 0, 0, -1)))
580 VAL( 1, 0, 0) - VAL(-1, 0, 0),((vgrad[0])[0] = (VAL( 1, 0, 0) - VAL(-1, 0, 0)), (vgrad[0])[
1] = (VAL( 0, 1, 0) - VAL( 0, -1, 0)), (vgrad[0])[2] = (VAL( 0
, 0, 1) - VAL( 0, 0, -1)))
581 VAL( 0, 1, 0) - VAL( 0, -1, 0),((vgrad[0])[0] = (VAL( 1, 0, 0) - VAL(-1, 0, 0)), (vgrad[0])[
1] = (VAL( 0, 1, 0) - VAL( 0, -1, 0)), (vgrad[0])[2] = (VAL( 0
, 0, 1) - VAL( 0, 0, -1)))
582 VAL( 0, 0, 1) - VAL( 0, 0, -1))((vgrad[0])[0] = (VAL( 1, 0, 0) - VAL(-1, 0, 0)), (vgrad[0])[
1] = (VAL( 0, 1, 0) - VAL( 0, -1, 0)), (vgrad[0])[2] = (VAL( 0
, 0, 1) - VAL( 0, 0, -1)));
583 ELL_3V_SET(vgrad[1],((vgrad[1])[0] = (VAL( 2, 0, 0) - VAL( 0, 0, 0)), (vgrad[1])[
1] = (VAL( 1, 1, 0) - VAL( 1, -1, 0)), (vgrad[1])[2] = (VAL( 1
, 0, 1) - VAL( 1, 0, -1)))
584 VAL( 2, 0, 0) - VAL( 0, 0, 0),((vgrad[1])[0] = (VAL( 2, 0, 0) - VAL( 0, 0, 0)), (vgrad[1])[
1] = (VAL( 1, 1, 0) - VAL( 1, -1, 0)), (vgrad[1])[2] = (VAL( 1
, 0, 1) - VAL( 1, 0, -1)))
585 VAL( 1, 1, 0) - VAL( 1, -1, 0),((vgrad[1])[0] = (VAL( 2, 0, 0) - VAL( 0, 0, 0)), (vgrad[1])[
1] = (VAL( 1, 1, 0) - VAL( 1, -1, 0)), (vgrad[1])[2] = (VAL( 1
, 0, 1) - VAL( 1, 0, -1)))
586 VAL( 1, 0, 1) - VAL( 1, 0, -1))((vgrad[1])[0] = (VAL( 2, 0, 0) - VAL( 0, 0, 0)), (vgrad[1])[
1] = (VAL( 1, 1, 0) - VAL( 1, -1, 0)), (vgrad[1])[2] = (VAL( 1
, 0, 1) - VAL( 1, 0, -1)));
587 ELL_3V_SET(vgrad[2],((vgrad[2])[0] = (VAL( 1, 1, 0) - VAL(-1, 1, 0)), (vgrad[2])[
1] = (VAL( 0, 2, 0) - VAL( 0, 0, 0)), (vgrad[2])[2] = (VAL( 0
, 1, 1) - VAL( 0, 1, -1)))
588 VAL( 1, 1, 0) - VAL(-1, 1, 0),((vgrad[2])[0] = (VAL( 1, 1, 0) - VAL(-1, 1, 0)), (vgrad[2])[
1] = (VAL( 0, 2, 0) - VAL( 0, 0, 0)), (vgrad[2])[2] = (VAL( 0
, 1, 1) - VAL( 0, 1, -1)))
589 VAL( 0, 2, 0) - VAL( 0, 0, 0),((vgrad[2])[0] = (VAL( 1, 1, 0) - VAL(-1, 1, 0)), (vgrad[2])[
1] = (VAL( 0, 2, 0) - VAL( 0, 0, 0)), (vgrad[2])[2] = (VAL( 0
, 1, 1) - VAL( 0, 1, -1)))
590 VAL( 0, 1, 1) - VAL( 0, 1, -1))((vgrad[2])[0] = (VAL( 1, 1, 0) - VAL(-1, 1, 0)), (vgrad[2])[
1] = (VAL( 0, 2, 0) - VAL( 0, 0, 0)), (vgrad[2])[2] = (VAL( 0
, 1, 1) - VAL( 0, 1, -1)));
591 ELL_3V_SET(vgrad[3],((vgrad[3])[0] = (VAL( 2, 1, 0) - VAL( 0, 1, 0)), (vgrad[3])[
1] = (VAL( 1, 2, 0) - VAL( 1, 0, 0)), (vgrad[3])[2] = (VAL( 1
, 1, 1) - VAL( 1, 1, -1)))
592 VAL( 2, 1, 0) - VAL( 0, 1, 0),((vgrad[3])[0] = (VAL( 2, 1, 0) - VAL( 0, 1, 0)), (vgrad[3])[
1] = (VAL( 1, 2, 0) - VAL( 1, 0, 0)), (vgrad[3])[2] = (VAL( 1
, 1, 1) - VAL( 1, 1, -1)))
593 VAL( 1, 2, 0) - VAL( 1, 0, 0),((vgrad[3])[0] = (VAL( 2, 1, 0) - VAL( 0, 1, 0)), (vgrad[3])[
1] = (VAL( 1, 2, 0) - VAL( 1, 0, 0)), (vgrad[3])[2] = (VAL( 1
, 1, 1) - VAL( 1, 1, -1)))
594 VAL( 1, 1, 1) - VAL( 1, 1, -1))((vgrad[3])[0] = (VAL( 2, 1, 0) - VAL( 0, 1, 0)), (vgrad[3])[
1] = (VAL( 1, 2, 0) - VAL( 1, 0, 0)), (vgrad[3])[2] = (VAL( 1
, 1, 1) - VAL( 1, 1, -1)));
595 ELL_3V_SET(vgrad[4],((vgrad[4])[0] = (VAL( 1, 0, 1) - VAL(-1, 0, 1)), (vgrad[4])[
1] = (VAL( 0, 1, 1) - VAL( 0, -1, 1)), (vgrad[4])[2] = (VAL( 0
, 0, 2) - VAL( 0, 0, 0)))
596 VAL( 1, 0, 1) - VAL(-1, 0, 1),((vgrad[4])[0] = (VAL( 1, 0, 1) - VAL(-1, 0, 1)), (vgrad[4])[
1] = (VAL( 0, 1, 1) - VAL( 0, -1, 1)), (vgrad[4])[2] = (VAL( 0
, 0, 2) - VAL( 0, 0, 0)))
597 VAL( 0, 1, 1) - VAL( 0, -1, 1),((vgrad[4])[0] = (VAL( 1, 0, 1) - VAL(-1, 0, 1)), (vgrad[4])[
1] = (VAL( 0, 1, 1) - VAL( 0, -1, 1)), (vgrad[4])[2] = (VAL( 0
, 0, 2) - VAL( 0, 0, 0)))
598 VAL( 0, 0, 2) - VAL( 0, 0, 0))((vgrad[4])[0] = (VAL( 1, 0, 1) - VAL(-1, 0, 1)), (vgrad[4])[
1] = (VAL( 0, 1, 1) - VAL( 0, -1, 1)), (vgrad[4])[2] = (VAL( 0
, 0, 2) - VAL( 0, 0, 0)));
599 ELL_3V_SET(vgrad[5],((vgrad[5])[0] = (VAL( 2, 0, 1) - VAL( 0, 0, 1)), (vgrad[5])[
1] = (VAL( 1, 1, 1) - VAL( 1, -1, 1)), (vgrad[5])[2] = (VAL( 1
, 0, 2) - VAL( 1, 0, 0)))
600 VAL( 2, 0, 1) - VAL( 0, 0, 1),((vgrad[5])[0] = (VAL( 2, 0, 1) - VAL( 0, 0, 1)), (vgrad[5])[
1] = (VAL( 1, 1, 1) - VAL( 1, -1, 1)), (vgrad[5])[2] = (VAL( 1
, 0, 2) - VAL( 1, 0, 0)))
601 VAL( 1, 1, 1) - VAL( 1, -1, 1),((vgrad[5])[0] = (VAL( 2, 0, 1) - VAL( 0, 0, 1)), (vgrad[5])[
1] = (VAL( 1, 1, 1) - VAL( 1, -1, 1)), (vgrad[5])[2] = (VAL( 1
, 0, 2) - VAL( 1, 0, 0)))
602 VAL( 1, 0, 2) - VAL( 1, 0, 0))((vgrad[5])[0] = (VAL( 2, 0, 1) - VAL( 0, 0, 1)), (vgrad[5])[
1] = (VAL( 1, 1, 1) - VAL( 1, -1, 1)), (vgrad[5])[2] = (VAL( 1
, 0, 2) - VAL( 1, 0, 0)));
603 ELL_3V_SET(vgrad[6],((vgrad[6])[0] = (VAL( 1, 1, 1) - VAL(-1, 1, 1)), (vgrad[6])[
1] = (VAL( 0, 2, 1) - VAL( 0, 0, 1)), (vgrad[6])[2] = (VAL( 0
, 1, 2) - VAL( 0, 1, 0)))
604 VAL( 1, 1, 1) - VAL(-1, 1, 1),((vgrad[6])[0] = (VAL( 1, 1, 1) - VAL(-1, 1, 1)), (vgrad[6])[
1] = (VAL( 0, 2, 1) - VAL( 0, 0, 1)), (vgrad[6])[2] = (VAL( 0
, 1, 2) - VAL( 0, 1, 0)))
605 VAL( 0, 2, 1) - VAL( 0, 0, 1),((vgrad[6])[0] = (VAL( 1, 1, 1) - VAL(-1, 1, 1)), (vgrad[6])[
1] = (VAL( 0, 2, 1) - VAL( 0, 0, 1)), (vgrad[6])[2] = (VAL( 0
, 1, 2) - VAL( 0, 1, 0)))
606 VAL( 0, 1, 2) - VAL( 0, 1, 0))((vgrad[6])[0] = (VAL( 1, 1, 1) - VAL(-1, 1, 1)), (vgrad[6])[
1] = (VAL( 0, 2, 1) - VAL( 0, 0, 1)), (vgrad[6])[2] = (VAL( 0
, 1, 2) - VAL( 0, 1, 0)));
607 ELL_3V_SET(vgrad[7],((vgrad[7])[0] = (VAL( 2, 1, 1) - VAL( 0, 1, 1)), (vgrad[7])[
1] = (VAL( 1, 2, 1) - VAL( 1, 0, 1)), (vgrad[7])[2] = (VAL( 1
, 1, 2) - VAL( 1, 1, 0)))
608 VAL( 2, 1, 1) - VAL( 0, 1, 1),((vgrad[7])[0] = (VAL( 2, 1, 1) - VAL( 0, 1, 1)), (vgrad[7])[
1] = (VAL( 1, 2, 1) - VAL( 1, 0, 1)), (vgrad[7])[2] = (VAL( 1
, 1, 2) - VAL( 1, 1, 0)))
609 VAL( 1, 2, 1) - VAL( 1, 0, 1),((vgrad[7])[0] = (VAL( 2, 1, 1) - VAL( 0, 1, 1)), (vgrad[7])[
1] = (VAL( 1, 2, 1) - VAL( 1, 0, 1)), (vgrad[7])[2] = (VAL( 1
, 1, 2) - VAL( 1, 1, 0)))
610 VAL( 1, 1, 2) - VAL( 1, 1, 0))((vgrad[7])[0] = (VAL( 2, 1, 1) - VAL( 0, 1, 1)), (vgrad[7])[
1] = (VAL( 1, 2, 1) - VAL( 1, 0, 1)), (vgrad[7])[2] = (VAL( 1
, 1, 2) - VAL( 1, 1, 0)));
611}
612#undef VAL
613
614static void
615vvalIsoSet(seekContext *sctx, baggage *bag, double vval[8],
616 unsigned int xi, unsigned int yi) {
617 unsigned int sx, si, spi, vi;
618
619 AIR_UNUSED(bag)(void)(bag);
620 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
621 si = xi + sx*yi;
622 spi = (xi+1) + (sx+2)*(yi+1);
623
624 /* learn voxel values */
625 /* X Y Z */
626 vval[0] = sctx->sclv[4*(0 + 0*(sx+2) + spi) + 1];
627 vval[1] = sctx->sclv[4*(1 + 0*(sx+2) + spi) + 1];
628 vval[2] = sctx->sclv[4*(0 + 1*(sx+2) + spi) + 1];
629 vval[3] = sctx->sclv[4*(1 + 1*(sx+2) + spi) + 1];
630 vval[4] = sctx->sclv[4*(0 + 0*(sx+2) + spi) + 2];
631 vval[5] = sctx->sclv[4*(1 + 0*(sx+2) + spi) + 2];
632 vval[6] = sctx->sclv[4*(0 + 1*(sx+2) + spi) + 2];
633 vval[7] = sctx->sclv[4*(1 + 1*(sx+2) + spi) + 2];
634 if (sctx->strengthUse) {
635 double s, w, ssum, wsum;
636 /* Z X Y */
637 ssum = wsum = 0;
638#define ACCUM(vi) w = AIR_ABS(1.0/vval[vi])((1.0/vval[vi]) > 0.0f ? (1.0/vval[vi]) : -(1.0/vval[vi])); ssum += w*s; wsum += w
639 s = sctx->stng[0 + 2*(0 + 0*sx + si)]; ACCUM(0);
640 s = sctx->stng[0 + 2*(1 + 0*sx + si)]; ACCUM(1);
641 s = sctx->stng[0 + 2*(0 + 1*sx + si)]; ACCUM(2);
642 s = sctx->stng[0 + 2*(1 + 1*sx + si)]; ACCUM(3);
643 s = sctx->stng[1 + 2*(0 + 0*sx + si)]; ACCUM(4);
644 s = sctx->stng[1 + 2*(1 + 0*sx + si)]; ACCUM(5);
645 s = sctx->stng[1 + 2*(0 + 1*sx + si)]; ACCUM(6);
646 s = sctx->stng[1 + 2*(1 + 1*sx + si)]; ACCUM(7);
647#undef ACCUM
648 if (ssum/wsum < sctx->strength) {
649 for (vi=0; vi<8; vi++) {
650 vval[vi] = 0;
651 }
652 }
653 }
654 return;
655}
656
657
658static void
659vvalSurfSet(seekContext *sctx, baggage *bag, double vval[8],
660 unsigned int xi, unsigned int yi) {
661 /* static const char me[]="vvalSurfSet"; */
662 double evec[8][3], grad[8][3], stng[8], maxStrength=0;
663 signed char flip[12]={0,0,0,0,0,0,0,0,0,0,0,0}, flipProd;
664 unsigned int sx, si, vi, ei, vrti[8];
665
666 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
667 si = xi + sx*yi;
668 vrti[0] = 0 + 2*(xi + 0 + sx*(yi + 0));
669 vrti[1] = 0 + 2*(xi + 1 + sx*(yi + 0));
670 vrti[2] = 0 + 2*(xi + 0 + sx*(yi + 1));
671 vrti[3] = 0 + 2*(xi + 1 + sx*(yi + 1));
672 vrti[4] = 1 + 2*(xi + 0 + sx*(yi + 0));
673 vrti[5] = 1 + 2*(xi + 1 + sx*(yi + 0));
674 vrti[6] = 1 + 2*(xi + 0 + sx*(yi + 1));
675 vrti[7] = 1 + 2*(xi + 1 + sx*(yi + 1));
676
677 /* Our strategy is to create all triangles of which at least some
678 * part meets the strength criterion, and to trim them in a
679 * post-process. This avoids ragged boundaries */
680 for (vi=0; vi<8; vi++) {
681 ELL_3V_COPY(grad[vi], sctx->grad + 3*vrti[vi])((grad[vi])[0] = (sctx->grad + 3*vrti[vi])[0], (grad[vi])[
1] = (sctx->grad + 3*vrti[vi])[1], (grad[vi])[2] = (sctx->
grad + 3*vrti[vi])[2]);
682 ELL_3V_COPY(evec[vi], sctx->evec + 3*(bag->esIdx + 3*vrti[vi]))((evec[vi])[0] = (sctx->evec + 3*(bag->esIdx + 3*vrti[vi
]))[0], (evec[vi])[1] = (sctx->evec + 3*(bag->esIdx + 3
*vrti[vi]))[1], (evec[vi])[2] = (sctx->evec + 3*(bag->esIdx
+ 3*vrti[vi]))[2]);
683 if (sctx->strengthUse) {
684 stng[vi] = sctx->stng[vrti[vi]];
685 if (!vi) {
686 maxStrength = stng[vi];
687 } else {
688 maxStrength = AIR_MAX(maxStrength, stng[vi])((maxStrength) > (stng[vi]) ? (maxStrength) : (stng[vi]));
689 }
690 }
691 }
692 flipProd = 1;
693 if (sctx->type!=seekTypeRidgeSurfaceOP &&
694 sctx->type!=seekTypeValleySurfaceOP) {
695 for (ei=0; ei<12; ei++) {
696 flip[ei] = sctx->flip[bag->evti[ei] + 5*si];
697 flipProd *= flip[ei];
698 }
699 }
700
701 if ((sctx->strengthUse && maxStrength < sctx->strength)
702 || !flipProd) {
703 /* either the corners this voxel don't meet strength,
704 or something else is funky */
705 for (vi=0; vi<8; vi++) {
706 vval[vi] = 0;
707 }
708 } else {
709 if (sctx->type==seekTypeRidgeSurfaceOP ||
710 sctx->type==seekTypeValleySurfaceOP) {
711 /* find orientation based on outer product rule */
712 double outer[9];
713 double outerevals[3],outerevecs[9];
714 ELL_3MV_OUTER(outer,evec[0],evec[0])((((outer)+0)[0] = ((evec[0])[0])*((evec[0]))[0], ((outer)+0)
[1] = ((evec[0])[0])*((evec[0]))[1], ((outer)+0)[2] = ((evec[
0])[0])*((evec[0]))[2]), (((outer)+3)[0] = ((evec[0])[1])*((evec
[0]))[0], ((outer)+3)[1] = ((evec[0])[1])*((evec[0]))[1], ((outer
)+3)[2] = ((evec[0])[1])*((evec[0]))[2]), (((outer)+6)[0] = (
(evec[0])[2])*((evec[0]))[0], ((outer)+6)[1] = ((evec[0])[2])
*((evec[0]))[1], ((outer)+6)[2] = ((evec[0])[2])*((evec[0]))[
2]));
715 for (vi=1; vi<8; ++vi) {
716 ELL_3MV_OUTER_INCR(outer,evec[vi],evec[vi])((((outer)+0)[0] += ((evec[vi])[0])*((evec[vi]))[0], ((outer)
+0)[1] += ((evec[vi])[0])*((evec[vi]))[1], ((outer)+0)[2] += (
(evec[vi])[0])*((evec[vi]))[2]), (((outer)+3)[0] += ((evec[vi
])[1])*((evec[vi]))[0], ((outer)+3)[1] += ((evec[vi])[1])*((evec
[vi]))[1], ((outer)+3)[2] += ((evec[vi])[1])*((evec[vi]))[2])
, (((outer)+6)[0] += ((evec[vi])[2])*((evec[vi]))[0], ((outer
)+6)[1] += ((evec[vi])[2])*((evec[vi]))[1], ((outer)+6)[2] +=
((evec[vi])[2])*((evec[vi]))[2]));
717 }
718 ell_3m_eigensolve_d(outerevals, outerevecs, outer, AIR_TRUE1);
719 for (vi=0; vi<8; ++vi) {
720 if (ELL_3V_DOT(evec[vi],outerevecs)((evec[vi])[0]*(outerevecs)[0] + (evec[vi])[1]*(outerevecs)[1
] + (evec[vi])[2]*(outerevecs)[2])<0)
721 ELL_3V_SCALE(evec[vi], -1.0, evec[vi])((evec[vi])[0] = (-1.0)*(evec[vi])[0], (evec[vi])[1] = (-1.0)
*(evec[vi])[1], (evec[vi])[2] = (-1.0)*(evec[vi])[2]);
722 }
723 } else {
724 ELL_3V_SCALE(evec[1], flip[0], evec[1])((evec[1])[0] = (flip[0])*(evec[1])[0], (evec[1])[1] = (flip[
0])*(evec[1])[1], (evec[1])[2] = (flip[0])*(evec[1])[2]);
725 ELL_3V_SCALE(evec[2], flip[1], evec[2])((evec[2])[0] = (flip[1])*(evec[2])[0], (evec[2])[1] = (flip[
1])*(evec[2])[1], (evec[2])[2] = (flip[1])*(evec[2])[2]);
726 ELL_3V_SCALE(evec[3], flip[0]*flip[2], evec[3])((evec[3])[0] = (flip[0]*flip[2])*(evec[3])[0], (evec[3])[1] =
(flip[0]*flip[2])*(evec[3])[1], (evec[3])[2] = (flip[0]*flip
[2])*(evec[3])[2]);
727 ELL_3V_SCALE(evec[4], flip[4], evec[4])((evec[4])[0] = (flip[4])*(evec[4])[0], (evec[4])[1] = (flip[
4])*(evec[4])[1], (evec[4])[2] = (flip[4])*(evec[4])[2]);
728 ELL_3V_SCALE(evec[5], flip[4]*flip[8], evec[5])((evec[5])[0] = (flip[4]*flip[8])*(evec[5])[0], (evec[5])[1] =
(flip[4]*flip[8])*(evec[5])[1], (evec[5])[2] = (flip[4]*flip
[8])*(evec[5])[2]);
729 ELL_3V_SCALE(evec[6], flip[4]*flip[9], evec[6])((evec[6])[0] = (flip[4]*flip[9])*(evec[6])[0], (evec[6])[1] =
(flip[4]*flip[9])*(evec[6])[1], (evec[6])[2] = (flip[4]*flip
[9])*(evec[6])[2]);
730 ELL_3V_SCALE(evec[7], flip[4]*flip[8]*flip[10], evec[7])((evec[7])[0] = (flip[4]*flip[8]*flip[10])*(evec[7])[0], (evec
[7])[1] = (flip[4]*flip[8]*flip[10])*(evec[7])[1], (evec[7])[
2] = (flip[4]*flip[8]*flip[10])*(evec[7])[2]);
731 }
732 for (vi=0; vi<8; vi++) {
733 vval[vi] = ELL_3V_DOT(grad[vi], evec[vi])((grad[vi])[0]*(evec[vi])[0] + (grad[vi])[1]*(evec[vi])[1] + (
grad[vi])[2]*(evec[vi])[2]);
734 }
735 }
736 return;
737}
738
739static int
740triangulate(seekContext *sctx, baggage *bag, limnPolyData *lpld) {
741 /* static const char me[]="triangulate"; */
742 unsigned xi, yi, sx, sy, si, spi;
743 /* ========================================================== */
744 /* NOTE: these things must agree with information in tables.c */
745 int e2v[12][2] = { /* maps edge index to corner vertex indices */
746 {0, 1}, /* 0 */
747 {0, 2}, /* 1 */
748 {1, 3}, /* 2 */
749 {2, 3}, /* 3 */
750 {0, 4}, /* 4 */
751 {1, 5}, /* 5 */
752 {2, 6}, /* 6 */
753 {3, 7}, /* 7 */
754 {4, 5}, /* 8 */
755 {4, 6}, /* 9 */
756 {5, 7}, /* 10 */
757 {6, 7} /* 11 */
758 };
759 double vccoord[8][3] = { /* vertex corner coordinates */
760 {0, 0, 0}, /* 0 */
761 {1, 0, 0}, /* 1 */
762 {0, 1, 0}, /* 2 */
763 {1, 1, 0}, /* 3 */
764 {0, 0, 1}, /* 4 */
765 {1, 0, 1}, /* 5 */
766 {0, 1, 1}, /* 6 */
767 {1, 1, 1} /* 7 */
768 };
769 /* ========================================================== */
770
771 sx = AIR_CAST(unsigned int, sctx->sx)((unsigned int)(sctx->sx));
772 sy = AIR_CAST(unsigned int, sctx->sy)((unsigned int)(sctx->sy));
773
774 for (yi=0; yi<sy-1; yi++) {
775 double vval[8], vgrad[8][3], vert[3], tvertA[4], tvertB[4], ww;
776 unsigned char vcase;
777 int ti, vi, ei, vi0, vi1, ecase;
778 const int *tcase;
779 unsigned int vii[3];
780 for (xi=0; xi<sx-1; xi++) {
781 si = xi + sx*yi;
782 spi = (xi+1) + (sx+2)*(yi+1);
783 switch (sctx->type) {
784 case seekTypeIsocontour:
785 vvalIsoSet(sctx, bag, vval, xi, yi);
786 break;
787 case seekTypeRidgeSurface:
788 case seekTypeValleySurface:
789 case seekTypeMaximalSurface:
790 case seekTypeMinimalSurface:
791 case seekTypeRidgeSurfaceOP:
792 case seekTypeValleySurfaceOP:
793 vvalSurfSet(sctx, bag, vval, xi, yi);
794 break;
795 }
796 /* determine voxel and edge case */
797 vcase = 0;
798 for (vi=0; vi<8; vi++) {
799 vcase |= (vval[vi] > 0) << vi;
800 }
801 if (0 == vcase || 255 == vcase) {
802 /* no triangles added here */
803 continue;
804 }
805 /* set voxel corner gradients */
806 if (seekTypeIsocontour == sctx->type
807 && sctx->normalsFind
808 && !sctx->normAns) {
809 voxelGrads(vgrad, sctx->sclv, sx, spi);
810 }
811 sctx->voxNum++;
812 ecase = seekContour3DTopoHackEdge[vcase];
813 /* create new vertices as needed */
814 for (ei=0; ei<12; ei++) {
815 if ((ecase & (1 << ei))
816 && -1 == sctx->vidx[bag->evti[ei] + 5*si]) {
817 int ovi;
818 double tvec[3], grad[3], tlen;
819 /* this edge is needed for triangulation,
820 and, we haven't already created a vertex for it */
821 vi0 = e2v[ei][0];
822 vi1 = e2v[ei][1];
823 ww = vval[vi0]/(vval[vi0] - vval[vi1]);
824 ELL_3V_LERP(vert, ww, vccoord[vi0], vccoord[vi1])((vert)[0] = (((ww))*(((vccoord[vi1])[0]) - ((vccoord[vi0])[0
])) + ((vccoord[vi0])[0])), (vert)[1] = (((ww))*(((vccoord[vi1
])[1]) - ((vccoord[vi0])[1])) + ((vccoord[vi0])[1])), (vert)[
2] = (((ww))*(((vccoord[vi1])[2]) - ((vccoord[vi0])[2])) + ((
vccoord[vi0])[2])));
825 ELL_4V_SET(tvertA, vert[0] + xi, vert[1] + yi, vert[2] + bag->zi, 1)((tvertA)[0] = (vert[0] + xi), (tvertA)[1] = (vert[1] + yi), (
tvertA)[2] = (vert[2] + bag->zi), (tvertA)[3] = (1));
826 ELL_4MV_MUL(tvertB, sctx->txfIdx, tvertA)((tvertB)[0]=(sctx->txfIdx)[ 0]*(tvertA)[0]+(sctx->txfIdx
)[ 1]*(tvertA)[1]+(sctx->txfIdx)[ 2]*(tvertA)[2]+(sctx->
txfIdx)[ 3]*(tvertA)[3], (tvertB)[1]=(sctx->txfIdx)[ 4]*(tvertA
)[0]+(sctx->txfIdx)[ 5]*(tvertA)[1]+(sctx->txfIdx)[ 6]*
(tvertA)[2]+(sctx->txfIdx)[ 7]*(tvertA)[3], (tvertB)[2]=(sctx
->txfIdx)[ 8]*(tvertA)[0]+(sctx->txfIdx)[ 9]*(tvertA)[1
]+(sctx->txfIdx)[10]*(tvertA)[2]+(sctx->txfIdx)[11]*(tvertA
)[3], (tvertB)[3]=(sctx->txfIdx)[12]*(tvertA)[0]+(sctx->
txfIdx)[13]*(tvertA)[1]+(sctx->txfIdx)[14]*(tvertA)[2]+(sctx
->txfIdx)[15]*(tvertA)[3]);
827 /* tvertB is now in input index space */
828 ELL_4MV_MUL(tvertA, sctx->shape->ItoW, tvertB)((tvertA)[0]=(sctx->shape->ItoW)[ 0]*(tvertB)[0]+(sctx->
shape->ItoW)[ 1]*(tvertB)[1]+(sctx->shape->ItoW)[ 2]
*(tvertB)[2]+(sctx->shape->ItoW)[ 3]*(tvertB)[3], (tvertA
)[1]=(sctx->shape->ItoW)[ 4]*(tvertB)[0]+(sctx->shape
->ItoW)[ 5]*(tvertB)[1]+(sctx->shape->ItoW)[ 6]*(tvertB
)[2]+(sctx->shape->ItoW)[ 7]*(tvertB)[3], (tvertA)[2]=(
sctx->shape->ItoW)[ 8]*(tvertB)[0]+(sctx->shape->
ItoW)[ 9]*(tvertB)[1]+(sctx->shape->ItoW)[10]*(tvertB)[
2]+(sctx->shape->ItoW)[11]*(tvertB)[3], (tvertA)[3]=(sctx
->shape->ItoW)[12]*(tvertB)[0]+(sctx->shape->ItoW
)[13]*(tvertB)[1]+(sctx->shape->ItoW)[14]*(tvertB)[2]+(
sctx->shape->ItoW)[15]*(tvertB)[3]);
829 /* tvertA is now in input world space */
830 ELL_4V_HOMOG(tvertA, tvertA)((tvertA)[0] = (tvertA)[0]/(tvertA)[3], (tvertA)[1] = (tvertA
)[1]/(tvertA)[3], (tvertA)[2] = (tvertA)[2]/(tvertA)[3], (tvertA
)[3] = 1.0);
831 ELL_4V_HOMOG(tvertB, tvertB)((tvertB)[0] = (tvertB)[0]/(tvertB)[3], (tvertB)[1] = (tvertB
)[1]/(tvertB)[3], (tvertB)[2] = (tvertB)[2]/(tvertB)[3], (tvertB
)[3] = 1.0);
832 ovi = sctx->vidx[bag->evti[ei] + 5*si] =
833 airArrayLenIncr(bag->xyzwArr, 1);
834 ELL_4V_SET_TT(lpld->xyzw + 4*ovi, float,((lpld->xyzw + 4*ovi)[0] = ((float)((tvertA[0]))), (lpld->
xyzw + 4*ovi)[1] = ((float)((tvertA[1]))), (lpld->xyzw + 4
*ovi)[2] = ((float)((tvertA[2]))), (lpld->xyzw + 4*ovi)[3]
= ((float)((1.0))))
835 tvertA[0], tvertA[1], tvertA[2], 1.0)((lpld->xyzw + 4*ovi)[0] = ((float)((tvertA[0]))), (lpld->
xyzw + 4*ovi)[1] = ((float)((tvertA[1]))), (lpld->xyzw + 4
*ovi)[2] = ((float)((tvertA[2]))), (lpld->xyzw + 4*ovi)[3]
= ((float)((1.0))));
836 /*
837 fprintf(stderr, "!%s: vert %u: %g %g %g\n", me, ovi,
838 tvertA[0], tvertA[1], tvertA[2]);
839 */
840 if (sctx->normalsFind) {
841 airArrayLenIncr(bag->normArr, 1);
842 if (sctx->normAns) {
843 gageProbe(sctx->gctx, tvertB[0], tvertB[1], tvertB[2]);
844 ELL_3V_SCALE_TT(lpld->norm + 3*ovi, float, -1, sctx->normAns)((lpld->norm + 3*ovi)[0] = ((float)((-1)*(sctx->normAns
)[0])), (lpld->norm + 3*ovi)[1] = ((float)((-1)*(sctx->
normAns)[1])), (lpld->norm + 3*ovi)[2] = ((float)((-1)*(sctx
->normAns)[2])));
845 if (sctx->reverse) {
846 ELL_3V_SCALE(lpld->norm + 3*ovi, -1, lpld->norm + 3*ovi)((lpld->norm + 3*ovi)[0] = (-1)*(lpld->norm + 3*ovi)[0]
, (lpld->norm + 3*ovi)[1] = (-1)*(lpld->norm + 3*ovi)[1
], (lpld->norm + 3*ovi)[2] = (-1)*(lpld->norm + 3*ovi)[
2]);
847 }
848 } else {
849 ELL_3V_LERP(grad, ww, vgrad[vi0], vgrad[vi1])((grad)[0] = (((ww))*(((vgrad[vi1])[0]) - ((vgrad[vi0])[0])) +
((vgrad[vi0])[0])), (grad)[1] = (((ww))*(((vgrad[vi1])[1]) -
((vgrad[vi0])[1])) + ((vgrad[vi0])[1])), (grad)[2] = (((ww))
*(((vgrad[vi1])[2]) - ((vgrad[vi0])[2])) + ((vgrad[vi0])[2]))
);
850 ELL_3MV_MUL(tvec, sctx->txfNormal, grad)((tvec)[0] = (sctx->txfNormal)[0]*(grad)[0] + (sctx->txfNormal
)[1]*(grad)[1] + (sctx->txfNormal)[2]*(grad)[2], (tvec)[1]
= (sctx->txfNormal)[3]*(grad)[0] + (sctx->txfNormal)[4
]*(grad)[1] + (sctx->txfNormal)[5]*(grad)[2], (tvec)[2] = (
sctx->txfNormal)[6]*(grad)[0] + (sctx->txfNormal)[7]*(grad
)[1] + (sctx->txfNormal)[8]*(grad)[2]);
851 ELL_3V_NORM_TT(lpld->norm + 3*ovi, float, tvec, tlen)(tlen = ((float)((sqrt((((tvec))[0]*((tvec))[0] + ((tvec))[1]
*((tvec))[1] + ((tvec))[2]*((tvec))[2]))))), ((lpld->norm +
3*ovi)[0] = ((float)((1.0/tlen)*(tvec)[0])), (lpld->norm +
3*ovi)[1] = ((float)((1.0/tlen)*(tvec)[1])), (lpld->norm +
3*ovi)[2] = ((float)((1.0/tlen)*(tvec)[2]))));
852 }
853 }
854 sctx->vertNum++;
855 /*
856 fprintf(stderr, "%s: vert %d (edge %d) of (%d,%d,%d) "
857 "at %g %g %g\n",
858 me, sctx->vidx[bag->evti[ei] + 5*si], ei, xi, yi, zi,
859 vert[0] + xi, vert[1] + yi, vert[2] + bag->zi);
860 */
861 }
862 }
863 /* add triangles */
864 ti = 0;
865 tcase = seekContour3DTopoHackTriangle[vcase];
866 while (-1 != tcase[0 + 3*ti]) {
867 unsigned iii;
868 ELL_3V_SET(vii,((vii)[0] = (sctx->vidx[bag->evti[tcase[0 + 3*ti]] + 5*
si]), (vii)[1] = (sctx->vidx[bag->evti[tcase[1 + 3*ti]]
+ 5*si]), (vii)[2] = (sctx->vidx[bag->evti[tcase[2 + 3
*ti]] + 5*si]))
869 sctx->vidx[bag->evti[tcase[0 + 3*ti]] + 5*si],((vii)[0] = (sctx->vidx[bag->evti[tcase[0 + 3*ti]] + 5*
si]), (vii)[1] = (sctx->vidx[bag->evti[tcase[1 + 3*ti]]
+ 5*si]), (vii)[2] = (sctx->vidx[bag->evti[tcase[2 + 3
*ti]] + 5*si]))
870 sctx->vidx[bag->evti[tcase[1 + 3*ti]] + 5*si],((vii)[0] = (sctx->vidx[bag->evti[tcase[0 + 3*ti]] + 5*
si]), (vii)[1] = (sctx->vidx[bag->evti[tcase[1 + 3*ti]]
+ 5*si]), (vii)[2] = (sctx->vidx[bag->evti[tcase[2 + 3
*ti]] + 5*si]))
871 sctx->vidx[bag->evti[tcase[2 + 3*ti]] + 5*si])((vii)[0] = (sctx->vidx[bag->evti[tcase[0 + 3*ti]] + 5*
si]), (vii)[1] = (sctx->vidx[bag->evti[tcase[1 + 3*ti]]
+ 5*si]), (vii)[2] = (sctx->vidx[bag->evti[tcase[2 + 3
*ti]] + 5*si]));
872 if (sctx->reverse) {
873 int tmpi;
874 tmpi = vii[1]; vii[1] = vii[2]; vii[2] = tmpi;
875 }
876 iii = airArrayLenIncr(bag->indxArr, 3);
877 ELL_3V_COPY(lpld->indx + iii, vii)((lpld->indx + iii)[0] = (vii)[0], (lpld->indx + iii)[1
] = (vii)[1], (lpld->indx + iii)[2] = (vii)[2]);
878 /*
879 fprintf(stderr, "!%s: tri %u: %u %u %u\n",
880 me, iii/3, vii[0], vii[1], vii[2]);
881 */
882 lpld->icnt[0] += 3;
883 sctx->faceNum++;
884 ti++;
885 }
886 }
887 }
888 return 0;
889}
890
891static int
892surfaceExtract(seekContext *sctx, limnPolyData *lpld) {
893 static const char me[]="surfaceExtract";
894 char done[AIR_STRLEN_SMALL(128+1)];
895 unsigned int zi, sz;
896 baggage *bag;
897
898 bag = baggageNew(sctx);
6
←
Calling 'baggageNew'
→
11
←
Returned allocated memory
→
899 sz = AIR_CAST(unsigned int, sctx->sz)((unsigned int)(sctx->sz));
900
901 /* this creates the airArrays in bag */
902 if (outputInit(sctx, bag, lpld)) {
12
←
Taking true branch
→
903 biffAddf(SEEKseekBiffKey, "%s: trouble", me);
13
←
Potential leak of memory pointed to by 'bag'
904 return 1;
905 }
906
907 if (sctx->verbose > 2) {
908 fprintf(stderr__stderrp, "%s: extracting ... ", me);
909 }
910 for (zi=0; zi<sz-1; zi++) {
911 char trouble=0;
912 if (sctx->verbose > 2) {
913 fprintf(stderr__stderrp, "%s", airDoneStr(0, zi, sz-2, done));
914 fflush(stderr__stderrp);
915 }
916 bag->zi = zi;
917 if (sctx->type==seekTypeRidgeSurfaceT ||
918 sctx->type==seekTypeValleySurfaceT) {
919 if (_seekShuffleProbeT(sctx, bag) ||
920 _seekTriangulateT(sctx, bag, lpld))
921 trouble = 1;
922 } else {
923 if (shuffleProbe(sctx, bag) ||
924 triangulate(sctx, bag, lpld))
925 trouble = 1;
926 }
927 if (trouble) {
928 biffAddf(SEEKseekBiffKey, "%s: trouble on zi = %u", me, zi);
929 return 1;
930 }
931 }
932 if (sctx->verbose > 2) {
933 fprintf(stderr__stderrp, "%s\n", airDoneStr(0, zi, sz-2, done));
934 }
935
936 /* this cleans up the airArrays in bag */
937 baggageNix(bag);
938
939 return 0;
940}
941
942int
943seekExtract(seekContext *sctx, limnPolyData *lpld) {
944 static const char me[]="seekExtract";
945 double time0;
946 int E;
947
948 if (!( sctx && lpld )) {
1
Taking false branch
→
949 biffAddf(SEEKseekBiffKey, "%s: got NULL pointer", me);
950 return 1;
951 }
952
953 if (seekTypeIsocontour == sctx->type) {
2
←
Taking false branch
→
954 if (!AIR_EXISTS(sctx->isovalue)(((int)(!((sctx->isovalue) - (sctx->isovalue)))))) {
955 biffAddf(SEEKseekBiffKey, "%s: didn't seem to ever set isovalue (now %g)", me,
956 sctx->isovalue);
957 return 1;
958 }
959 }
960
961 if (sctx->verbose) {
3
←
Taking false branch
→
962 fprintf(stderr__stderrp, "%s: --------------------\n", me);
963 fprintf(stderr__stderrp, "%s: flagResult = %d\n", me,
964 sctx->flag[flagResult]);
965 }
966
967 /* reset max strength seen */
968 sctx->strengthSeenMax = AIR_NAN(airFloatQNaN.f);
969
970 /* start time */
971 time0 = airTime();
972
973 switch(sctx->type) {
4
←
Control jumps to 'case seekTypeRidgeSurface:' at line 975
→
974 case seekTypeIsocontour:
975 case seekTypeRidgeSurface:
976 case seekTypeValleySurface:
977 case seekTypeMinimalSurface:
978 case seekTypeMaximalSurface:
979 case seekTypeRidgeSurfaceOP:
980 case seekTypeRidgeSurfaceT:
981 case seekTypeValleySurfaceOP:
982 case seekTypeValleySurfaceT:
983 E = surfaceExtract(sctx, lpld);
5
←
Calling 'surfaceExtract'
→
984 break;
985 default:
986 biffAddf(SEEKseekBiffKey, "%s: sorry, %s extraction not implemented", me,
987 airEnumStr(seekType, sctx->type));
988 return 1;
989 break;
990 }
991 if (E) {
992 biffAddf(SEEKseekBiffKey, "%s: trouble", me);
993 return 1;
994 }
995
996 /* end time */
997 sctx->time = airTime() - time0;
998
999 sctx->flag[flagResult] = AIR_FALSE0;
1000
1001 return 0;
1002}