[Bf-blender-cvs] [b40acb2] soc-2014-nurbs: Added new eval code (computes first and second surface partials). Still buggy near multi-knots.
Jonathan deWerd
noreply at git.blender.org
Sat Jul 19 12:25:55 CEST 2014
Commit: b40acb2d98f2d95c7d2300b966798cf40b7b9eca
Author: Jonathan deWerd
Date: Thu Jul 17 11:25:23 2014 -0400
https://developer.blender.org/rBb40acb2d98f2d95c7d2300b966798cf40b7b9eca
Added new eval code (computes first and second surface partials). Still buggy near multi-knots.
===================================================================
M source/blender/blenkernel/BKE_curve.h
M source/blender/blenkernel/intern/curve.cpp
M source/blender/editors/curve/editcurve_add.c
M source/blender/makesdna/DNA_curve_types.h
===================================================================
diff --git a/source/blender/blenkernel/BKE_curve.h b/source/blender/blenkernel/BKE_curve.h
index e40ed00..8420b86 100644
--- a/source/blender/blenkernel/BKE_curve.h
+++ b/source/blender/blenkernel/BKE_curve.h
@@ -53,6 +53,8 @@ typedef struct CurveCache {
struct Path *path;
} CurveCache;
+#define NURBS_MAX_ORDER 10
+
#define KNOTSU(nu) ( (nu)->orderu + (nu)->pntsu + (((nu)->flagu & CU_NURB_CYCLIC) ? ((nu)->orderu - 1) : 0) )
#define KNOTSV(nu) ( (nu)->orderv + (nu)->pntsv + (((nu)->flagv & CU_NURB_CYCLIC) ? ((nu)->orderv - 1) : 0) )
diff --git a/source/blender/blenkernel/intern/curve.cpp b/source/blender/blenkernel/intern/curve.cpp
index 7b51ff5..26faac0 100644
--- a/source/blender/blenkernel/intern/curve.cpp
+++ b/source/blender/blenkernel/intern/curve.cpp
@@ -872,132 +872,226 @@ void BKE_nurb_bezt_calc_plane(struct Nurb *nu, struct BezTriple *bezt, float r_p
}
/* ~~~~~~~~~~~~~~~~~~~~Non Uniform Rational B Spline calculations ~~~~~~~~~~~ */
+float nurbs_eps = 1e-5;
-
-static void calcknots(float *knots, const int pnts, const short order, const short flag)
+/* uv: 1=u 2=v */
+static void makeknots(Nurb *nu, short uv)
{
- /* knots: number of pnts NOT corrected for cyclic */
- const int pnts_order = pnts + order;
- float k;
- int a;
-
- switch (flag & (CU_NURB_ENDPOINT | CU_NURB_BEZIER)) {
- case CU_NURB_ENDPOINT:
- k = 0.0;
- for (a = 1; a <= pnts_order; a++) {
- knots[a - 1] = k;
- if (a >= order && a <= pnts)
- k += 1.0f;
- }
- break;
- case CU_NURB_BEZIER:
- /* Warning, the order MUST be 2 or 4,
- * if this is not enforced, the displist will be corrupt */
- if (order == 4) {
- k = 0.34;
- for (a = 0; a < pnts_order; a++) {
- knots[a] = floorf(k);
- k += (1.0f / 3.0f);
- }
- }
- else if (order == 3) {
- k = 0.6f;
- for (a = 0; a < pnts_order; a++) {
- if (a >= order && a <= pnts)
- k += 0.5f;
- knots[a] = floorf(k);
- }
- }
- else {
- printf("bez nurb curve order is not 3 or 4, should never happen\n");
- }
- break;
- default:
- for (a = 0; a < pnts_order; a++) {
- knots[a] = (float)a;
+ bool u = uv==1;
+ int flags=0, pnts=0, order=0;
+ if (u) {
+ if (!BKE_nurb_check_valid_u(nu)) return;
+ if (nu->knotsu) MEM_freeN(nu->knotsu);
+ flags = nu->flagu;
+ pnts = nu->pntsu;
+ order = nu->orderu;
+ } else {
+ if (!BKE_nurb_check_valid_v(nu)) return;
+ if (nu->knotsv) MEM_freeN(nu->knotsv);
+ flags = nu->flagv;
+ pnts = nu->pntsv;
+ order = nu->orderv;
+ }
+ printf("p:%i o:%i\n",pnts,order);
+ int num_knots = pnts + order;
+ if (flags & CU_NURB_CYCLIC) num_knots += order-1;
+ float *knots = (float*)MEM_callocN(sizeof(float)*num_knots, "makeknots");
+ if (flags & CU_NURB_BEZIER) {
+ /* Previous versions of blender supported "Bezier" knot vectors. These
+ * are useless to the user. *All* NURBS surfaces can be transformed into
+ * Bezier quilts (grids of connected Bezier surfaces). The "Bezier" knot
+ * vector is only an intermediate mathematical step in this process.
+ * Bezier quilts may be exposed to the user but there is no point in having
+ * the option to use Bezier-style knot vectors without Bezier-style controls.
+ * |------------------ pnts+order -------------|
+ * |--ord--||-ord-1-||-ord-1-||-ord-1-||--ord--|
+ * Bezier: { 0 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 4 }
+ */
+ int v=0;
+ for (int i=0; i<order; i++)
+ knots[i] = v;
+ for (int i=order,reps=0; i<pnts; i++) {
+ if (reps==0) {
+ v += 1;
+ reps = order-1;
}
- break;
+ knots[i] = v;
+ }
+ for (int i=pnts; i<pnts+order; i++)
+ knots[i] = v;
+ } else if (flags & CU_NURB_ENDPOINT) {
+ /* "Endpoint" knot vectors ensure that the first and last knots are
+ * repeated $order times so that the valid NURBS domain actually touches
+ * the edges of the control polygon. These are the default.
+ * |------- pnts+order ------------|
+ * |-order-||--pnts-order-||-order-|
+ * { 0 0 0 0 1 2 3 4 5 6 7 8 9 9 9 9 }
+ */
+ for (int i=0; i<order; i++)
+ knots[i] = 0;
+ int v=1;
+ for (int i=order; i<pnts; i++)
+ knots[i] = v++;
+ for (int i=pnts; i<pnts+order; i++)
+ knots[i] = v;
+ } else if ((flags&CU_NURB_CYCLIC) || true) {
+ /* Uniform knot vectors are the simplest mathematically but they create
+ * the annoying problem where the edges of the surface do not reach the
+ * edges of the control polygon which makes positioning them difficult.
+ * |------ pnts+order --------|
+ * |-----pnts----||---order---|
+ * { 0 1 2 3 4 5 6 7 8 9 10 11 }
+ */
+ /* Cyclic knot vectors are equivalent to uniform knot vectors over
+ * a control polygon of pnts+(order-1) points for a total of
+ * pnts+(order-1)+order knots. The additional (order-1) control points
+ * are repeats of the first (order-1) control points. This guarantees
+ * that all derivatives match at the join point.
+ * |----- (pnts+order-1)+order --------|
+ * |-----pnts----||--reps-||---order---|
+ * { 0 1 2 3 4 5 6 7 8 9 10 11 12 13 }
+ */
+ for (int i=0; i<num_knots; i++)
+ knots[i]=i;
}
+ /*
+ printf("Knots%s: ",u?"u":"v");
+ for (int i=0; i<num_knots; i++)
+ printf((i!=num_knots-1)?"%f, ":"%f\n", knots[i]);
+ */
+
+ if (u)
+ nu->knotsu=knots;
+ else
+ nu->knotsv=knots;
}
-static void makecyclicknots(float *knots, int pnts, short order)
-/* pnts, order: number of pnts NOT corrected for cyclic */
+void BKE_nurb_knot_calc_u(Nurb *nu)
{
- int a, b, order2, c;
-
- if (knots == NULL)
- return;
-
- order2 = order - 1;
-
- /* do first long rows (order -1), remove identical knots at endpoints */
- if (order > 2) {
- b = pnts + order2;
- for (a = 1; a < order2; a++) {
- if (knots[b] != knots[b - a])
- break;
- }
- if (a == order2)
- knots[pnts + order - 2] += 1.0f;
- }
-
- b = order;
- c = pnts + order + order2;
- for (a = pnts + order2; a < c; a++) {
- knots[a] = knots[a - 1] + (knots[b] - knots[b - 1]);
- b--;
- }
+ makeknots(nu, 1);
}
+void BKE_nurb_knot_calc_v(Nurb *nu)
+{
+ makeknots(nu, 2);
+}
-
-static void makeknots(Nurb *nu, short uv)
+/* Points on the surface of a NURBS curve are defined as an affine combination
+ * (sum of coeffs is 1) of points from the control polygon. HOWEVER, the
+ * locality property of NURBS dictates that at most $order consecutive
+ * control points are nonzero at a given curve coordinate u (or v, but we assume
+ * u WLOG for the purposes of naming the variable). Therefore
+ * C(u) = \sum_{j=0}^n Njp(u)*Pj
+ * = \sum_{j=i-p}^j Njp(u)*Pj
+ * = N(i-p)p*P(i-p) + N(i-p+1)p*P(i-p+1) + ... + Nip*Pi
+ * for u in knot range [u_i, u_{i+1}) given the m+1 knots
+ * U[] = {u0, u1, ..., um}
+ * Arguments:
+ * uv = 1:u, 2:v
+ * u = the curve parameter, as above (u is actually v if uv==2)
+ * returns: i such that basis functions i-p,i-p+1,...,i are possibly nonzero
+ */
+static int nurbs_nz_basis_range(Nurb *nu, int uv, float u)
{
- if (nu->type == CU_NURBS) {
- if (uv == 1) {
- if (nu->knotsu)
- MEM_freeN(nu->knotsu);
- if (BKE_nurb_check_valid_u(nu)) {
- nu->knotsu = (float*)MEM_callocN(4 + sizeof(float) * KNOTSU(nu), "makeknots");
- if (nu->flagu & CU_NURB_CYCLIC) {
- calcknots(nu->knotsu, nu->pntsu, nu->orderu, 0); /* cyclic should be uniform */
- makecyclicknots(nu->knotsu, nu->pntsu, nu->orderu);
- }
- else {
- calcknots(nu->knotsu, nu->pntsu, nu->orderu, nu->flagu);
- }
+ int num_knots = (uv==1)? KNOTSU(nu) : KNOTSV(nu);
+ float *knots = (uv==1)? nu->knotsu : nu->knotsv;
+ int order = (uv==1)? nu->orderu : nu->orderv;
+ int p = order-1; /* p is the NURBS degree */
+ int n = num_knots-p; /* = number of control points + 1 */
+
+ if (u>=knots[n+1])
+ return n;
+ if (u<=knots[p])
+ return p;
+ int low=p, high=n+1, mid=(low+high)/2;
+ while (u<knots[mid] || u>=knots[mid+1]) {
+ if (u<knots[mid]) high=mid;
+ else low=mid;
+ mid = (low+high)/2;
+ }
+ return mid; /* called i in nurbs_basis_eval */
+}
+
+/* Computes the p+1=order nonvanishing NURBS basis functions at coordinate u.
+ * Arguments:
+ * uv = 1:u, 2:v
+ * i = the return value of nurbs_nz_basis_range or -1 to compute automatically
+ * u = the curve parameter, as above (u is actually v if uv==2)
+ * nd= the number of derivatives to calculate (0 = just regular basis funcs)
+ * out = an array to put N(i-p),N(i-p+1),...,N(i) and their derivatives into
+ * out[0][0]=N(i-p) out[0][1]=N(i-p+1) ... out[0][p]=N(i)
+ * out[1][0]=N'(i-p) out[1][1]=N'(i-p+1) ... out[1][p]=N'(i)
+ * ...
+ * out[nd-1][0]=N'''(i-p) out[nd-1][1]=N'''(i-p+1) ... out[nd-1][p]=N'''(i)
+ * ^ let ''' stand for differentiation nd-1 times
+ * Adapted from Piegl&Tiller 1995
+ */
+static void nurbs_basis_eval(Nurb *nu, int uv, float u, int i, float out[][NURBS_MAX_ORDER], int nd) {
+ if (.99<u && u<1.01) {
+ i=4;
+ }
+ int num_knots = (uv==1)? KNOTSU(nu) : KNOTSV(nu);
+ float *U = (uv==1)? nu->knotsu : nu->knotsv;
+ int order = (uv==1)? nu->orderu : nu->orderv;
+ int p = order-1; /* p is the NURBS degree */
+ int n = num_knots-p; /* = number of control points + 1 */
+ if (i==-1) /* index st N(i-p),N(i-p+1),...,N(i) are nonzero */
+ i = nurbs_nz_basis_range(nu, uv, u);
+ double left[NURBS_MAX_ORDER], right[NURBS_MAX_ORDER];
+ double ndu[NURBS_MAX_ORDER][NURBS_MAX_ORDER];
+ double a[2][NURBS_MAX_ORDER];
+ double saved,temp;
+
+ /* First, compute the 0th derivatives of the basis functions. */
+ ndu[0][0] = 1.0;
+ for (int j=1; j<=p; j++) {
+ left[j] = u-U[i+1-j];
+ right[j] = U[i+j]-u;
+ saved = 0;
+ /* Upper and lower triangles of Piegl&Tiller eval grid */
+ for (int r=0; r<j; r++) {
+ ndu[j][r] = right[r+1]+left[j-r];
+ temp = ndu[r][j-1]/ndu[j][r];
+ ndu[r][j] = saved+right[r+1]*temp;
+ saved = left[j-r]*temp;
+ }
+ ndu[j][j] = saved;
+ }
+ for (int j=0; j<=p; j++)
+ out[0][j] = ndu[j][p];
+
+ /* Now compute the higher nd derivatives */
+ for (int r=0; r<=p; r++) {
+ int s1=0, s2=1, j1=0, j2=0;
+ a[0][0] = 1.0;
+ for (int k=1; k<=nd && k<=n; k++) {
+ double d = 0.0;
+ int rk=r-k, pk=p-k, j=0;
+ if (r>=k) {
+ a[s2][0] = a[s1][0]/ndu[pk+1][rk];
+ d = a[s2][0]*ndu[
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list