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