1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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) 2011 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 limnPolyData *
limnPolyDataNew(void)29 limnPolyDataNew(void) {
30 limnPolyData *pld;
31
32 pld = AIR_CALLOC(1, limnPolyData);
33 if (pld) {
34 pld->xyzw = NULL;
35 pld->xyzwNum = 0;
36 pld->rgba = NULL;
37 pld->rgbaNum = 0;
38 pld->norm = NULL;
39 pld->normNum = 0;
40 pld->tex2 = NULL;
41 pld->tex2Num = 0;
42 pld->tang = NULL;
43 pld->tangNum = 0;
44 pld->indx = NULL;
45 pld->indxNum = 0;
46 pld->primNum = 0;
47 pld->type = NULL;
48 pld->icnt = NULL;
49 }
50 return pld;
51 }
52
53 limnPolyData *
limnPolyDataNix(limnPolyData * pld)54 limnPolyDataNix(limnPolyData *pld) {
55
56 if (pld) {
57 airFree(pld->xyzw);
58 airFree(pld->rgba);
59 airFree(pld->norm);
60 airFree(pld->tex2);
61 airFree(pld->tang);
62 airFree(pld->indx);
63 airFree(pld->type);
64 airFree(pld->icnt);
65 }
66 airFree(pld);
67 return NULL;
68 }
69
70 /*
71 ** doesn't set pld->xyzwNum, only the per-attribute xxxNum variables
72 */
73 int
_limnPolyDataInfoAlloc(limnPolyData * pld,unsigned int infoBitFlag,unsigned int vertNum)74 _limnPolyDataInfoAlloc(limnPolyData *pld, unsigned int infoBitFlag,
75 unsigned int vertNum) {
76 static const char me[]="_limnPolyDataInfoAlloc";
77
78 if (vertNum != pld->rgbaNum
79 && ((1 << limnPolyDataInfoRGBA) & infoBitFlag)) {
80 pld->rgba = (unsigned char *)airFree(pld->rgba);
81 if (vertNum) {
82 pld->rgba = AIR_CALLOC(vertNum*4, unsigned char);
83 if (!pld->rgba) {
84 biffAddf(LIMN, "%s: couldn't allocate %u rgba", me, vertNum);
85 return 1;
86 }
87 }
88 pld->rgbaNum = vertNum;
89 }
90
91 if (vertNum != pld->normNum
92 && ((1 << limnPolyDataInfoNorm) & infoBitFlag)) {
93 pld->norm = (float *)airFree(pld->norm);
94 if (vertNum) {
95 pld->norm = AIR_CALLOC(vertNum*4, float);
96 if (!pld->norm) {
97 biffAddf(LIMN, "%s: couldn't allocate %u norm", me, vertNum);
98 return 1;
99 }
100 }
101 pld->normNum = vertNum;
102 }
103
104 if (vertNum != pld->tex2Num
105 && ((1 << limnPolyDataInfoTex2) & infoBitFlag)) {
106 pld->tex2 = (float *)airFree(pld->tex2);
107 if (vertNum) {
108 pld->tex2 = AIR_CALLOC(vertNum*2, float);
109 if (!pld->tex2) {
110 biffAddf(LIMN, "%s: couldn't allocate %u tex2", me, vertNum);
111 return 1;
112 }
113 }
114 pld->tex2Num = vertNum;
115 }
116
117 if (vertNum != pld->tangNum
118 && ((1 << limnPolyDataInfoTang) & infoBitFlag)) {
119 pld->tang = (float *)airFree(pld->tang);
120 if (vertNum) {
121 pld->tang = AIR_CALLOC(vertNum*3, float);
122 if (!pld->tang) {
123 biffAddf(LIMN, "%s: couldn't allocate %u tang", me, vertNum);
124 return 1;
125 }
126 }
127 pld->tangNum = vertNum;
128 }
129
130 return 0;
131 }
132
133 unsigned int
limnPolyDataInfoBitFlag(const limnPolyData * pld)134 limnPolyDataInfoBitFlag(const limnPolyData *pld) {
135 unsigned int ret;
136
137 ret = 0;
138 if (pld) {
139 if (pld->rgba && pld->rgbaNum == pld->xyzwNum) {
140 ret |= (1 << limnPolyDataInfoRGBA);
141 }
142 if (pld->norm && pld->normNum == pld->xyzwNum) {
143 ret |= (1 << limnPolyDataInfoNorm);
144 }
145 if (pld->tex2 && pld->tex2Num == pld->xyzwNum) {
146 ret |= (1 << limnPolyDataInfoTex2);
147 }
148 if (pld->tang && pld->tangNum == pld->xyzwNum) {
149 ret |= (1 << limnPolyDataInfoTang);
150 }
151 }
152 return ret;
153 }
154
155 int
limnPolyDataAlloc(limnPolyData * pld,unsigned int infoBitFlag,unsigned int vertNum,unsigned int indxNum,unsigned int primNum)156 limnPolyDataAlloc(limnPolyData *pld,
157 unsigned int infoBitFlag,
158 unsigned int vertNum,
159 unsigned int indxNum,
160 unsigned int primNum) {
161 static const char me[]="limnPolyDataAlloc";
162
163 if (!pld) {
164 biffAddf(LIMN, "%s: got NULL pointer", me);
165 return 1;
166 }
167 if (vertNum != pld->xyzwNum) {
168 pld->xyzw = (float *)airFree(pld->xyzw);
169 if (vertNum) {
170 pld->xyzw = AIR_CALLOC(vertNum*4, float);
171 if (!pld->xyzw) {
172 biffAddf(LIMN, "%s: couldn't allocate %u xyzw", me, vertNum);
173 return 1;
174 }
175 }
176 pld->xyzwNum = vertNum;
177 }
178 if (_limnPolyDataInfoAlloc(pld, infoBitFlag, vertNum)) {
179 biffAddf(LIMN, "%s: couldn't allocate info", me);
180 return 1;
181 }
182 if (indxNum != pld->indxNum) {
183 pld->indx = (unsigned int *)airFree(pld->indx);
184 if (indxNum) {
185 pld->indx = AIR_CALLOC(indxNum, unsigned int);
186 if (!pld->indx) {
187 biffAddf(LIMN, "%s: couldn't allocate %u indices", me, indxNum);
188 return 1;
189 }
190 }
191 pld->indxNum = indxNum;
192 }
193 if (primNum != pld->primNum) {
194 pld->type = (unsigned char *)airFree(pld->type);
195 pld->icnt = (unsigned int *)airFree(pld->icnt);
196 if (primNum) {
197 pld->type = AIR_CALLOC(primNum, unsigned char);
198 pld->icnt = AIR_CALLOC(primNum, unsigned int);
199 if (!(pld->type && pld->icnt)) {
200 biffAddf(LIMN, "%s: couldn't allocate %u primitives", me, primNum);
201 return 1;
202 }
203 }
204 pld->primNum = primNum;
205 }
206 return 0;
207 }
208
209 size_t
limnPolyDataSize(const limnPolyData * pld)210 limnPolyDataSize(const limnPolyData *pld) {
211 size_t ret = 0;
212
213 if (pld) {
214 ret += pld->xyzwNum*sizeof(float)*4;
215 if (pld->rgba) {
216 ret += pld->rgbaNum*sizeof(unsigned char)*4;
217 }
218 if (pld->norm) {
219 ret += pld->normNum*sizeof(float)*3;
220 }
221 if (pld->tex2) {
222 ret += pld->tex2Num*sizeof(float)*2;
223 }
224 if (pld->tang) {
225 ret += pld->tangNum*sizeof(float)*3;
226 }
227 ret += pld->indxNum*sizeof(unsigned int);
228 ret += pld->primNum*sizeof(signed char);
229 ret += pld->primNum*sizeof(unsigned int);
230 }
231 return ret;
232 }
233
234 int
limnPolyDataCopy(limnPolyData * pldB,const limnPolyData * pldA)235 limnPolyDataCopy(limnPolyData *pldB, const limnPolyData *pldA) {
236 static const char me[]="limnPolyDataCopy";
237
238 if (!( pldB && pldA )) {
239 biffAddf(LIMN, "%s: got NULL pointer", me);
240 return 1;
241 }
242 if (limnPolyDataAlloc(pldB, limnPolyDataInfoBitFlag(pldA),
243 pldA->xyzwNum, pldA->indxNum, pldA->primNum)) {
244 biffAddf(LIMN, "%s: couldn't allocate output", me);
245 return 1;
246 }
247 memcpy(pldB->xyzw, pldA->xyzw, pldA->xyzwNum*sizeof(float)*4);
248 if (pldA->rgba) {
249 memcpy(pldB->rgba, pldA->rgba, pldA->rgbaNum*sizeof(unsigned char)*4);
250 }
251 if (pldA->norm) {
252 memcpy(pldB->norm, pldA->norm, pldA->normNum*sizeof(float)*3);
253 }
254 if (pldA->tex2) {
255 memcpy(pldB->tex2, pldA->tex2, pldA->tex2Num*sizeof(float)*2);
256 }
257 if (pldA->tang) {
258 memcpy(pldB->tang, pldA->tang, pldA->tangNum*sizeof(float)*3);
259 }
260 memcpy(pldB->indx, pldA->indx, pldA->indxNum*sizeof(unsigned int));
261 memcpy(pldB->type, pldA->type, pldA->primNum*sizeof(signed char));
262 memcpy(pldB->icnt, pldA->icnt, pldA->primNum*sizeof(unsigned int));
263 return 0;
264 }
265
266 int
limnPolyDataCopyN(limnPolyData * pldB,const limnPolyData * pldA,unsigned int num)267 limnPolyDataCopyN(limnPolyData *pldB, const limnPolyData *pldA,
268 unsigned int num) {
269 static const char me[]="limnPolyDataCopyN";
270 unsigned int ii, jj, size;
271
272 if (!( pldB && pldA )) {
273 biffAddf(LIMN, "%s: got NULL pointer", me);
274 return 1;
275 }
276 if (limnPolyDataAlloc(pldB, limnPolyDataInfoBitFlag(pldA),
277 num*pldA->xyzwNum,
278 num*pldA->indxNum,
279 num*pldA->primNum)) {
280 biffAddf(LIMN, "%s: couldn't allocate output", me);
281 return 1;
282 }
283 for (ii=0; ii<num; ii++) {
284 /* fprintf(stderr, "!%s: ii = %u/%u\n", me, ii, num); */
285 size = pldA->xyzwNum*4;
286 /*
287 char *_beg = (char *)(pldB->xyzw + ii*size);
288 char *_end = _beg + size - 1;
289 fprintf(stderr, "!%s: memcpy(%p+%u=%p,%u) --> [%p,%p] inside: %d %d\n", me,
290 pldB->xyzw, ii*size, pldB->xyzw + ii*size, size,
291 _beg, _end, AIR_IN_CL(_xyzwBeg, _beg, _xyzwEnd),
292 AIR_IN_CL(_xyzwBeg, _end, _xyzwEnd));
293 */
294 memcpy(pldB->xyzw + ii*size, pldA->xyzw, size*sizeof(float));
295 for (jj=0; jj<pldA->indxNum; jj++) {
296 (pldB->indx + ii*pldA->indxNum)[jj] = pldA->indx[jj] + ii*pldA->xyzwNum;
297 }
298 size = pldA->primNum;
299 memcpy(pldB->type + ii*size, pldA->type, size*sizeof(unsigned char));
300 memcpy(pldB->icnt + ii*size, pldA->icnt, size*sizeof(unsigned int));
301 if (pldA->rgba) {
302 size = pldA->rgbaNum*4;
303 memcpy(pldB->rgba + ii*size, pldA->rgba, size*sizeof(unsigned char));
304 }
305 if (pldA->norm) {
306 size = pldA->normNum*3;
307 memcpy(pldB->norm + ii*size, pldA->norm, size*sizeof(float));
308 }
309 if (pldA->tex2) {
310 size = pldA->tex2Num*2;
311 memcpy(pldB->tex2 + ii*size, pldA->tex2, size*sizeof(float));
312 }
313 if (pldA->tang) {
314 size = pldA->tangNum*3;
315 memcpy(pldB->tang + ii*size, pldA->tang, size*sizeof(float));
316 }
317 }
318 return 0;
319 }
320
321 /*
322 ******** limnPolyDataTransform_f, limnPolyDataTransform_d
323 **
324 ** transforms a surface (both positions, and normals (if set))
325 ** by given homogenous transform
326 */
327 void
limnPolyDataTransform_f(limnPolyData * pld,const float homat[16])328 limnPolyDataTransform_f(limnPolyData *pld,
329 const float homat[16]) {
330 double hovec[4], mat[9], inv[9], norm[3], tang[3], nmat[9], tlen;
331 unsigned int vertIdx;
332
333 if (pld && homat) {
334 ELL_34M_EXTRACT(mat, homat);
335 if (pld->norm) {
336 ell_3m_inv_d(inv, mat);
337 ELL_3M_TRANSPOSE(nmat, inv);
338 } else {
339 ELL_3M_IDENTITY_SET(nmat); /* hush unused value compiler warnings */
340 }
341 for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
342 ELL_4MV_MUL(hovec, homat, pld->xyzw + 4*vertIdx);
343 ELL_4V_COPY_TT(pld->xyzw + 4*vertIdx, float, hovec);
344 if (pld->norm) {
345 ELL_3MV_MUL(norm, nmat, pld->norm + 3*vertIdx);
346 ELL_3V_NORM(norm, norm, tlen);
347 ELL_3V_COPY_TT(pld->norm + 3*vertIdx, float, norm);
348 }
349 if (pld->tang) { /* HEY: not tested */
350 ELL_3MV_MUL(tang, mat, pld->tang + 3*vertIdx);
351 ELL_3V_NORM(tang, tang, tlen);
352 ELL_3V_COPY_TT(pld->tang + 3*vertIdx, float, tang);
353 }
354 }
355 }
356 return;
357 }
358
359 /* !!! COPY AND PASTE !!! (except for double homat[16]) */
360 void
limnPolyDataTransform_d(limnPolyData * pld,const double homat[16])361 limnPolyDataTransform_d(limnPolyData *pld, const double homat[16]) {
362 double hovec[4], mat[9], inv[9], norm[3], nmat[9];
363 unsigned int vertIdx;
364
365 if (pld && homat) {
366 if (pld->norm) {
367 ELL_34M_EXTRACT(mat, homat);
368 ell_3m_inv_d(inv, mat);
369 ELL_3M_TRANSPOSE(nmat, inv);
370 } else {
371 ELL_3M_IDENTITY_SET(nmat); /* hush unused value compiler warnings */
372 }
373 for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
374 ELL_4MV_MUL(hovec, homat, pld->xyzw + 4*vertIdx);
375 ELL_4V_COPY_TT(pld->xyzw + 4*vertIdx, float, hovec);
376 if (pld->norm) {
377 ELL_3MV_MUL(norm, nmat, pld->norm + 3*vertIdx);
378 ELL_3V_COPY_TT(pld->norm + 3*vertIdx, float, norm);
379 }
380 }
381 }
382 return;
383 }
384
385 unsigned int
limnPolyDataPolygonNumber(const limnPolyData * pld)386 limnPolyDataPolygonNumber(const limnPolyData *pld) {
387 unsigned int ret, primIdx;
388
389 ret = 0;
390 if (pld) {
391 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
392 switch(pld->type[primIdx]) {
393 case limnPrimitiveNoop:
394 /* no increment */
395 break;
396 case limnPrimitiveTriangles:
397 ret += pld->icnt[primIdx]/3;
398 break;
399 case limnPrimitiveTriangleStrip:
400 case limnPrimitiveTriangleFan:
401 ret += pld->icnt[primIdx] - 2;
402 break;
403 case limnPrimitiveQuads:
404 ret += pld->icnt[primIdx]/4;
405 break;
406 }
407 }
408 }
409 return ret;
410 }
411
412 static int
limnPolyDataVertexNormals_(limnPolyData * pld,int nonorient)413 limnPolyDataVertexNormals_(limnPolyData *pld, int nonorient) {
414 static const char me[]="limnPolyDataVertexNormals_";
415 unsigned int infoBitFlag, primIdx, triIdx, normIdx, baseVertIdx;
416 double len, *matrix=NULL;
417 airArray *mop;
418
419 infoBitFlag = limnPolyDataInfoBitFlag(pld);
420 if (limnPolyDataAlloc(pld,
421 infoBitFlag | (1 << limnPolyDataInfoNorm),
422 pld->xyzwNum,
423 pld->indxNum,
424 pld->primNum)) {
425 biffAddf(LIMN, "%s: couldn't alloc polydata normals", me);
426 return 1;
427 }
428
429 mop = airMopNew();
430 if (nonorient) { /* we will accumulate outer products */
431 matrix = (double *) malloc(sizeof(double)*9*pld->xyzwNum);
432 if (matrix==NULL) {
433 biffAddf(LIMN, "%s: couldn't allocate matrix buffer", me);
434 return 1;
435 }
436 airMopAdd(mop, matrix, airFree, airMopAlways);
437 for (normIdx=0; normIdx<pld->normNum; normIdx++) {
438 ELL_3M_ZERO_SET(matrix + 9*normIdx);
439 }
440 } else {
441 for (normIdx=0; normIdx<pld->normNum; normIdx++) {
442 ELL_3V_SET(pld->norm + 3*normIdx, 0, 0, 0);
443 }
444 }
445
446 baseVertIdx = 0;
447 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
448 unsigned int triNum, indxLine[3]={0,0,0}, ii;
449 float pos[3][3], edgeA[3], edgeB[3], norm[3];
450
451 switch (pld->type[primIdx]) {
452 case limnPrimitiveTriangles:
453 triNum = pld->icnt[primIdx]/3;
454 break;
455 case limnPrimitiveTriangleStrip:
456 case limnPrimitiveTriangleFan:
457 triNum = pld->icnt[primIdx]-2;
458 break;
459 case limnPrimitiveNoop:
460 triNum = 0;
461 break;
462 default:
463 biffAddf(LIMN, "%s: came across unsupported limnPrimitive \"%s\"", me,
464 airEnumStr(limnPrimitive, pld->type[primIdx]));
465 airMopError(mop);
466 return 1;
467 }
468
469 if (limnPrimitiveNoop != pld->type[primIdx]) {
470 for (triIdx=0; triIdx<triNum; triIdx++) {
471 switch (pld->type[primIdx]) {
472 case limnPrimitiveTriangles:
473 ELL_3V_COPY(indxLine, pld->indx + baseVertIdx + 3*triIdx);
474 break;
475 case limnPrimitiveTriangleStrip:
476 if (triIdx%2==0) {
477 ELL_3V_COPY(indxLine, pld->indx+baseVertIdx+triIdx);
478 } else {
479 ELL_3V_SET(indxLine, pld->indx[baseVertIdx+triIdx+1],
480 pld->indx[baseVertIdx+triIdx],
481 pld->indx[baseVertIdx+triIdx+2]);
482 }
483 break;
484 case limnPrimitiveTriangleFan:
485 ELL_3V_SET(indxLine, pld->indx[baseVertIdx],
486 pld->indx[baseVertIdx+triIdx+1],
487 pld->indx[baseVertIdx+triIdx+2]);
488 break;
489 }
490 for (ii=0; ii<3; ii++) {
491 ELL_34V_HOMOG(pos[ii], pld->xyzw + 4*indxLine[ii]);
492 }
493 ELL_3V_SUB(edgeA, pos[1], pos[0]);
494 ELL_3V_SUB(edgeB, pos[2], pos[0]);
495 ELL_3V_CROSS(norm, edgeA, edgeB);
496 /* Adding cross products without any normalization is
497 * equivalent to weighting by triangle area, as proposed
498 * (among others) by G. Taubin ("Estimating the tensor of
499 * curvature of a surface from a polyhedral approximation",
500 * ICCV 1995). This is efficient, avoids trouble with
501 * degenerate triangles and gives reasonable results in
502 * practice. */
503 if (nonorient) {
504 double outer[9];
505 ELL_3MV_OUTER(outer, norm, norm);
506 len = ELL_3V_LEN(norm);
507 /* re-scale so that we don't weight by square of area */
508 for (ii=0; ii<3; ii++) {
509 ELL_3M_SCALE_INCR(matrix+9*indxLine[ii], 1.0/len, outer);
510 }
511 } else {
512 for (ii=0; ii<3; ii++) {
513 ELL_3V_INCR(pld->norm + 3*indxLine[ii], norm);
514 }
515 }
516 }
517 }
518 baseVertIdx += pld->icnt[primIdx];
519 }
520
521 if (nonorient) {
522 double eval[3], evec[9];
523 for (normIdx=0; normIdx<pld->normNum; normIdx++) {
524 ell_3m_eigensolve_d(eval, evec, matrix+9*normIdx, AIR_TRUE);
525 ELL_3V_COPY_TT(pld->norm + 3*normIdx, float, evec);
526 }
527 } else {
528 for (normIdx=0; normIdx<pld->normNum; normIdx++) {
529 ELL_3V_NORM_TT(pld->norm + 3*normIdx, float, pld->norm + 3*normIdx, len);
530 }
531 }
532
533 airMopOkay(mop);
534 return 0;
535 }
536
537 int
limnPolyDataVertexNormals(limnPolyData * pld)538 limnPolyDataVertexNormals(limnPolyData *pld) {
539 return limnPolyDataVertexNormals_(pld, 0);
540 }
541
542 /* counterpart for non-orientable surfaces */
543 int
limnPolyDataVertexNormalsNO(limnPolyData * pld)544 limnPolyDataVertexNormalsNO(limnPolyData *pld) {
545 return limnPolyDataVertexNormals_(pld, 1);
546 }
547
548 unsigned int
limnPolyDataPrimitiveTypes(const limnPolyData * pld)549 limnPolyDataPrimitiveTypes(const limnPolyData *pld) {
550 unsigned int ret, primIdx;
551
552 ret = 0;
553 if (pld) {
554 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
555 ret |= (1 << pld->type[primIdx]);
556 }
557 }
558 return ret;
559 }
560
561 int
limnPolyDataPrimitiveVertexNumber(Nrrd * nout,limnPolyData * pld)562 limnPolyDataPrimitiveVertexNumber(Nrrd *nout, limnPolyData *pld) {
563 static const char me[]="limnPolyDataPrimitiveVertexNumber";
564 unsigned int *vnum, pidx;
565
566 if (!(nout && pld)) {
567 biffAddf(LIMN, "%s: got NULL pointer", me);
568 return 1;
569 }
570 if (nrrdMaybeAlloc_va(nout, nrrdTypeUInt, 1,
571 AIR_CAST(size_t, pld->primNum))) {
572 biffMovef(LIMN, NRRD, "%s: couldn't allocate output", me);
573 return 1;
574 }
575
576 vnum = AIR_CAST(unsigned int *, nout->data);
577 for (pidx=0; pidx<pld->primNum; pidx++) {
578 vnum[pidx] = pld->icnt[pidx];
579 }
580
581 return 0;
582 }
583
584 int
limnPolyDataPrimitiveArea(Nrrd * nout,limnPolyData * pld)585 limnPolyDataPrimitiveArea(Nrrd *nout, limnPolyData *pld) {
586 static const char me[]="limnPolyDataPrimitiveArea";
587 unsigned int primIdx, baseVertIdx;
588 unsigned int triNum, triIdx, *indx, ii;
589 float vert[3][3], edgeA[3], edgeB[3], cross[3];
590 double *area;
591
592 if (!(nout && pld)) {
593 biffAddf(LIMN, "%s: got NULL pointer", me);
594 return 1;
595 }
596 if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 1,
597 AIR_CAST(size_t, pld->primNum))) {
598 biffMovef(LIMN, NRRD, "%s: couldn't allocate output", me);
599 return 1;
600 }
601
602 area = AIR_CAST(double *, nout->data);
603 baseVertIdx = 0;
604 for (primIdx=0; primIdx<pld->primNum; primIdx++) {
605 area[primIdx] = 0;
606 switch (pld->type[primIdx]) {
607 case limnPrimitiveNoop:
608 area[primIdx] = 0.0;
609 break;
610 case limnPrimitiveTriangles:
611 triNum = pld->icnt[primIdx]/3;
612 for (triIdx=0; triIdx<triNum; triIdx++) {
613 indx = pld->indx + baseVertIdx + 3*triIdx;
614 for (ii=0; ii<3; ii++) {
615 ELL_34V_HOMOG(vert[ii], pld->xyzw + 4*indx[ii]);
616 }
617 ELL_3V_SUB(edgeA, vert[1], vert[0]);
618 ELL_3V_SUB(edgeB, vert[2], vert[0]);
619 ELL_3V_CROSS(cross, edgeA, edgeB);
620 area[primIdx] += ELL_3V_LEN(cross)/2;
621 }
622 break;
623 case limnPrimitiveTriangleStrip:
624 case limnPrimitiveTriangleFan:
625 case limnPrimitiveQuads:
626 biffAddf(LIMN,
627 "%s: sorry, haven't implemented area(prim[%u]=%s) yet", me,
628 primIdx, airEnumStr(limnPrimitive, pld->type[primIdx]));
629 return 1;
630 break;
631 case limnPrimitiveLineStrip:
632 /* lines have no area */
633 break;
634 }
635 baseVertIdx += pld->icnt[primIdx];
636 }
637
638 return 0;
639 }
640
641 /*
642 ** I may regret making this only be axis-aligned ...
643 */
644 int
limnPolyDataRasterize(Nrrd * nout,limnPolyData * pld,double min[3],double max[3],size_t size[3],int type)645 limnPolyDataRasterize(Nrrd *nout, limnPolyData *pld,
646 double min[3], double max[3],
647 size_t size[3], int type) {
648 static const char me[]="limnPolyDataRasterize";
649 size_t xi, yi, zi;
650 unsigned int vertIdx;
651 double (*ins)(void *, size_t, double);
652
653 if (!(nout && pld && min && max && size)) {
654 biffAddf(LIMN, "%s: got NULL pointer", me);
655 return 1;
656 }
657 if (airEnumValCheck(nrrdType, type)) {
658 biffAddf(LIMN, "%s: got invalid %s %d", me, nrrdType->name, type);
659 return 1;
660 }
661 if (nrrdTypeBlock == type) {
662 biffAddf(LIMN, "%s: can't use output type %s", me,
663 airEnumStr(nrrdType, type));
664 return 1;
665 }
666 if (!( min[0] < max[0] &&
667 min[1] < max[1] &&
668 min[2] < max[2] )) {
669 biffAddf(LIMN, "%s min (%g,%g,%g) not < max (%g,%g,%g)", me,
670 min[0], min[1], min[2], max[0], max[1], max[2]);
671 return 1;
672 }
673
674 if (nrrdMaybeAlloc_nva(nout, type, 3, size)) {
675 biffMovef(LIMN, NRRD, "%s: trouble allocating output", me);
676 return 1;
677 }
678 ins = nrrdDInsert[type];
679
680 for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
681 double xyz[3];
682
683 ELL_34V_HOMOG(xyz, pld->xyzw + 4*vertIdx);
684 if (!( AIR_IN_OP(min[0], xyz[0], max[0]) &&
685 AIR_IN_OP(min[1], xyz[1], max[1]) &&
686 AIR_IN_OP(min[2], xyz[2], max[2]) )) {
687 continue;
688 }
689 xi = airIndex(min[0], xyz[0], max[0], size[0]);
690 yi = airIndex(min[1], xyz[1], max[1], size[1]);
691 zi = airIndex(min[2], xyz[2], max[2], size[2]);
692 ins(nout->data, xi + size[0]*(yi + size[1]*zi), 1.0);
693 }
694
695 nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMin, min);
696 nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMax, max);
697
698 return 0;
699 }
700
701 void
limnPolyDataColorSet(limnPolyData * pld,unsigned char RR,unsigned char GG,unsigned char BB,unsigned char AA)702 limnPolyDataColorSet(limnPolyData *pld,
703 unsigned char RR, unsigned char GG,
704 unsigned char BB, unsigned char AA) {
705 unsigned int vertIdx;
706
707 if (pld && ((1 << limnPolyDataInfoRGBA) & limnPolyDataInfoBitFlag(pld))) {
708 for (vertIdx=0; vertIdx<pld->rgbaNum; vertIdx++) {
709 ELL_4V_SET(pld->rgba + 4*vertIdx, RR, GG, BB, AA);
710 }
711 }
712 return;
713 }
714