[Bf-blender-cvs] [f8afdc971f4] temp-sculpt-roll-mapping: temp-sculpt-roll-mappings: Experimental quadratic spline code
Joseph Eagar
noreply at git.blender.org
Wed Nov 16 09:45:20 CET 2022
Commit: f8afdc971f47f2a4fd6bd94fd9df558ee87f6fe6
Author: Joseph Eagar
Date: Wed Nov 16 00:40:46 2022 -0800
Branches: temp-sculpt-roll-mapping
https://developer.blender.org/rBf8afdc971f47f2a4fd6bd94fd9df558ee87f6fe6
temp-sculpt-roll-mappings: Experimental quadratic spline code
Doesn't quite work yet. The idea is to directly solve for
curve texture coordinates using a quadratic spline.
===================================================================
M source/blender/blenlib/BLI_even_spline.hh
M source/blender/editors/sculpt_paint/paint_intern.h
M source/blender/editors/sculpt_paint/paint_stroke.cc
===================================================================
diff --git a/source/blender/blenlib/BLI_even_spline.hh b/source/blender/blenlib/BLI_even_spline.hh
index 4b80df0202c..6d6e547ef4c 100644
--- a/source/blender/blenlib/BLI_even_spline.hh
+++ b/source/blender/blenlib/BLI_even_spline.hh
@@ -7,7 +7,9 @@
#include "BLI_math_vec_types.hh"
#include "BLI_vector.hh"
+#include <algorithm>
#include <cstdio>
+#include <type_traits>
#include <utility>
//#define FINITE_DIFF
@@ -34,6 +36,755 @@ template<typename Float> class Curve {
};
*/
+/** Quadratic curves */
+
+/*
+comment: Reduce algebra script;
+
+on factor;
+off period;
+
+procedure bez(a, b);
+ a + (b - a) * t;
+
+lin := bez(k1, k2);
+quad := bez(lin, sub(k2=k3, k1=k2, lin));
+
+dquad := df(quad, t);
+iquad := int(quad, t);
+
+x1 := 0;
+y1 := 0;
+
+dx := sub(k1=x1, k2=x2, k3=x3, dquad);
+dy := sub(k1=y1, k2=y2, k3=y3, dquad);
+darc := sqrt(dx**2 + dy**2);
+
+arcstep := darc*dt + 0.5*df(darc, t)*dt*dt;
+
+d2x := df(dx / darc, t);
+d2y := df(dy / darc, t);
+
+gentran
+begin
+declare <<
+x1,x2,x3,x4 : float;
+y1,y2,y3,y4 : float;
+dt,t : float;
+>>;
+return eval(arcstep)
+end;
+
+on fort;
+quad;
+dquad;
+iquad;
+arcstep;
+d2x;
+d2y;
+off fort;
+
+*/
+template<typename Float, int axes = 2, int table_size = 512> class QuadBezier {
+ using Vector = vec_base<Float, axes>;
+
+ public:
+ Vector ps[3];
+
+ static const int TableSize = table_size;
+
+ QuadBezier(Vector a, Vector b, Vector c)
+ {
+ ps[0] = a;
+ ps[1] = b;
+ ps[2] = c;
+
+ deleted = false;
+ _arc_to_t = new Float[table_size];
+ _t_to_arc = new Float[table_size];
+ }
+
+ ~QuadBezier()
+ {
+ deleted = true;
+
+ if (_arc_to_t) {
+ delete[] _arc_to_t;
+ _arc_to_t = nullptr;
+ }
+
+ if (_t_to_arc) {
+ delete[] _t_to_arc;
+ _t_to_arc = nullptr;
+ }
+ }
+
+ QuadBezier()
+ {
+ deleted = false;
+ _arc_to_t = new Float[table_size];
+ _t_to_arc = new Float[table_size];
+ }
+
+ QuadBezier(const QuadBezier &b)
+ {
+ _arc_to_t = new Float[table_size];
+ _t_to_arc = new Float[table_size];
+
+ *this = b;
+ deleted = false;
+ }
+
+ QuadBezier &operator=(const QuadBezier &b)
+ {
+ ps[0] = b.ps[0];
+ ps[1] = b.ps[1];
+ ps[2] = b.ps[2];
+
+ length = b.length;
+
+ if (!_arc_to_t) {
+ _arc_to_t = new Float[table_size];
+ }
+ if (!_arc_to_t) {
+ _t_to_arc = new Float[table_size];
+ }
+
+ if (b._arc_to_t) {
+ for (int i = 0; i < table_size; i++) {
+ _arc_to_t[i] = b._arc_to_t[i];
+ }
+ }
+
+ if (b._t_to_arc) {
+ for (int i = 0; i < table_size; i++) {
+ _t_to_arc[i] = b._t_to_arc[i];
+ }
+ }
+
+ return *this;
+ }
+
+#if 1
+ QuadBezier(QuadBezier &&b)
+ {
+ *this = b;
+ }
+
+ Vector closest_point(const Vector p, Float &r_s, Vector &r_tan, Float &r_dis)
+ {
+# if 0
+ r_dis = 0.01;
+ r_s = 0.0;
+ r_tan = Vector(0.0, 1.0, 0.0);
+
+ return p;
+# endif
+
+# if 0
+ const int steps = 16;
+ Float s = 0.0, ds = length / (steps - 1);
+ Vector retco;
+
+ r_dis = FLT_MAX;
+
+ for (int i = 0; i < steps; i++, s += ds) {
+ Vector co = evaluate(s);
+ Vector vec = co - p;
+
+ Float dist = _dot(vec, vec);
+
+ if (dist < r_dis) {
+ retco = co;
+ r_dis = dist;
+ r_s = s;
+ r_tan = derivative(s);
+ }
+ }
+
+ r_dis = std::sqrt(r_dis);
+
+ return retco;
+#elif 0
+ float weights[3];
+ float n[3];
+
+ //normal_tri_v3(ps[0], ps[1], ps[2]);
+ float p2[3];
+
+ closest_on_tri_to_point_v3(p2, p, ps[0], ps[1], ps[2]);
+
+ //barycentric_weights(ps[0], ps[1], ps[2], p, n, weights);
+ resolve_tri_uv_v3(weights, p2, ps[0], ps[1], ps[2]);
+
+ Float u = weights[0];
+ Float v = weights[1];
+
+ Float t1 = v / (2 * u + v);
+ Float t2 = (2 * (v - 1 + u)) / (2 * u + v - 2);
+
+ //v = min_ff(max_ff(v, 0.0f), 1.0f);
+
+ //v = std::sqrt(v);
+
+
+ t1 = min_ff(max_ff(t1, 0.0), 1.0);
+ t2 = min_ff(max_ff(t2, 0.0), 1.0);
+
+ // t1 = std::min(std::max(t1, 0.0), 1.0);
+ // t2 = std::min(std::max(t2, 0.0), 1.0);
+
+ //t1 = t2 = 0.5;
+
+ Float s1 = t_to_arc(t1);
+ Float s2 = t_to_arc(t2);
+
+ Vector a = evaluate(s1);
+ Vector b = evaluate(s2);
+
+ a -= p;
+ b -= p;
+
+ Float dis1 = _dot(a, a);
+ Float dis2 = _dot(b, b);
+
+ if (dis1 < dis2) {
+ r_s = s1;
+ r_dis = std::sqrt(dis1);
+ r_tan = derivative(s1);
+
+ return a + p;
+ }
+ else {
+ r_s = s2;
+ r_dis = std::sqrt(dis2);
+ r_tan = derivative(s2);
+
+ return b + p;
+ }
+ //r_dis = v;
+
+ //closest_on_tri_to_point_v3
+ //resolve_tri_uv_v3
+# else
+ /*
+ on factor;
+ off period;
+
+ procedure bez(a, b);
+ a + (b - a) * t;
+
+ lin := bez(k1, k2);
+ quad := bez(lin, sub(k2=k3, k1=k2, lin));
+
+ x := sub(k1=x1, k2=x2, k3=x3, quad);
+ y := sub(k1=y1, k2=y2, k3=y3, quad);
+ z := sub(k1=z1, k2=z2, k3=z3, quad);
+
+ comment: origin is at x1, y1, z1;
+
+ x1 := 0;
+ y1 := 0;
+ z1 := 0;
+
+ comment: rotate to align with control triangle normal;
+ comment: z2 := 0;
+ comment: z3 := 0;
+
+ dx := df(x, t);
+ dy := df(y, t);
+ dz := df(z, t);
+
+ tx := x - px;
+ ty := y - py;
+ tz := z - pz;
+
+ x1 := 0.0;
+ y1 := 0.0;
+ z1 := 0.0;
+
+ x2 := 0.5;
+ y2 := 0.0;
+ z2 := 0.0;
+
+ x3 := 0.5;
+ y3 := 0.5;
+ z3 := 0.0;
+
+ f := df(tx**2 + ty**2 + tz**2, t, 1);
+ f := tx*dx + ty*dy + tz*dz;
+
+ ff := solve(f, t);
+
+ on fort;
+ part(ff, 1, 2);
+ part(ff, 2, 2);
+ off fort;
+
+ Computer Aided Geometric Design
+ Thomas W. Sederberg
+ page 203
+ https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=1000&context=facpub
+
+ xs := {0.0, 1.0, 0.0};
+ ys := {0.0, 0.0, 1.0};
+
+ degree := 2;
+ procedure l(i, j, x, y);
+ binomial(degree, i)*binomial(degree, j)*(
+ det(mat(
+ (x, y, 1),
+ (part(xs, i+1), part(ys, i+1), 1),
+ (part(xs, j+1), part(ys, j+1), 1)
+ ))
+ );
+
+
+ det mat(
+ (l(2, 1, x, y), l(2, 0, x, y)),
+ (l(2, 0, x, y), l(1, 0, x, y))
+ );
+
+ t1 := l(2, 0, x, y) / (l(2, 0, x, y) - l(2, 1, x, y));
+ t2 := l(1, 0, x, y) / (l(1, 0, x, y) - l(2, 0, x, y));
+
+ t1 := sub(x=part(xs, 1)*u + part(xs, 2)*v + part(xs, 3)*(1.0 - u - v), t1);
+ t1 := sub(y=part(ys, 1)*u + part(ys, 2)*v + part(ys, 3)*(1.0 - u - v), t1);
+
+ t2 := sub(x=part(xs, 1)*u + part(xs, 2)*v + part(xs, 3)*(1.0 - u - v), t2);
+ t2 := sub(y=part(ys, 1)*u + part(ys, 2)*v + part(ys, 3)*(1.0 - u - v), t2);
+ */
+
+
+ Float x1 = ps[0][0], y1 = ps[0][1], z1 = ps[0][2];
+
+ Float px = p[0] - x1;
+ Float py = p[1] - y1;
+ Float pz = p[2] - z1;
+
+ Float x2 = ps[1][0] - x1, y2 = ps[1][1] - y1, z2 = ps[1][2] - z1;
+ Float x3 = ps[2][0] - x1, y3 = ps[2][1] - y1, z3 = ps[2][2] - z1;
+
+ Float x22 = x2 * x2, y22 = y2 * y2, z22 = z2 * z2;
+ Float x32 = x3 * x3, y32 = y3 * y3, z32 = z3 * z3;
+
+ Float B =
+ (((2 * (4 * z22 - 2 * z2 * z3 - z32 - y32) + x32 + 4 * (2 * y2 - y3) * y2) * x2 -
+ 2 * (((2 * y2 - 3 * y3) * y2 + (2 * z2 - 3 * z3) * z2) * x3 - 2 * (x2 - x3) * x22)) *
+ x2 -
+ (x32 + y32 + (2 * z2 - z3) * (2 * z2 - z3) + 4 * (y2 - y3) * y2 + 4 * (x2 - x3) * x2) *
+ ((2 * x2 - x3) * px + (2 * y2 - y3) * py) -
+ (x32 + y32 + (2 * z2 - z3) * (2 * z2 - z3) + 4 * (y2 - y3) * y2 + 4 * (x2 - x3) * x2) *
+ (2 * z2 - z3) * pz +
+ 2 * (2 * (y2 - y3) * y22 - (2 * z2 - 3 * z3) * y3 * z2) * y2 +
+ (2 * (4 * z22 - 2 * z2 * z3 - z32) + y32) * y22 +
+ ((2 * z2 - z3) * (2 * z2 - z3) - 2 * y32) * z22 - 2 * (y22 + z22) * x32);
+
+ Float C = (3 * (x32 + y32 + (2 * z2 - z3) * (2 * z2 - z3) + 4 * (y2 - y3) * y2) +
+ 12 * (x2 - x3) * x2);
+ const Float eps = 0.000001;
+
+ if (B < 0.0 && B > -eps) {
+ B = 0.0;
+ }
+
+ if (C == 0.0 || B < 0.0) {
+ Vector p2 = evaluate(0.0);
+ Vector vec = p2 - p;
+
+ r_dis = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
+ r_tan = derivative(0.0);
+ r_s = 0.0;
+
+ return p2;
+ }
+
+ B = sqrt(B);
+
+ Float t1 = (3 * (2 * z2 - z3) * z2 + B * sqrt(3) + 3 * (2 * y2 - y3) * y2 +
+ 3 * (2 * x2 - x3) * x2) /
+ C;
+
+ Float t2 = (3 * (2 * z2 - z3) * z2 - B * sqrt(3) + 3 * (2 * y2 - y3) * y2 +
+ 3 * (2 * x2 - x3) * x2) /
+ C;
+
+ t1 = min_ff(max_ff(t1, 0.0), 1.0);
+ t2 = min_ff(max_ff(t2, 0.0), 1.0);
+
+ // t1 = std::min(std::max(t1, 0.0), 1.0);
+ // t2 = std::min(std::max(t2, 0.0), 1.0);
+
+ Float s1 = t_to_arc(t1);
+ Float s2 = t_to_arc(t2);
+
+ Vector a = evaluate(s1);
+ Vector b = evaluate(s2);
+
+ a -= p;
+ b -= p;
+
+ Float dis1 = _dot(a, a);
+ Float dis2 = _dot(b, b);
+
+ if (dis1 < dis2) {
+ r_s = s1;
+ r_dis = std::sqrt(dis1);
+ r_tan = derivative(s1);
+
+ return a + p;
+ }
+ else {
+ r_s = s2;
+ r_dis = std::sqrt(dis2);
+ r_tan = derivative(s2);
+
+ return b + p;
+ }
+# endif
+ }
+
+ QuadBezier &operator=(QuadBezier &&b)
+ {
+ ps[0] = b.ps[0];
+ ps[1] = b.ps[1];
+ ps[2] = b.ps[2];
+
+ length = b.length;
+
+ if (b._arc_to_t) {
+ _arc_to_t = std::move(b._arc_to_t);
+ _t_to_arc = std::move(b._t_to_arc);
+
+ b._arc_to_t = nullptr;
+ b._t_to_arc = nullptr;
+ }
+ else {
+ _arc_to_t = new Float[table_size];
+ _t_to_arc = new Float[table_size];
+ }
+
+ return *this;
+ }
+#endif
+
+ Float length;
+
+ void update()
+ {
+ Float t = 0.0, dt = 1.0 / (Float)table_size;
+ Float s = 0.0;
+
+ if (!_arc_to_t) {
+ _arc_to_t = new Float[table_size];
+ _t_to_arc = new Float[table_size];
+ }
+
+ for (int i : IndexRange(table_size)) {
+ _arc_to_t[i] = -1.0;
+ }
+
+ length = 0.0;
+
+ for (int i = 0; i < table_size; i++, t += dt) {
+ Float dlen = 0.0;
+ for (int j = 0; j < axes; j++) {
+ float dv = dquad(ps[0][j], ps[1][j], ps[2][j], t);
+
+ dlen += dv * dv;
+ }
+
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list