[Bf-blender-cvs] [6103835] soc-2014-nurbs: BSpline basis (derivs 0, 1, 2), BSpline surface (derivs 0, 1) fix + verification in Mathematica
Jonathan deWerd
noreply at git.blender.org
Sat Jul 26 20:47:53 CEST 2014
Commit: 6103835784c456d1fc7d0bfb2ae6d6dde8335763
Author: Jonathan deWerd
Date: Mon Jul 21 23:44:13 2014 -0400
Branches: soc-2014-nurbs
https://developer.blender.org/rB6103835784c456d1fc7d0bfb2ae6d6dde8335763
BSpline basis (derivs 0,1,2), BSpline surface (derivs 0,1) fix + verification in Mathematica
===================================================================
M source/blender/blenkernel/BKE_curve.h
M source/blender/blenkernel/intern/curve.cpp
M source/blender/blenkernel/intern/curve_eval.cpp
M tests/interactive/curve_eval.h
M tests/interactive/nurbs_derivative_eval.cpp
M tests/interactive/nurbs_eval.nb
===================================================================
diff --git a/source/blender/blenkernel/BKE_curve.h b/source/blender/blenkernel/BKE_curve.h
index 2268cc6..378e3c0 100644
--- a/source/blender/blenkernel/BKE_curve.h
+++ b/source/blender/blenkernel/BKE_curve.h
@@ -123,8 +123,8 @@ void BKE_curve_forward_diff_bezier(float q0, float q1, float q2, float q3, float
void BKE_curve_rect_from_textbox(const struct Curve *cu, const struct TextBox *tb, struct rctf *r_rect);
/* ** Nurbs ** */
-int BKE_nurbs_nz_basis_range(float u, float *knots, int num_knots, int order);
-void BKE_nurbs_basis_eval(float u, int i, float *U, int num_knots, int order, int nd, float out[][NURBS_MAX_ORDER]);
+int BKE_bspline_nz_basis_range(float u, float *knots, int num_pts, int order);
+void BKE_bspline_basis_eval(float u, int i, float *U, int num_pts, int order, int nd, float out[][NURBS_MAX_ORDER]);
bool BKE_nurbList_index_get_co(struct ListBase *editnurb, const int index, float r_co[3]);
diff --git a/source/blender/blenkernel/intern/curve.cpp b/source/blender/blenkernel/intern/curve.cpp
index b14bdd7..7c6d559 100644
--- a/source/blender/blenkernel/intern/curve.cpp
+++ b/source/blender/blenkernel/intern/curve.cpp
@@ -4134,8 +4134,8 @@ void BKE_nurb_make_displist(struct Nurb *nu, struct DispList *dl) {
// Right now xyz = {u,v,0.0}. Apply the NURBS map to transform it to {x,y,z}
if (xyz[1]!=cached_v) {
/* only N(i-p)p, N(i-p+1)p, ..., Nip are nz*/
- iv = BKE_nurbs_nz_basis_range(xyz[1], nu->knotsv, KNOTSV(nu), orderv);
- BKE_nurbs_basis_eval(xyz[1], iv, nu->knotsv, KNOTSV(nu), orderv, 1, basisv);
+ iv = BKE_bspline_nz_basis_range(xyz[1], nu->knotsv, nu->pntsv, orderv);
+ BKE_bspline_basis_eval(xyz[1], iv, nu->knotsv, nu->pntsv, orderv, 1, basisv);
if (xyz[0]>0.75) {
printf("vcoeff(%f): ",xyz[1]);
for (int l=0; l<orderv; l++)
@@ -4143,8 +4143,8 @@ void BKE_nurb_make_displist(struct Nurb *nu, struct DispList *dl) {
}
cached_v = xyz[1];
}
- iu = BKE_nurbs_nz_basis_range(xyz[0], nu->knotsu, KNOTSU(nu), orderu);
- BKE_nurbs_basis_eval(xyz[0], iu, nu->knotsu, KNOTSU(nu), orderu, 1, basisu);
+ iu = BKE_bspline_nz_basis_range(xyz[0], nu->knotsu, nu->pntsu, orderu);
+ BKE_bspline_basis_eval(xyz[0], iu, nu->knotsu, nu->pntsu, orderu, 1, basisu);
if (xyz[0]>0.75) {
printf("ucoeff(%f): ",xyz[0]);
for (int l=0; l<orderu; l++)
diff --git a/source/blender/blenkernel/intern/curve_eval.cpp b/source/blender/blenkernel/intern/curve_eval.cpp
index 7e5099f..7c883c9 100644
--- a/source/blender/blenkernel/intern/curve_eval.cpp
+++ b/source/blender/blenkernel/intern/curve_eval.cpp
@@ -40,11 +40,36 @@ extern "C" {
}
#endif
+int BKE_nurb_binom_coeffs[NURBS_MAX_ORDER][NURBS_MAX_ORDER+1];
+int BKE_nurb_binom_coeff(int k, int i) {
+ static int inited=0;
+ if (!inited) {
+ inited=1;
+ for (int kk=0; kk<NURBS_MAX_ORDER; kk++) {
+ int last = 1;
+ BKE_nurb_binom_coeffs[kk][0] = last;
+ for (int ii=1; ii<=kk; ii++) {
+ last *= kk-ii+1;
+ last /= i;
+ BKE_nurb_binom_coeffs[kk][ii] = last;
+ }
+ }
+ }
+#ifndef NDEBUG
+ if (k>=NURBS_MAX_ORDER || i>k) {
+ fprintf(stderr, "Invalid binomial coeff %i-choose-%i.\n",k,i);
+ return 0;
+ }
+#endif
+ return BKE_nurb_binom_coeffs[k][i];
+}
+
+
/* ~~~~~~~~~~~~~~~~~~~~Non Uniform Rational B Spline calculations ~~~~~~~~~~~ */
float nurbs_eps = 1e-5;
-void BKE_nurb_knot_calc(int flags, int pnts, int order, float knots[])
+void BKE_bspline_knot_calc(int flags, int pnts, int order, float knots[])
{
int num_knots = pnts + order;
if (flags & CU_NURB_CYCLIC) num_knots += order-1;
@@ -113,11 +138,11 @@ void BKE_nurb_knot_calc(int flags, int pnts, int order, float knots[])
*/
}
-/* Points on the surface of a NURBS curve are defined as an affine combination
+/* Points on a NURBS surface 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
+ * control points have nonzero weight 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
@@ -128,13 +153,11 @@ void BKE_nurb_knot_calc(int flags, int pnts, int order, float knots[])
* 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
*/
-int BKE_nurbs_nz_basis_range(float u, float *knots, int num_knots, int order)
+int BKE_bspline_nz_basis_range(float u, float *knots, int num_pts, int order)
{
int p = order-1; /* p is the NURBS degree */
- int m = num_knots-1; /* 0-based index of the last knot */
- int n = m-p-1; /* = number of control points + 1 */
-
- /*
+ int n = num_pts-1;
+ // Binary Search
if (u>=knots[n+1])
return n;
if (u<=knots[p])
@@ -145,27 +168,45 @@ int BKE_nurbs_nz_basis_range(float u, float *knots, int num_knots, int order)
else low=mid;
mid = (low+high)/2;
}
+#ifndef NDEBUG
+ // Linear Search
+ int ls_ret;
+ for (int i=p; i<=n; i++) {
+ if (knots[i]<=u && u<knots[i+1]) {
+ ls_ret = i;
+ break;
+ }
+ }
+ if (mid!=ls_ret) {
+ fprintf(stderr, "BS%i[%f,%f) LS%i[%f,%f) MISMATCH!\n",mid,knots[mid],knots[mid+1], ls_ret,knots[ls_ret],knots[ls_ret+1]);
+ }
+ if (!(p<=ls_ret && ls_ret<=n+1)) {
+ fprintf(stderr,"NURBS tess: knot index 2 out of bounds\n");
+ }
+#endif
return mid;
- */
-
}
/* Computes the p+1=order nonvanishing NURBS basis functions at coordinate u.
* Arguments:
* u = the curve parameter, as above (u is actually v if uv==2)
- * i = the return value of nurbs_nz_basis_range(u,...)
+ * i = such that u_i<=u<u_i+1 = the return value of nurbs_nz_basis_range(u,...)
* U = the knot vector
- * num_knots = number of knots in U
- * 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
+ * num_pts = number of knots in U
+ * order = the order of the NURBS curve
+ * 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
+ * Let N(i,p,u) stand for the ith basis func of order p eval'd at u.
+ * Let N(i) be a shorthand when p and u are implicit from context.
+ * Efficient algorithm adapted from Piegl&Tiller 1995
*/
-void BKE_nurbs_basis_eval(float u, int i, float *U, int num_knots, int order, int nd, float out[][NURBS_MAX_ORDER]) {
+void BKE_bspline_basis_eval(float u, int i, float *U, int num_pts, int order, int nd, float out[][NURBS_MAX_ORDER]) {
+ int num_knots = num_pts + order;
int p = order-1; /* p is the NURBS degree */
int n = num_knots-p; /* = number of control points + 1 */
double left[NURBS_MAX_ORDER], right[NURBS_MAX_ORDER];
@@ -173,35 +214,48 @@ void BKE_nurbs_basis_eval(float u, int i, float *U, int num_knots, int order, in
double a[2][NURBS_MAX_ORDER];
double saved,temp;
if (u<U[p]) u=U[p];
- if (u>U[num_knots-p-2]) u=U[num_knots-p-2];
+ if (u>U[num_knots-p-1]) u=U[num_knots-p-1];
if (i==-1) {
- i = BKE_nurbs_nz_basis_range(u, U, num_knots, order);
+ i = BKE_bspline_nz_basis_range(u, U, num_pts, order);
+ }
+#ifndef NDEBUG
+ if (!(U[i]<=u && u<=U[i+1])) {
+ fprintf(stderr, "NURBS tess error: curve argument out of bounds\n");
+ return;
+ }
+ if (!(p<=i && i<=num_knots-p-1 && i<=10)) {
+ fprintf(stderr, "NURBS tess error: knot index i out of bounds\n");
+ return;
}
+#endif
/* First, compute the 0th derivatives of the basis functions. */
ndu[0][0] = 1.0;
for (int j=1; j<=p; j++) {
+ /* Compute degree j basis functions eval'd at u */
+ /* N(i-j,j,u) N(i-j+1,j,u) ... N(i,j,u) */
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++) {
+ /* Compute N(i-j+r, j, u) */
ndu[j][r] = right[r+1]+left[j-r];
- printf("ndu[%i][%i]=%f for i=%i\n",j,r,ndu[j][r],i);
temp = ndu[r][j-1]/ndu[j][r];
- ndu[r][j] = saved+right[r+1]*temp;
+ ndu[r][j] = saved+right[r+1]*temp; /* N(i-p+r,j,u) */
saved = left[j-r]*temp;
}
- ndu[j][j] = saved;
+ /* Edge case: prev loop didn't store N(i-j+j,j,u)==N(i,j,u) */
+ ndu[j][j] = saved; /* N(i,j,u) */
}
for (int j=0; j<=p; j++)
out[0][j] = ndu[j][p];
- return;
+
/* 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++) {
+ /* compute kth derivative */
double d = 0.0;
int rk=r-k, pk=p-k, j=0;
if (r>=k) {
@@ -222,12 +276,114 @@ void BKE_nurbs_basis_eval(float u, int i, float *U, int num_knots, int order, in
j=s1; s1=s2; s2=j;
}
}
+ /* fix the normalization of derivatives */
int r=p;
- for (int k=1; k<=n; k++) {
+ for (int k=1; k<=n && k<=nd; k++) {
for (int j=0; j<=p; j++) {
out[k][j] *= r;
}
- r *=(p-k);
+ r *= (p-k);
}
}
+/* Evaluates a B-Spline curve (and nd of its derivatives) at point u.
+ * NOTE: *DOES NOT* perform perspective divide on output 4-vectors. Some algorithms
+ * need the information stored in the weights.
+ *
+ * u: the point to eval the curve at
+ * U: the knot vector
+ * num_pts: number of control points
+ * order: order of the NURBS curve
+ * P,stride: P[0], P[stride], ..., P[(num_pts-1)*stride] are the ctrl points
+ * nd: number of derivatives to calculate
+ * out: out[0], out[1], ..., out[nd] are the evaluated points/derivatives
+ */
+void BKE_bspline_curve_eval(float u, float *U, int num_pts, int order, BPoint *P, int stride, int nd, BPoint *out)
+{
+ int p = order-1;
+ if (!(
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list