Nurbs.cpp

///--------------------------------------------------------------------------------------
/// Non-Uniform Rational B-Splines, NuPatches
///--------------------------------------------------------------------------------------
///
/// Nurbs evaluation funcs are an adaptation of 'nurb_polyg' code from GraphicsGemsIV.
///
/// Currently avoids actual splitting and just adjusts clipping parameters.
/// Bound calculation is probably also not quite correct, currently done by evaluating
/// the control points of the patch.

#include "Nurbs.h"
#include "Bound.h"
#include "Transform.h"
#include "Mathutil.h"
#include "Framework.h"
#include "Hider.h"
#include "MicroPolygonGrid.h"
#include "Attributes.h"

__BEGIN_QDRENDER

//------------------------------------------------------------------------------

/*
 * Return the current knot the parameter u is less than or equal to.
 * Find this "breakpoint" allows the evaluation routines to concentrate on
 * only those control points actually effecting the curve around u.]
 *
 *    m   is the number of points on the curve (or surface direction)
 *    k   is the order of the curve (or surface direction)
 *    kv  is the knot vector ([0..m+k-1]) to find the break point in.
 */
inline int FindBreakPoint(float u, const fsArray_t<float>& kv, int m, int k)
{
      if (u == kv[m + 1]) return m; // Special case for closed interval
      /*
      int i = m + k;
      while ((u < kv[i]) && (i > 0)) i--;
      return i;
      */
      // do bin.search instead of above linear, might be faster (for larger arrays at least)
      int lo = 0, hi = m+k;
      while (lo != hi) {
            const int md = (hi + lo) >> 1;
            if (u >= kv[md]) lo = md + 1; else hi = md;
      }
      return ((hi == 0) ? 0 : (hi - 1));
      
}

/*
 * Compute Bi,k(u), for i = 0..k.
 *  u       is the parameter of the spline to find the basis functions for
 *  brkPoint      is the start of the knot interval ("segment")
 *  kv            is the knot vector
 *  k       is the order of the curve
 *  bvals   is the array of returned basis values.
 *
 * (From Bartels, Beatty & Barsky, p.387)
 */
inline void BasisFunctions(float u, int brkPoint, const fsArray_t<float>& kv, int k, fsArray_t<float>& bvals)
{
      bvals[0] = 1.f;
      for (int r=2; r<=k; ++r) {
            bvals[r - 1] = 0.f;
            for (int s=(r-2), i=(brkPoint-r+2); s>=0; --s, ++i) {
                  const float dv = kv[i + r - 1] - kv[i];
                  const float omega = (i < 0) ? 0.f : ((u - kv[i]) / ((dv==0.f) ? 1.f : dv));
                  bvals[s + 1] += (1 - omega)*bvals[s];
                  bvals[s] *= omega;
            }
      }
}

// Compute derivatives of the basis functions Bi,k(u)'
inline void BasisDerivatives(float u, int brkPoint, const fsArray_t<float>& kv, int k, fsArray_t<float>& dvals)
{
      BasisFunctions(u, brkPoint, kv, k-1, dvals);
      dvals[k-1] = 0.f;     // BasisFunctions misses this
      const float knotScale = kv[brkPoint + 1] - kv[brkPoint];
      for (int s=(k-2), i=(brkPoint-k+2); s>=0; --s, ++i) {
            const float dv = kv[i + k - 1] - kv[i];
            const float omega = knotScale*(k - 1)/((dv==0.f) ? 1.f : dv);
            dvals[s + 1] += -omega*dvals[s];
            dvals[s] *= omega;
      }
}

/*
 * Calculate a point p on NurbSurface n at a specific u, v using the tensor product.
 * If utan and vtan are not nil, compute the u and v tangents as well.
 *
 * Note the valid parameter range for u and v is
 * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV])
 */
void NuPatch::eval(float u, float v, RtPoint P, RtVector dPdu, RtVector dPdv)
{
      // clamp u & v to make sure it stays within range (precision errors cause havoc)
      u = CLAMP(u, umin, umax);
      v = CLAMP(v, vmin, vmax);
      // Evaluate non-uniform basis functions (and derivatives)
      const int ubrkPoint = FindBreakPoint(u, nudata->uknot, nudata->nu-1, nudata->uorder);
      const int ufirst = ubrkPoint - nudata->uorder + 1;
      BasisFunctions(u, ubrkPoint, nudata->uknot, nudata->uorder, nudata->bu);
      if (dPdu) BasisDerivatives(u, ubrkPoint, nudata->uknot, nudata->uorder, nudata->buprime);

      const int vbrkPoint = FindBreakPoint(v, nudata->vknot, nudata->nv-1, nudata->vorder);
      const int vfirst = vbrkPoint - nudata->vorder + 1;
      BasisFunctions(v, vbrkPoint, nudata->vknot, nudata->vorder, nudata->bv);
      if (dPdv) BasisDerivatives(v, vbrkPoint, nudata->vknot, nudata->vorder, nudata->bvprime);

      // Weight control points against the basis functions
      RtHpoint r = {0,0,0,0}, rutan = {0,0,0,0}, rvtan = {0,0,0,0};
      for (int i=0; i<nudata->vorder; ++i) {
            RtHpoint* pw = &nudata->Pw[(i+vfirst)*nudata->nu + ufirst];
            for (int j=0; j<nudata->uorder; ++j) {
                  const int ri = nudata->vorder - 1 - i, rj = nudata->uorder - 1 - j;
                  float tmp = nudata->bu[rj] * nudata->bv[ri];
                  r[0] += pw[j][0] * tmp;
                  r[1] += pw[j][1] * tmp;
                  r[2] += pw[j][2] * tmp;
                  r[3] += pw[j][3] * tmp;
                  if (dPdu) {
                        tmp = nudata->buprime[rj] * nudata->bv[ri];
                        rutan[0] += pw[j][0] * tmp;
                        rutan[1] += pw[j][1] * tmp;
                        rutan[2] += pw[j][2] * tmp;
                        rutan[3] += pw[j][3] * tmp;
                  }
                  if (dPdv) {
                        tmp = nudata->bu[rj] * nudata->bvprime[ri];
                        rvtan[0] += pw[j][0] * tmp;
                        rvtan[1] += pw[j][1] * tmp;
                        rvtan[2] += pw[j][2] * tmp;
                        rvtan[3] += pw[j][3] * tmp;
                  }
            }
      }

      // Project tangents, using the quotient rule for differentiation
      const float dv = (r[3]==0.f) ? 1.f : (1.f/r[3]);
      const float wsqrdiv = dv*dv;
      if (dPdu && dPdv) {
            dPdu[0] = (r[3]*rutan[0] - rutan[3]*r[0])*wsqrdiv;
            dPdu[1] = (r[3]*rutan[1] - rutan[3]*r[1])*wsqrdiv;
            dPdu[2] = (r[3]*rutan[2] - rutan[3]*r[2])*wsqrdiv;
            dPdv[0] = (r[3]*rvtan[0] - rvtan[3]*r[0])*wsqrdiv;
            dPdv[1] = (r[3]*rvtan[1] - rvtan[3]*r[1])*wsqrdiv;
            dPdv[2] = (r[3]*rvtan[2] - rvtan[3]*r[2])*wsqrdiv;
      }
      P[0] = r[0]*dv, P[1] = r[1]*dv, P[2] = r[2]*dv;
}

//--------------------------------------------------------------------------
// NuPatch

NuPatch::NuPatch(RtInt _nu, RtInt _uorder, RtFloat _uknot[], RtFloat _umin, RtFloat _umax,
                 RtInt _nv, RtInt _vorder, RtFloat _vknot[], RtFloat _vmin, RtFloat _vmax,
                 RtInt n, RtToken tokens[], RtPointer parms[])
                 : nudata(NULL), umin(_umin), vmin(_vmin), umax(_umax), vmax(_vmax)
{
      // all checking for correctness done in ri/ribInterface
      nudata = new NuData(_nu, _uorder, _uknot, _nv, _vorder, _vknot);
      const int nump = _nu * _nv;
      RtHpoint* Pw = nudata->Pw = new RtHpoint[nump];
      // get control points, P/Pw, in case of 'P', w = 1
      for (int i=0; i<n; ++i) {
            if (!strcmp(tokens[i], RI_P)) {
                  const RtPoint* rP = reinterpret_cast<RtPoint*>(parms[i]);
                  for (int j=0; j<nump; ++j)
                        Pw[j][0] = rP[j][0], Pw[j][1] = rP[j][1], Pw[j][2] = rP[j][2], Pw[j][3] = 1.f;
                  break;
            }
            else if (!strcmp(tokens[i], RI_PW)) {
                  memcpy(Pw, parms[i], sizeof(RtHpoint)*nump);
                  break;
            }
      }
      // dont' forget reference count in duplicateSelf (must do this differently similar to Parametric)
      const int nuseg = _nu - _uorder + 1, nvseg = _nv - _vorder + 1;
      Primitive::initPrimVars(n, tokens, parms, nuseg*nvseg, (nuseg+1)*(nvseg+1), nump, nuseg*nvseg*4);
}

NuPatch::~NuPatch()
{
      if (nudata && (--nudata->refc == 0)) {
            delete nudata;
            nudata = NULL;
      }
}

void NuPatch::post_init()
{
      // only done after initial prim insert in state
      RtHpoint* Pw = nudata->Pw;
      const RtMatrix* const ri_mtx = xform->getRtMatrixPtr();
      const int nump = nudata->nu * nudata->nv;
      for (int i=0; i<nump; ++i)
            mulPMP4(Pw[i], *ri_mtx, Pw[i]);
}

Bound NuPatch::bound()
{
      Bound b;
      // evaluate at cp.verts for bound.
      // NOTE while it seems to work sofar, it can't be right for all cases,
      // surface might have a 'bump' or 'dent' inbetween control points...
      const float du = (umax - umin)/(float)(nudata->nu - 1), dv = (vmax - vmin)/(float)(nudata->nv - 1);
      int ug, vg;
      float u, v;
      RtPoint P;
      for (vg=0, v=vmin; vg<nudata->nv; ++vg, v+=dv)
            for (ug=0, u=umin; ug<nudata->nu; ++ug, u+=du) {
                  eval(u, v, P);
                  b.include(Point3(P[0], P[1], P[2]));
            }
      b.addEpsilon();
      return b;
}

// splitting is done here simply by only adjusting the clipping coords...
void NuPatch::split(const Framework &f, bool usplit, bool vsplit, splitbprims_t* spb)
{
      //printf("SPLITTING : %d %d\n", usplit, vsplit);
      if (usplit && vsplit) {
            // Left to bottom & top 1
            NuPatch* LT = static_cast<NuPatch*>(duplicateSelf());
            NuPatch *LB = static_cast<NuPatch*>(duplicateSelf());
            const float umid = (umin + umax)*0.5f;
            const float vmid = (vmin + vmax)*0.5f;
            LB->umax = LT->umax = umid;
            LB->vmax = LT->vmin = vmid;
            f.insert(LT);
            f.insert(LB);
            // Right to bottom & top 2
            NuPatch *RT = static_cast<NuPatch*>(duplicateSelf());
            NuPatch *RB = static_cast<NuPatch*>(duplicateSelf());
            RB->umin = RT->umin = umid;
            RB->vmax = RT->vmin = vmid;
            f.insert(RT);
            f.insert(RB);
      }
      else {
            if (usplit) {
                  NuPatch *L = static_cast<NuPatch*>(duplicateSelf());
                  NuPatch *R = static_cast<NuPatch*>(duplicateSelf());
                  const float umid = (umin + umax)*0.5f;
                  L->umax = umid, R->umin = umid;
                  f.insert(L);
                  f.insert(R);
            }
            else if (vsplit) {
                  NuPatch *T = static_cast<NuPatch*>(duplicateSelf());
                  NuPatch *B = static_cast<NuPatch*>(duplicateSelf());
                  const float vmid = (vmin + vmax)*0.5f;
                  T->vmax = vmid, B->vmin = vmid;
                  f.insert(T);
                  f.insert(B);
            }
      }
}

bool NuPatch::diceable(MicroPolygonGrid &g, Hider &h, bool &usplit, bool &vsplit)
{
      // TODO: make TESTSIZE user option
      #define TESTSIZE 8

      MicroPolygonGrid testGrid;
      testGrid.setDim(TESTSIZE, TESTSIZE, this);
      dice(testGrid);
      RtPoint* P_grid = (RtPoint*)testGrid.findVariable("P");
      float maxUDist=-1e10f, maxVDist=-1e10f;
      for (int vg=0; vg<TESTSIZE; ++vg) {
            const int Tvg = (TESTSIZE+1)*vg;
            for (int ug=0; ug<TESTSIZE; ++ug) {
                  const int ugTvg = ug + Tvg;
                  Point3 p1(P_grid[ugTvg][0], P_grid[ugTvg][1], P_grid[ugTvg][2]);
                  maxUDist = MAX2(h.rasterEstimate(p1, Point3(P_grid[ugTvg + 1][0], P_grid[ugTvg + 1][1], P_grid[ugTvg + 1][2])), maxUDist);
                  maxVDist = MAX2(h.rasterEstimate(p1, Point3(P_grid[ugTvg + (TESTSIZE+1)][0], P_grid[ugTvg + (TESTSIZE+1)][1], P_grid[ugTvg + (TESTSIZE+1)][2])), maxVDist);
            }
      }

      const Attributes* attr = getAttributeReference();
      // 2x2 mp minimum
      unsigned int xdim = (attr->flags & AF_DICE_BINARY) ? (1 << MAX2(1, int(0.5f + (float)log(MAX2(1e-7f, TESTSIZE*maxUDist) / attr->effectiveShadingRate)*(float)M_LOG2E)))
                                                         : MAX2(2, int(0.5f + (TESTSIZE*maxUDist) / attr->effectiveShadingRate));
      unsigned int ydim = (attr->flags & AF_DICE_BINARY) ? (1 << MAX2(1, int(0.5f + (float)log(MAX2(1e-7f, TESTSIZE*maxVDist) / attr->effectiveShadingRate)*(float)M_LOG2E)))
                                                         : MAX2(2, int(0.5f + (TESTSIZE*maxVDist) / attr->effectiveShadingRate));

      g.setDim(xdim, ydim, this);
      if ((xdim<h.gridSize()) && (ydim<h.gridSize()) && ((xdim*ydim)<h.gridSize())) return true;
      usplit = (xdim > ydim);
      vsplit = !usplit;
      return false;
      #undef TESTSIZE
}

void NuPatch::dice(MicroPolygonGrid &g, bool Pclose)
{
      int ug, vg, idx = 0;
      float u, v;
      RtPoint* P_grid = (RtPoint*)(Pclose ? g.addVariable("=Pclose") : g.addVariable("P"));
      RtVector *dPdu_grid = (RtVector*)g.addVariable("dPdu"), *dPdv_grid = (RtVector*)g.addVariable("dPdv");
      RtFloat *u_grid = g.addVariable("u", 1), *v_grid = g.addVariable("v", 1);
      const int xdim = g.get_xdim(), ydim = g.get_ydim();

      // Dice P, dPdu & dPdv
      // note: du/dv adjusted for clipping,
      // this will result in minute precision errors that make u and/or v
      // less or greater than the valid u/v min/max range causing lots of chaos...
      // To prevent this, even though it happens less (or not at all) with just using 0-1 range, eval() clamps the values
      const float du = (umax - umin)/(float)xdim, dv = (vmax - vmin)/(float)ydim;
      for (vg=0, v=vmin; vg<=ydim; ++vg, v+=dv)
            for (ug=0, u=umin; ug<=xdim; ++ug, u+=du, ++idx) {
                  eval(u, v, P_grid[idx], dPdu_grid[idx], dPdv_grid[idx]);
                  u_grid[idx] = u, v_grid[idx] = v;
            }

      // does not deal with 'vertex' type yet (here does same as varying), TODO

      // dice other primvars, if any
      g.initVars(std_dice);
      Primitive::linear_dice(g);
}

Primitive* NuPatch::duplicateSelf()
{
      NuPatch* np = new NuPatch();
      memcpy(np, this, sizeof(NuPatch));
      if (nudata) nudata->refc++;   // after copying, not before!!!
      if (primvars) primvars->incRefCount();
      return np;
}

__END_QDRENDER