Skip to content

Commit

Permalink
Add interpolating splines: both periodic splines (that form a
Browse files Browse the repository at this point in the history 
loop), and open-ended splines, with their tangents specified at
their endpoints.

Also change constraint solver matrix size to 1024, on the theory
that a power of two will generate better array indexing, and
replace fabs() with my own function that for some reason is
faster.

[git-p4: depot-paths = "//depot/solvespace/": change = 2055]
  • Loading branch information
@jwesthues
jwesthues committed Oct 21, 2009
1 parent f6bb680 commit 2ca5334
Show file tree
Hide file tree
Showing 18 changed files with 388 additions and 55 deletions.
8 changes: 4 additions & 4 deletions constraint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -591,12 +591,12 @@ void Constraint::MenuConstrain(int id) {
}
Vector l0 = SK.GetEntity(line->point[0])->PointGetNum(),
l1 = SK.GetEntity(line->point[1])->PointGetNum();
Vector a0 = SK.GetEntity(cubic->point[0])->PointGetNum(),
a3 = SK.GetEntity(cubic->point[3])->PointGetNum();
Vector as = cubic->CubicGetStartNum(),
af = cubic->CubicGetFinishNum();

if(l0.Equals(a0) || l1.Equals(a0)) {
if(l0.Equals(as) || l1.Equals(as)) {
c.other = false;
} else if(l0.Equals(a3) || l1.Equals(a3)) {
} else if(l0.Equals(af) || l1.Equals(af)) {
c.other = true;
} else {
Error("The tangent cubic and line segment must share an "
Expand Down
12 changes: 6 additions & 6 deletions constrainteq.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -639,12 +639,12 @@ void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) {
EntityBase *cubic = SK.GetEntity(entityA);
EntityBase *line = SK.GetEntity(entityB);

ExprVector endpoint =
SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetExprs();
ExprVector ctrlpoint =
SK.GetEntity(cubic->point[other ? 2 : 1])->PointGetExprs();

ExprVector a = endpoint.Minus(ctrlpoint);
ExprVector a;
if(other) {
a = cubic->CubicGetFinishTangentExprs();
} else {
a = cubic->CubicGetStartTangentExprs();
}

ExprVector b = line->VectorGetExprs();

Expand Down
13 changes: 12 additions & 1 deletion draw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,7 @@ void GraphicsWindow::GroupSelection(void) {
case Entity::WORKPLANE: (gs.workplanes)++; break;
case Entity::LINE_SEGMENT: (gs.lineSegments)++; break;
case Entity::CUBIC: (gs.cubics)++; break;
// but don't count periodic cubics

case Entity::ARC_OF_CIRCLE:
(gs.circlesOrArcs)++;
Expand Down Expand Up @@ -203,7 +204,17 @@ void GraphicsWindow::HitTestMakeSelection(Point2d mp) {
for(i = 0; i < SK.entity.n; i++) {
Entity *e = &(SK.entity.elem[i]);
// Don't hover whatever's being dragged.
if(e->h.request().v == pending.point.request().v) continue;
if(e->h.request().v == pending.point.request().v) {
// The one exception is when we're creating a new cubic; we
// want to be able to hover the first point, because that's
// how we turn it into a periodic spline.
if(!e->IsPoint()) continue;
if(!e->h.isFromRequest()) continue;
Request *r = SK.GetRequest(e->h.request());
if(r->type != Request::CUBIC) continue;
if(r->extraPoints < 2) continue;
if(e->h.v != r->h.entity(1).v) continue;
}

d = e->GetDistance(mp);
if(d < 10 && d < dmin) {
Expand Down
4 changes: 2 additions & 2 deletions drawconstraint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -699,8 +699,8 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) {
}

Entity *cubic = SK.GetEntity(entityA);
Vector p =
SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetNum();
Vector p = other ? cubic->CubicGetFinishNum() :
cubic->CubicGetStartNum();
Vector dir = SK.GetEntity(entityB)->VectorGetNum();
Vector out = n.Cross(dir);
textAt = p.Plus(out.WithMagnitude(14/SS.GW.scale));
Expand Down
152 changes: 144 additions & 8 deletions drawentity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,6 +215,143 @@ bool Entity::PointIsFromReferences(void) {
return h.request().IsFromReferences();
}

//-----------------------------------------------------------------------------
// Compute a cubic, second derivative continuous, interpolating spline. Same
// routine for periodic splines (in a loop) or open splines (with specified
// end tangents).
//-----------------------------------------------------------------------------
void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) {
static const int MAX_N = BandedMatrix::MAX_UNKNOWNS;
int ep = extraPoints;

// The number of unknowns to solve for.
int n = periodic ? 3 + ep : ep;
if(n >= MAX_N) oops();
// The number of on-curve points, one more than the number of segments.
int pts = periodic ? 4 + ep : 2 + ep;

int i, j, a;

// The starting and finishing control points that define our end tangents
// (if the spline isn't periodic), and the on-curve points.
Vector ctrl_s, ctrl_f, pt[MAX_N+4];
if(periodic) {
for(i = 0; i < ep + 3; i++) {
pt[i] = SK.GetEntity(point[i])->PointGetNum();
}
pt[i++] = SK.GetEntity(point[0])->PointGetNum();
} else {
ctrl_s = SK.GetEntity(point[1])->PointGetNum();
ctrl_f = SK.GetEntity(point[ep+2])->PointGetNum();
j = 0;
pt[j++] = SK.GetEntity(point[0])->PointGetNum();
for(i = 2; i <= ep + 1; i++) {
pt[j++] = SK.GetEntity(point[i])->PointGetNum();
}
pt[j++] = SK.GetEntity(point[ep+3])->PointGetNum();
}

// The unknowns that we will be solving for, a set for each coordinate.
double Xx[MAX_N], Xy[MAX_N], Xz[MAX_N];
// For a cubic Bezier section f(t) as t goes from 0 to 1,
// f' (0) = 3*(P1 - P0)
// f' (1) = 3*(P3 - P2)
// f''(0) = 6*(P0 - 2*P1 + P2)
// f''(1) = 6*(P3 - 2*P2 + P1)
for(a = 0; a < 3; a++) {
BandedMatrix bm;
ZERO(&bm);
bm.n = n;

for(i = 0; i < n; i++) {
int im, it, ip;
if(periodic) {
im = WRAP(i - 1, n);
it = i;
ip = WRAP(i + 1, n);
} else {
im = i;
it = i + 1;
ip = i + 2;
}
// All of these are expressed in terms of a constant part, and
// of X[i-1], X[i], and X[i+1]; so let these be the four
// components of that vector;
Vector4 A, B, C, D, E;
// The on-curve interpolated point
C = Vector4::From((pt[it]).Element(a), 0, 0, 0);
// control point one back, C - X[i]
B = C.Plus(Vector4::From(0, 0, -1, 0));
// control point one forward, C + X[i]
D = C.Plus(Vector4::From(0, 0, 1, 0));
// control point two back
if(i == 0 && !periodic) {
A = Vector4::From(ctrl_s.Element(a), 0, 0, 0);
} else {
// pt[im] + X[i-1]
A = Vector4::From(pt[im].Element(a), 1, 0, 0);
}
// control point two forward
if(i == (n - 1) && !periodic) {
E = Vector4::From(ctrl_f.Element(a), 0, 0, 0);
} else {
// pt[ip] - X[i+1]
E = Vector4::From((pt[ip]).Element(a), 0, 0, -1);
}
// Write the second derivatives of each segment, dropping constant
Vector4 fprev_pp = (C.Minus(B.ScaledBy(2))).Plus(A),
fnext_pp = (C.Minus(D.ScaledBy(2))).Plus(E),
eq = fprev_pp.Minus(fnext_pp);

bm.B[i] = -eq.w;
if(periodic) {
bm.A[i][WRAP(i-2, n)] = eq.x;
bm.A[i][WRAP(i-1, n)] = eq.y;
bm.A[i][i] = eq.z;
} else {
// The wrapping would work, except when n = 1 and everything
// wraps to zero...
if(i > 0) bm.A[i][i - 1] = eq.x;
bm.A[i][i] = eq.y;
if(i < (n-1)) bm.A[i][i + 1] = eq.z;
}
}
bm.Solve();
double *X = (a == 0) ? Xx :
(a == 1) ? Xy :
Xz;
memcpy(X, bm.X, n*sizeof(double));
}

for(i = 0; i < pts - 1; i++) {
Vector p0, p1, p2, p3;
if(periodic) {
p0 = pt[i];
int iw = WRAP(i - 1, n);
p1 = p0.Plus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
} else if(i == 0) {
p0 = pt[0];
p1 = ctrl_s;
} else {
p0 = pt[i];
p1 = p0.Plus(Vector::From(Xx[i-1], Xy[i-1], Xz[i-1]));
}
if(periodic) {
p3 = pt[i+1];
int iw = WRAP(i, n);
p2 = p3.Minus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
} else if(i == (pts - 2)) {
p3 = pt[pts-1];
p2 = ctrl_f;
} else {
p3 = pt[i+1];
p2 = p3.Minus(Vector::From(Xx[i], Xy[i], Xz[i]));
}
SBezier sb = SBezier::From(p0, p1, p2, p3);
sbl->l.Add(&sb);
}
}

void Entity::GenerateBezierCurves(SBezierList *sbl) {
SBezier sb;

Expand All @@ -228,15 +365,13 @@ void Entity::GenerateBezierCurves(SBezierList *sbl) {
sbl->l.Add(&sb);
break;
}
case CUBIC: {
Vector p0 = SK.GetEntity(point[0])->PointGetNum();
Vector p1 = SK.GetEntity(point[1])->PointGetNum();
Vector p2 = SK.GetEntity(point[2])->PointGetNum();
Vector p3 = SK.GetEntity(point[3])->PointGetNum();
sb = SBezier::From(p0, p1, p2, p3);
sbl->l.Add(&sb);
case CUBIC:
ComputeInterpolatingSpline(sbl, false);
break;

case CUBIC_PERIODIC:
ComputeInterpolatingSpline(sbl, true);
break;
}

case CIRCLE:
case ARC_OF_CIRCLE: {
Expand Down Expand Up @@ -476,6 +611,7 @@ void Entity::DrawOrGetDistance(void) {
case CIRCLE:
case ARC_OF_CIRCLE:
case CUBIC:
case CUBIC_PERIODIC:
case TTF_TEXT:
// Nothing but the curve(s).
break;
Expand Down
14 changes: 14 additions & 0 deletions dsc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -367,4 +367,18 @@ class NameStr {
}
};

class BandedMatrix {
public:
static const int MAX_UNKNOWNS = 16;
static const int RIGHT_OF_DIAG = 1;
static const int LEFT_OF_DIAG = 2;

double A[MAX_UNKNOWNS][MAX_UNKNOWNS];
double B[MAX_UNKNOWNS];
double X[MAX_UNKNOWNS];
int n;

void Solve(void);
};

#endif
17 changes: 17 additions & 0 deletions entity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,23 @@ void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) {
while(*dtheta > (2*PI)) *dtheta -= 2*PI;
}

Vector EntityBase::CubicGetStartNum(void) {
return SK.GetEntity(point[0])->PointGetNum();
}
Vector EntityBase::CubicGetFinishNum(void) {
return SK.GetEntity(point[3+extraPoints])->PointGetNum();
}
ExprVector EntityBase::CubicGetStartTangentExprs(void) {
ExprVector pon = SK.GetEntity(point[0])->PointGetExprs(),
poff = SK.GetEntity(point[1])->PointGetExprs();
return (pon.Minus(poff));
}
ExprVector EntityBase::CubicGetFinishTangentExprs(void) {
ExprVector pon = SK.GetEntity(point[3+extraPoints])->PointGetExprs(),
poff = SK.GetEntity(point[2+extraPoints])->PointGetExprs();
return (pon.Minus(poff));
}

bool EntityBase::IsWorkplane(void) {
return (type == WORKPLANE);
}
Expand Down
10 changes: 10 additions & 0 deletions file.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {

{ 'r', "Request.h.v", 'x', &(SS.sv.r.h.v) },
{ 'r', "Request.type", 'd', &(SS.sv.r.type) },
{ 'r', "Request.extraPoints", 'd', &(SS.sv.r.extraPoints) },
{ 'r', "Request.workplane.v", 'x', &(SS.sv.r.workplane.v) },
{ 'r', "Request.group.v", 'x', &(SS.sv.r.group.v) },
{ 'r', "Request.construction", 'b', &(SS.sv.r.construction) },
Expand All @@ -120,6 +121,15 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {
{ 'e', "Entity.point[1].v", 'x', &(SS.sv.e.point[1].v) },
{ 'e', "Entity.point[2].v", 'x', &(SS.sv.e.point[2].v) },
{ 'e', "Entity.point[3].v", 'x', &(SS.sv.e.point[3].v) },
{ 'e', "Entity.point[4].v", 'x', &(SS.sv.e.point[4].v) },
{ 'e', "Entity.point[5].v", 'x', &(SS.sv.e.point[5].v) },
{ 'e', "Entity.point[6].v", 'x', &(SS.sv.e.point[6].v) },
{ 'e', "Entity.point[7].v", 'x', &(SS.sv.e.point[7].v) },
{ 'e', "Entity.point[8].v", 'x', &(SS.sv.e.point[8].v) },
{ 'e', "Entity.point[9].v", 'x', &(SS.sv.e.point[9].v) },
{ 'e', "Entity.point[10].v", 'x', &(SS.sv.e.point[10].v) },
{ 'e', "Entity.point[11].v", 'x', &(SS.sv.e.point[11].v) },
{ 'e', "Entity.extraPoints", 'd', &(SS.sv.e.extraPoints) },
{ 'e', "Entity.normal.v", 'x', &(SS.sv.e.normal.v) },
{ 'e', "Entity.distance.v", 'x', &(SS.sv.e.distance.v) },
{ 'e', "Entity.workplane.v", 'x', &(SS.sv.e.workplane.v) },
Expand Down
19 changes: 14 additions & 5 deletions group.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -636,6 +636,7 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
Entity en;
memset(&en, 0, sizeof(en));
en.type = ep->type;
en.extraPoints = ep->extraPoints;
en.h = Remap(ep->h, remap);
en.timesApplied = timesApplied;
en.group = h;
Expand All @@ -652,12 +653,20 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
en.point[1] = Remap(ep->point[1], remap);
break;

case Entity::CUBIC:
en.point[0] = Remap(ep->point[0], remap);
en.point[1] = Remap(ep->point[1], remap);
en.point[2] = Remap(ep->point[2], remap);
en.point[3] = Remap(ep->point[3], remap);
case Entity::CUBIC: {
int i;
for(i = 0; i < 4 + ep->extraPoints; i++) {
en.point[i] = Remap(ep->point[i], remap);
}
break;
}
case Entity::CUBIC_PERIODIC: {
int i;
for(i = 0; i < 3 + ep->extraPoints; i++) {
en.point[i] = Remap(ep->point[i], remap);
}
break;
}

case Entity::CIRCLE:
en.point[0] = Remap(ep->point[0], remap);
Expand Down
4 changes: 2 additions & 2 deletions modify.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,7 @@ hEntity GraphicsWindow::SplitCubic(hEntity he, Vector pinter) {
b01.ClosestPointTo(pinter, &t, true);
b01.SplitAt(t, &b0i, &bi1);

// Add the two line segments this one gets split into.
// Add the two cubic segments this one gets split into.
hRequest r0i = AddRequest(Request::CUBIC, false),
ri1 = AddRequest(Request::CUBIC, false);
// Don't get entities till after adding, realloc issues
Expand Down Expand Up @@ -343,7 +343,7 @@ hEntity GraphicsWindow::SplitEntity(hEntity he, Vector pinter) {
ret = SplitCircle(he, pinter);
} else if(e->type == Entity::LINE_SEGMENT) {
ret = SplitLine(he, pinter);
} else if(e->type == Entity::CUBIC) {
} else if(e->type == Entity::CUBIC && e->extraPoints == 0) {
ret = SplitCubic(he, pinter);
} else {
Error("Couldn't split this entity; lines, circles, or cubics only.");
Expand Down
Loading

0 comments on commit 2ca5334

Please sign in to comment.