Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/limn/polydata.c	Lines:	0	361	0.0 %
Date:	2017-05-26	Branches:	0	250	0.0 %


/*
  Teem: Tools to process and visualize scientific data and images             .
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah
  Copyright (C) 2011  Thomas Schultz

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/


#include "limn.h"

limnPolyData *
limnPolyDataNew(void) {
  limnPolyData *pld;

  pld = AIR_CALLOC(1, limnPolyData);
  if (pld) {
    pld->xyzw = NULL;
    pld->xyzwNum = 0;
    pld->rgba = NULL;
    pld->rgbaNum = 0;
    pld->norm = NULL;
    pld->normNum = 0;
    pld->tex2 = NULL;
    pld->tex2Num = 0;
    pld->tang = NULL;
    pld->tangNum = 0;
    pld->indx = NULL;
    pld->indxNum = 0;
    pld->primNum = 0;
    pld->type = NULL;
    pld->icnt = NULL;
  }
  return pld;
}

limnPolyData *
limnPolyDataNix(limnPolyData *pld) {

  if (pld) {
    airFree(pld->xyzw);
    airFree(pld->rgba);
    airFree(pld->norm);
    airFree(pld->tex2);
    airFree(pld->tang);
    airFree(pld->indx);
    airFree(pld->type);
    airFree(pld->icnt);
  }
  airFree(pld);
  return NULL;
}

/*
** doesn't set pld->xyzwNum, only the per-attribute xxxNum variables
*/
int
_limnPolyDataInfoAlloc(limnPolyData *pld, unsigned int infoBitFlag,
                       unsigned int vertNum) {
  static const char me[]="_limnPolyDataInfoAlloc";

  if (vertNum != pld->rgbaNum
      && ((1 << limnPolyDataInfoRGBA) & infoBitFlag)) {
    pld->rgba = (unsigned char *)airFree(pld->rgba);
    if (vertNum) {
      pld->rgba = AIR_CALLOC(vertNum*4, unsigned char);
      if (!pld->rgba) {
        biffAddf(LIMN, "%s: couldn't allocate %u rgba", me, vertNum);
        return 1;
      }
    }
    pld->rgbaNum = vertNum;
  }

  if (vertNum != pld->normNum
      && ((1 << limnPolyDataInfoNorm) & infoBitFlag)) {
    pld->norm = (float *)airFree(pld->norm);
    if (vertNum) {
      pld->norm = AIR_CALLOC(vertNum*4, float);
      if (!pld->norm) {
        biffAddf(LIMN, "%s: couldn't allocate %u norm", me, vertNum);
        return 1;
      }
    }
    pld->normNum = vertNum;
  }

  if (vertNum != pld->tex2Num
      && ((1 << limnPolyDataInfoTex2) & infoBitFlag)) {
    pld->tex2 = (float *)airFree(pld->tex2);
    if (vertNum) {
      pld->tex2 = AIR_CALLOC(vertNum*2, float);
      if (!pld->tex2) {
        biffAddf(LIMN, "%s: couldn't allocate %u tex2", me, vertNum);
        return 1;
      }
    }
    pld->tex2Num = vertNum;
  }

  if (vertNum != pld->tangNum
      && ((1 << limnPolyDataInfoTang) & infoBitFlag)) {
    pld->tang = (float *)airFree(pld->tang);
    if (vertNum) {
      pld->tang = AIR_CALLOC(vertNum*3, float);
      if (!pld->tang) {
        biffAddf(LIMN, "%s: couldn't allocate %u tang", me, vertNum);
        return 1;
      }
    }
    pld->tangNum = vertNum;
  }

  return 0;
}

unsigned int
limnPolyDataInfoBitFlag(const limnPolyData *pld) {
  unsigned int ret;

  ret = 0;
  if (pld) {
    if (pld->rgba && pld->rgbaNum == pld->xyzwNum) {
      ret |= (1 << limnPolyDataInfoRGBA);
    }
    if (pld->norm && pld->normNum == pld->xyzwNum) {
      ret |= (1 << limnPolyDataInfoNorm);
    }
    if (pld->tex2 && pld->tex2Num == pld->xyzwNum) {
      ret |= (1 << limnPolyDataInfoTex2);
    }
    if (pld->tang && pld->tangNum == pld->xyzwNum) {
      ret |= (1 << limnPolyDataInfoTang);
    }
  }
  return ret;
}

int
limnPolyDataAlloc(limnPolyData *pld,
                  unsigned int infoBitFlag,
                  unsigned int vertNum,
                  unsigned int indxNum,
                  unsigned int primNum) {
  static const char me[]="limnPolyDataAlloc";

  if (!pld) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (vertNum != pld->xyzwNum) {
    pld->xyzw = (float *)airFree(pld->xyzw);
    if (vertNum) {
      pld->xyzw = AIR_CALLOC(vertNum*4, float);
      if (!pld->xyzw) {
        biffAddf(LIMN, "%s: couldn't allocate %u xyzw", me, vertNum);
        return 1;
      }
    }
    pld->xyzwNum = vertNum;
  }
  if (_limnPolyDataInfoAlloc(pld, infoBitFlag, vertNum)) {
    biffAddf(LIMN, "%s: couldn't allocate info", me);
    return 1;
  }
  if (indxNum != pld->indxNum) {
    pld->indx = (unsigned int *)airFree(pld->indx);
    if (indxNum) {
      pld->indx = AIR_CALLOC(indxNum, unsigned int);
      if (!pld->indx) {
        biffAddf(LIMN, "%s: couldn't allocate %u indices", me, indxNum);
        return 1;
      }
    }
    pld->indxNum = indxNum;
  }
  if (primNum != pld->primNum) {
    pld->type = (unsigned char *)airFree(pld->type);
    pld->icnt = (unsigned int *)airFree(pld->icnt);
    if (primNum) {
      pld->type = AIR_CALLOC(primNum, unsigned char);
      pld->icnt = AIR_CALLOC(primNum, unsigned int);
      if (!(pld->type && pld->icnt)) {
        biffAddf(LIMN, "%s: couldn't allocate %u primitives", me, primNum);
        return 1;
      }
    }
    pld->primNum = primNum;
  }
  return 0;
}

size_t
limnPolyDataSize(const limnPolyData *pld) {
  size_t ret = 0;

  if (pld) {
    ret += pld->xyzwNum*sizeof(float)*4;
    if (pld->rgba) {
      ret += pld->rgbaNum*sizeof(unsigned char)*4;
    }
    if (pld->norm) {
      ret += pld->normNum*sizeof(float)*3;
    }
    if (pld->tex2) {
      ret += pld->tex2Num*sizeof(float)*2;
    }
    if (pld->tang) {
      ret += pld->tangNum*sizeof(float)*3;
    }
    ret += pld->indxNum*sizeof(unsigned int);
    ret += pld->primNum*sizeof(signed char);
    ret += pld->primNum*sizeof(unsigned int);
  }
  return ret;
}

int
limnPolyDataCopy(limnPolyData *pldB, const limnPolyData *pldA) {
  static const char me[]="limnPolyDataCopy";

  if (!( pldB && pldA )) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (limnPolyDataAlloc(pldB, limnPolyDataInfoBitFlag(pldA),
                        pldA->xyzwNum, pldA->indxNum, pldA->primNum)) {
    biffAddf(LIMN, "%s: couldn't allocate output", me);
    return 1;
  }
  memcpy(pldB->xyzw, pldA->xyzw, pldA->xyzwNum*sizeof(float)*4);
  if (pldA->rgba) {
    memcpy(pldB->rgba, pldA->rgba, pldA->rgbaNum*sizeof(unsigned char)*4);
  }
  if (pldA->norm) {
    memcpy(pldB->norm, pldA->norm, pldA->normNum*sizeof(float)*3);
  }
  if (pldA->tex2) {
    memcpy(pldB->tex2, pldA->tex2, pldA->tex2Num*sizeof(float)*2);
  }
  if (pldA->tang) {
    memcpy(pldB->tang, pldA->tang, pldA->tangNum*sizeof(float)*3);
  }
  memcpy(pldB->indx, pldA->indx, pldA->indxNum*sizeof(unsigned int));
  memcpy(pldB->type, pldA->type, pldA->primNum*sizeof(signed char));
  memcpy(pldB->icnt, pldA->icnt, pldA->primNum*sizeof(unsigned int));
  return 0;
}

int
limnPolyDataCopyN(limnPolyData *pldB, const limnPolyData *pldA,
                  unsigned int num) {
  static const char me[]="limnPolyDataCopyN";
  unsigned int ii, jj, size;

  if (!( pldB && pldA )) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (limnPolyDataAlloc(pldB, limnPolyDataInfoBitFlag(pldA),
                        num*pldA->xyzwNum,
                        num*pldA->indxNum,
                        num*pldA->primNum)) {
    biffAddf(LIMN, "%s: couldn't allocate output", me);
    return 1;
  }
  for (ii=0; ii<num; ii++) {
    /* fprintf(stderr, "!%s: ii = %u/%u\n", me, ii, num); */
    size = pldA->xyzwNum*4;
    /*
    char *_beg = (char *)(pldB->xyzw + ii*size);
    char *_end = _beg + size - 1;
    fprintf(stderr, "!%s: memcpy(%p+%u=%p,%u) --> [%p,%p] inside: %d %d\n", me,
            pldB->xyzw, ii*size, pldB->xyzw + ii*size, size,
            _beg, _end, AIR_IN_CL(_xyzwBeg, _beg, _xyzwEnd),
            AIR_IN_CL(_xyzwBeg, _end, _xyzwEnd));
    */
    memcpy(pldB->xyzw + ii*size, pldA->xyzw, size*sizeof(float));
    for (jj=0; jj<pldA->indxNum; jj++) {
      (pldB->indx + ii*pldA->indxNum)[jj] = pldA->indx[jj] + ii*pldA->xyzwNum;
    }
    size = pldA->primNum;
    memcpy(pldB->type + ii*size, pldA->type, size*sizeof(unsigned char));
    memcpy(pldB->icnt + ii*size, pldA->icnt, size*sizeof(unsigned int));
    if (pldA->rgba) {
      size = pldA->rgbaNum*4;
      memcpy(pldB->rgba + ii*size, pldA->rgba, size*sizeof(unsigned char));
    }
    if (pldA->norm) {
      size = pldA->normNum*3;
      memcpy(pldB->norm + ii*size, pldA->norm, size*sizeof(float));
    }
    if (pldA->tex2) {
      size = pldA->tex2Num*2;
      memcpy(pldB->tex2 + ii*size, pldA->tex2, size*sizeof(float));
    }
    if (pldA->tang) {
      size = pldA->tangNum*3;
      memcpy(pldB->tang + ii*size, pldA->tang, size*sizeof(float));
    }
  }
  return 0;
}

/*
******** limnPolyDataTransform_f, limnPolyDataTransform_d
**
** transforms a surface (both positions, and normals (if set))
** by given homogenous transform
*/
void
limnPolyDataTransform_f(limnPolyData *pld,
                        const float homat[16]) {
  double hovec[4], mat[9], inv[9], norm[3], tang[3], nmat[9], tlen;
  unsigned int vertIdx;

  if (pld && homat) {
    ELL_34M_EXTRACT(mat, homat);
    if (pld->norm) {
      ell_3m_inv_d(inv, mat);
      ELL_3M_TRANSPOSE(nmat, inv);
    } else {
      ELL_3M_IDENTITY_SET(nmat);  /* hush unused value compiler warnings */
    }
    for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
      ELL_4MV_MUL(hovec, homat, pld->xyzw + 4*vertIdx);
      ELL_4V_COPY_TT(pld->xyzw + 4*vertIdx, float, hovec);
      if (pld->norm) {
        ELL_3MV_MUL(norm, nmat, pld->norm + 3*vertIdx);
        ELL_3V_NORM(norm, norm, tlen);
        ELL_3V_COPY_TT(pld->norm + 3*vertIdx, float, norm);
      }
      if (pld->tang) { /* HEY: not tested */
        ELL_3MV_MUL(tang, mat, pld->tang + 3*vertIdx);
        ELL_3V_NORM(tang, tang, tlen);
        ELL_3V_COPY_TT(pld->tang + 3*vertIdx, float, tang);
      }
    }
  }
  return;
}

/* !!! COPY AND PASTE !!!  (except for double homat[16]) */
void
limnPolyDataTransform_d(limnPolyData *pld, const double homat[16]) {
  double hovec[4], mat[9], inv[9], norm[3], nmat[9];
  unsigned int vertIdx;

  if (pld && homat) {
    if (pld->norm) {
      ELL_34M_EXTRACT(mat, homat);
      ell_3m_inv_d(inv, mat);
      ELL_3M_TRANSPOSE(nmat, inv);
    } else {
      ELL_3M_IDENTITY_SET(nmat);  /* hush unused value compiler warnings */
    }
    for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
      ELL_4MV_MUL(hovec, homat, pld->xyzw + 4*vertIdx);
      ELL_4V_COPY_TT(pld->xyzw + 4*vertIdx, float, hovec);
      if (pld->norm) {
        ELL_3MV_MUL(norm, nmat, pld->norm + 3*vertIdx);
        ELL_3V_COPY_TT(pld->norm + 3*vertIdx, float, norm);
      }
    }
  }
  return;
}

unsigned int
limnPolyDataPolygonNumber(const limnPolyData *pld) {
  unsigned int ret, primIdx;

  ret = 0;
  if (pld) {
    for (primIdx=0; primIdx<pld->primNum; primIdx++) {
      switch(pld->type[primIdx]) {
      case limnPrimitiveNoop:
        /* no increment */
        break;
      case limnPrimitiveTriangles:
        ret += pld->icnt[primIdx]/3;
        break;
      case limnPrimitiveTriangleStrip:
      case limnPrimitiveTriangleFan:
        ret += pld->icnt[primIdx] - 2;
        break;
      case limnPrimitiveQuads:
        ret += pld->icnt[primIdx]/4;
        break;
      }
    }
  }
  return ret;
}

static int
limnPolyDataVertexNormals_(limnPolyData *pld, int nonorient) {
  static const char me[]="limnPolyDataVertexNormals_";
  unsigned int infoBitFlag, primIdx, triIdx, normIdx, baseVertIdx;
  double len, *matrix=NULL;
  airArray *mop;

  infoBitFlag = limnPolyDataInfoBitFlag(pld);
  if (limnPolyDataAlloc(pld,
                        infoBitFlag | (1 << limnPolyDataInfoNorm),
                        pld->xyzwNum,
                        pld->indxNum,
                        pld->primNum)) {
    biffAddf(LIMN, "%s: couldn't alloc polydata normals", me);
    return 1;
  }

  mop = airMopNew();
  if (nonorient) { /* we will accumulate outer products */
    matrix = (double *) malloc(sizeof(double)*9*pld->xyzwNum);
    if (matrix==NULL) {
      biffAddf(LIMN, "%s: couldn't allocate matrix buffer", me);
      return 1;
    }
    airMopAdd(mop, matrix, airFree, airMopAlways);
    for (normIdx=0; normIdx<pld->normNum; normIdx++) {
      ELL_3M_ZERO_SET(matrix + 9*normIdx);
    }
  } else {
    for (normIdx=0; normIdx<pld->normNum; normIdx++) {
      ELL_3V_SET(pld->norm + 3*normIdx, 0, 0, 0);
    }
  }

  baseVertIdx = 0;
  for (primIdx=0; primIdx<pld->primNum; primIdx++) {
    unsigned int triNum, indxLine[3]={0,0,0}, ii;
    float pos[3][3], edgeA[3], edgeB[3], norm[3];

    switch (pld->type[primIdx]) {
    case limnPrimitiveTriangles:
      triNum = pld->icnt[primIdx]/3;
      break;
    case limnPrimitiveTriangleStrip:
    case limnPrimitiveTriangleFan:
      triNum = pld->icnt[primIdx]-2;
      break;
    case limnPrimitiveNoop:
      triNum = 0;
      break;
    default:
      biffAddf(LIMN, "%s: came across unsupported limnPrimitive \"%s\"", me,
               airEnumStr(limnPrimitive, pld->type[primIdx]));
      airMopError(mop);
      return 1;
    }

    if (limnPrimitiveNoop != pld->type[primIdx]) {
      for (triIdx=0; triIdx<triNum; triIdx++) {
        switch (pld->type[primIdx]) {
        case limnPrimitiveTriangles:
          ELL_3V_COPY(indxLine, pld->indx + baseVertIdx + 3*triIdx);
          break;
        case limnPrimitiveTriangleStrip:
          if (triIdx%2==0) {
            ELL_3V_COPY(indxLine, pld->indx+baseVertIdx+triIdx);
          } else {
            ELL_3V_SET(indxLine, pld->indx[baseVertIdx+triIdx+1],
                       pld->indx[baseVertIdx+triIdx],
                       pld->indx[baseVertIdx+triIdx+2]);
          }
          break;
        case limnPrimitiveTriangleFan:
          ELL_3V_SET(indxLine, pld->indx[baseVertIdx],
                     pld->indx[baseVertIdx+triIdx+1],
                     pld->indx[baseVertIdx+triIdx+2]);
          break;
        }
        for (ii=0; ii<3; ii++) {
          ELL_34V_HOMOG(pos[ii], pld->xyzw + 4*indxLine[ii]);
        }
        ELL_3V_SUB(edgeA, pos[1], pos[0]);
        ELL_3V_SUB(edgeB, pos[2], pos[0]);
        ELL_3V_CROSS(norm, edgeA, edgeB);
        /* Adding cross products without any normalization is
         * equivalent to weighting by triangle area, as proposed
         * (among others) by G. Taubin ("Estimating the tensor of
         * curvature of a surface from a polyhedral approximation",
         * ICCV 1995). This is efficient, avoids trouble with
         * degenerate triangles and gives reasonable results in
         * practice. */
        if (nonorient) {
          double outer[9];
          ELL_3MV_OUTER(outer, norm, norm);
          len = ELL_3V_LEN(norm);
          /* re-scale so that we don't weight by square of area */
          for (ii=0; ii<3; ii++) {
            ELL_3M_SCALE_INCR(matrix+9*indxLine[ii], 1.0/len, outer);
          }
        } else {
          for (ii=0; ii<3; ii++) {
            ELL_3V_INCR(pld->norm + 3*indxLine[ii], norm);
          }
        }
      }
    }
    baseVertIdx += pld->icnt[primIdx];
  }

  if (nonorient) {
    double eval[3], evec[9];
    for (normIdx=0; normIdx<pld->normNum; normIdx++) {
      ell_3m_eigensolve_d(eval, evec, matrix+9*normIdx, AIR_TRUE);
      ELL_3V_COPY_TT(pld->norm + 3*normIdx, float, evec);
    }
  } else {
    for (normIdx=0; normIdx<pld->normNum; normIdx++) {
      ELL_3V_NORM_TT(pld->norm + 3*normIdx, float, pld->norm + 3*normIdx, len);
    }
  }

  airMopOkay(mop);
  return 0;
}

int
limnPolyDataVertexNormals(limnPolyData *pld) {
  return limnPolyDataVertexNormals_(pld, 0);
}

/* counterpart for non-orientable surfaces */
int
limnPolyDataVertexNormalsNO(limnPolyData *pld) {
  return limnPolyDataVertexNormals_(pld, 1);
}

unsigned int
limnPolyDataPrimitiveTypes(const limnPolyData *pld) {
  unsigned int ret, primIdx;

  ret = 0;
  if (pld) {
    for (primIdx=0; primIdx<pld->primNum; primIdx++) {
      ret |= (1 << pld->type[primIdx]);
    }
  }
  return ret;
}

int
limnPolyDataPrimitiveVertexNumber(Nrrd *nout, limnPolyData *pld) {
  static const char me[]="limnPolyDataPrimitiveVertexNumber";
  unsigned int *vnum, pidx;

  if (!(nout && pld)) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (nrrdMaybeAlloc_va(nout, nrrdTypeUInt, 1,
                        AIR_CAST(size_t, pld->primNum))) {
    biffMovef(LIMN, NRRD, "%s: couldn't allocate output", me);
    return 1;
  }

  vnum = AIR_CAST(unsigned int *, nout->data);
  for (pidx=0; pidx<pld->primNum; pidx++) {
    vnum[pidx] = pld->icnt[pidx];
  }

  return 0;
}

int
limnPolyDataPrimitiveArea(Nrrd *nout, limnPolyData *pld) {
  static const char me[]="limnPolyDataPrimitiveArea";
  unsigned int primIdx, baseVertIdx;
  unsigned int triNum, triIdx, *indx, ii;
  float vert[3][3], edgeA[3], edgeB[3], cross[3];
  double *area;

  if (!(nout && pld)) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 1,
                        AIR_CAST(size_t, pld->primNum))) {
    biffMovef(LIMN, NRRD, "%s: couldn't allocate output", me);
    return 1;
  }

  area = AIR_CAST(double *, nout->data);
  baseVertIdx = 0;
  for (primIdx=0; primIdx<pld->primNum; primIdx++) {
    area[primIdx] = 0;
    switch (pld->type[primIdx]) {
    case limnPrimitiveNoop:
      area[primIdx] = 0.0;
      break;
    case limnPrimitiveTriangles:
      triNum = pld->icnt[primIdx]/3;
      for (triIdx=0; triIdx<triNum; triIdx++) {
        indx = pld->indx + baseVertIdx + 3*triIdx;
        for (ii=0; ii<3; ii++) {
          ELL_34V_HOMOG(vert[ii], pld->xyzw + 4*indx[ii]);
        }
        ELL_3V_SUB(edgeA, vert[1], vert[0]);
        ELL_3V_SUB(edgeB, vert[2], vert[0]);
        ELL_3V_CROSS(cross, edgeA, edgeB);
        area[primIdx] += ELL_3V_LEN(cross)/2;
      }
      break;
    case limnPrimitiveTriangleStrip:
    case limnPrimitiveTriangleFan:
    case limnPrimitiveQuads:
      biffAddf(LIMN,
               "%s: sorry, haven't implemented area(prim[%u]=%s) yet", me,
               primIdx, airEnumStr(limnPrimitive, pld->type[primIdx]));
      return 1;
      break;
    case limnPrimitiveLineStrip:
      /* lines have no area */
      break;
    }
    baseVertIdx += pld->icnt[primIdx];
  }

  return 0;
}

/*
** I may regret making this only be axis-aligned ...
*/
int
limnPolyDataRasterize(Nrrd *nout, limnPolyData *pld,
                      double min[3], double max[3],
                      size_t size[3], int type) {
  static const char me[]="limnPolyDataRasterize";
  size_t xi, yi, zi;
  unsigned int vertIdx;
  double (*ins)(void *, size_t, double);

  if (!(nout && pld && min && max && size)) {
    biffAddf(LIMN, "%s: got NULL pointer", me);
    return 1;
  }
  if (airEnumValCheck(nrrdType, type)) {
    biffAddf(LIMN, "%s: got invalid %s %d", me, nrrdType->name, type);
    return 1;
  }
  if (nrrdTypeBlock == type) {
    biffAddf(LIMN, "%s: can't use output type %s", me,
             airEnumStr(nrrdType, type));
    return 1;
  }
  if (!( min[0] < max[0] &&
         min[1] < max[1] &&
         min[2] < max[2] )) {
    biffAddf(LIMN, "%s min (%g,%g,%g) not < max (%g,%g,%g)", me,
             min[0], min[1], min[2], max[0], max[1], max[2]);
    return 1;
  }

  if (nrrdMaybeAlloc_nva(nout, type, 3, size)) {
    biffMovef(LIMN, NRRD, "%s: trouble allocating output", me);
    return 1;
  }
  ins = nrrdDInsert[type];

  for (vertIdx=0; vertIdx<pld->xyzwNum; vertIdx++) {
    double xyz[3];

    ELL_34V_HOMOG(xyz, pld->xyzw + 4*vertIdx);
    if (!( AIR_IN_OP(min[0], xyz[0], max[0]) &&
           AIR_IN_OP(min[1], xyz[1], max[1]) &&
           AIR_IN_OP(min[2], xyz[2], max[2]) )) {
      continue;
    }
    xi = airIndex(min[0], xyz[0], max[0], size[0]);
    yi = airIndex(min[1], xyz[1], max[1], size[1]);
    zi = airIndex(min[2], xyz[2], max[2], size[2]);
    ins(nout->data, xi + size[0]*(yi + size[1]*zi), 1.0);
  }

  nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMin, min);
  nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMax, max);

  return 0;
}

void
limnPolyDataColorSet(limnPolyData *pld,
                     unsigned char RR, unsigned char GG,
                     unsigned char BB, unsigned char AA) {
  unsigned int vertIdx;

  if (pld && ((1 << limnPolyDataInfoRGBA) & limnPolyDataInfoBitFlag(pld))) {
    for (vertIdx=0; vertIdx<pld->rgbaNum; vertIdx++) {
      ELL_4V_SET(pld->rgba + 4*vertIdx, RR, GG, BB, AA);
    }
  }
  return;
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/*
2			Teem: Tools to process and visualize scientific data and images .
3			Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4			Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5			Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6			Copyright (C) 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 *
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 *
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
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
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
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
210			limnPolyDataSize(const limnPolyData *pld) {
211			size_t ret = 0;
212
213			if (pld) {
214			ret += pld->xyzwNumsizeof(float)4;
215			if (pld->rgba) {
216			ret += pld->rgbaNumsizeof(unsigned char)4;
217			}
218			if (pld->norm) {
219			ret += pld->normNumsizeof(float)3;
220			}
221			if (pld->tex2) {
222			ret += pld->tex2Numsizeof(float)2;
223			}
224			if (pld->tang) {
225			ret += pld->tangNumsizeof(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
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->xyzwNumsizeof(float)4);
248			if (pldA->rgba) {
249			memcpy(pldB->rgba, pldA->rgba, pldA->rgbaNumsizeof(unsigned char)4);
250			}
251			if (pldA->norm) {
252			memcpy(pldB->norm, pldA->norm, pldA->normNumsizeof(float)3);
253			}
254			if (pldA->tex2) {
255			memcpy(pldB->tex2, pldA->tex2, pldA->tex2Numsizeof(float)2);
256			}
257			if (pldA->tang) {
258			memcpy(pldB->tang, pldA->tang, pldA->tangNumsizeof(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
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, iisize, pldB->xyzw + iisize, size,
291			_beg, _end, AIR_IN_CL(_xyzwBeg, _beg, _xyzwEnd),
292			AIR_IN_CL(_xyzwBeg, _end, _xyzwEnd));
293			*/
294			memcpy(pldB->xyzw + iisize, pldA->xyzw, sizesizeof(float));
295			for (jj=0; jj<pldA->indxNum; jj++) {
296			(pldB->indx + iipldA->indxNum)[jj] = pldA->indx[jj] + iipldA->xyzwNum;
297			}
298			size = pldA->primNum;
299			memcpy(pldB->type + iisize, pldA->type, sizesizeof(unsigned char));
300			memcpy(pldB->icnt + iisize, pldA->icnt, sizesizeof(unsigned int));
301			if (pldA->rgba) {
302			size = pldA->rgbaNum*4;
303			memcpy(pldB->rgba + iisize, pldA->rgba, sizesizeof(unsigned char));
304			}
305			if (pldA->norm) {
306			size = pldA->normNum*3;
307			memcpy(pldB->norm + iisize, pldA->norm, sizesizeof(float));
308			}
309			if (pldA->tex2) {
310			size = pldA->tex2Num*2;
311			memcpy(pldB->tex2 + iisize, pldA->tex2, sizesizeof(float));
312			}
313			if (pldA->tang) {
314			size = pldA->tangNum*3;
315			memcpy(pldB->tang + iisize, pldA->tang, sizesizeof(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
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
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
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
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 + 3normIdx, float, pld->norm + 3normIdx, len);
530			}
531			}
532
533			airMopOkay(mop);
534			return 0;
535			}
536
537			int
538			limnPolyDataVertexNormals(limnPolyData *pld) {
539			return limnPolyDataVertexNormals_(pld, 0);
540			}
541
542			/* counterpart for non-orientable surfaces */
543			int
544			limnPolyDataVertexNormalsNO(limnPolyData *pld) {
545			return limnPolyDataVertexNormals_(pld, 1);
546			}
547
548			unsigned int
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
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
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
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
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			}