summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/knot.h
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/knot.h')
-rw-r--r--Build/source/utils/asymptote/knot.h533
1 files changed, 533 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/knot.h b/Build/source/utils/asymptote/knot.h
new file mode 100644
index 00000000000..e087f271e61
--- /dev/null
+++ b/Build/source/utils/asymptote/knot.h
@@ -0,0 +1,533 @@
+/*****
+ * knot.h
+ * Andy Hammerlindl 200/02/10
+ *
+ * Describes a knot, a point and its neighbouring specifiers, used as an
+ * intermediate structure in solving paths.
+ *****/
+
+#ifndef KNOT_H
+#define KNOT_H
+
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+#include "mod.h"
+#include "pair.h"
+#include "path.h"
+
+namespace camp {
+
+using mem::vector;
+
+// The choice of branch cuts of atan2 disguishes between y=+0.0 and y=-0.0 in
+// the case where x<0. This can lead to strange looking paths being
+// calculated from guides of the form a..b..cycle. To avoid these degenerate
+// cases, the niceAngle routine moves the branch cut so that the sign of a
+// zero won't matter.
+double niceAngle(pair z);
+
+// A cyclic vector: ie. a vector where the index is taken mod the size of the
+// vector.
+template <typename T>
+class cvector : public vector<T> {
+public:
+ cvector() {}
+ cvector(size_t n) : vector<T>(n) {}
+ cvector(size_t n, const T& t) : vector<T>(n,t) {}
+ cvector(const vector<T>& v) : vector<T>(v) {}
+
+ T& operator[](Int j) {
+ return vector<T>::operator[](imod(j,(Int) this->size()));
+ }
+ const T& operator[](Int j) const {
+ return vector<T>::operator[](imod(j,(Int) this->size()));
+ }
+};
+
+// Forward declaration.
+class knotlist;
+
+/* A linear equation (one of a set of equations to solve for direction through
+ knots in a path). The i-th equation is:
+
+ pre*theta[i-1] + piv*theta[i] + post*theta[i+1] = aug
+
+ where indices are taken mod n.
+*/
+struct eqn {
+ double pre,piv,post,aug;
+ eqn(double pre, double piv, double post, double aug)
+ : pre(pre), piv(piv), post(post), aug(aug) {}
+
+ friend ostream& operator<< (ostream& out, const eqn& e)
+ {
+ return out << e.pre << " * pre + "
+ << e.piv << " * piv + "
+ << e.post << " * post = "
+ << e.aug;
+ }
+};
+
+
+// A direction specifier, telling how the path behaves coming in or out of a
+// point. The base class represents the "open" specifier.
+class spec : public gc {
+public:
+ virtual ~spec() {}
+
+ // If the knot is open, it gives no restriction on the behavior of the
+ // path.
+ virtual bool open() { return true; }
+ virtual bool controlled() { return false; }
+ virtual pair control() {return pair(0.0,0.0);}
+ virtual double curl() { return -1.0; }
+ virtual pair dir() { return pair(0.0,0.0); }
+
+ // When a knot has a restriction on one side but is open on the other, the
+ // restriction implies a restriction on the other side. This is the partner
+ // restriction defined here, where the pair argument is for the location of
+ // the knot.
+ virtual spec *outPartner(pair) { return this; }
+ virtual spec *inPartner(pair) { return this; }
+
+ virtual void print(ostream&) const {}
+};
+
+inline ostream& operator<< (ostream& out, spec& s)
+{
+ s.print(out);
+ return out;
+}
+
+
+// Specifier used at an endpoint.
+class endSpec : public spec {
+public:
+ bool open() { return false; }
+
+ // Returns an equation used to solve for the thetas along the knot. These are
+ // called by eqnprop in the non-cyclic case for the first and last equations.
+ virtual eqn eqnOut(Int j, knotlist& l,
+ cvector<double>& d, cvector<double>& psi) = 0;
+ virtual eqn eqnIn (Int j, knotlist& l,
+ cvector<double>& d, cvector<double>& psi) = 0;
+};
+
+// A specifier with a given direction (in radians).
+class dirSpec : public endSpec {
+ double given;
+public:
+ // Direction should be given in the range [-PI,PI]
+ dirSpec(double given)
+ : given(given) {}
+ dirSpec(pair z)
+ : given(niceAngle(z)) {}
+
+ pair dir() { return expi(given); }
+
+ eqn eqnOut(Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+ eqn eqnIn (Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+
+ void print(ostream& out) const {
+ out << "{dir(" << degrees(given) << ")}";
+ }
+};
+
+// A curl specifier. The curvature at the end knot should be gamma times the
+// curvature at the neighbouring knot.
+class curlSpec : public endSpec {
+ double gamma;
+public:
+ // Gamma should be non-negative.
+ curlSpec(double gamma=1.0)
+ : gamma(gamma) {
+ if(gamma < 0)
+ reportError("curl cannot be less than 0");
+ }
+
+ double curl() { return gamma; }
+
+ eqn eqnOut(Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+ eqn eqnIn (Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+
+ void print(ostream& out) const {
+ out << "{curl " << gamma << "}";
+ }
+};
+
+
+// A specifier with a control point. All information for this portion of the
+// curve has been determined.
+class controlSpec : public spec {
+public:
+ pair cz;
+ bool straight;
+
+ controlSpec(pair cz, bool straight=false)
+ : cz(cz), straight(straight) {}
+
+ bool open() { return false; }
+ bool controlled() { return true; }
+ pair control() { return cz; }
+
+ // The partner spec will be a dirSpec in the same direction the specifier
+ // takes the path, unless the velocity is zero, then it uses a curl
+ // specifier.
+ spec *outPartner(pair);
+ spec *inPartner(pair);
+
+ void print(ostream& out) const {
+ // NOTE: This format cannot be read back in.
+ out << "{control " << cz << "}";
+ }
+};
+
+
+// The tension information for one side of a knot.
+struct tension {
+ double val;
+ bool atleast;
+
+ tension(double val=1.0, bool atleast=false)
+ : val(val), atleast(atleast) {
+ if(val < 0.75)
+ reportError("tension cannot be less than 3/4");
+ }
+};
+
+inline
+ostream& operator<<(ostream& out, tension t)
+{
+ return out << "tension" << (t.atleast ? " atleast " : " ") << t.val;
+}
+
+// A knot, a point with specifiers to double the path coming in and going out
+// of the knot.
+struct knot {
+ pair z;
+ spec *in;
+ spec *out;
+ tension tin, tout;
+
+ knot() {}
+ knot(pair z, spec *in, spec *out,
+ tension tin=tension(), tension tout=tension())
+ : z(z), in(in), out(out), tin(tin), tout(tout) {}
+
+ double alpha() { return 1.0/tout.val; }
+ double beta() { return 1.0/tin.val; }
+};
+
+ostream& operator<<(ostream& out, const knot& k);
+
+// Abstract base class for a section of a guide.
+class knotlist {
+public:
+ virtual ~knotlist() {}
+
+ virtual Int length() = 0;
+ virtual bool cyclic() = 0;
+
+ // Returns the number of knots.
+ Int size() {
+ return cyclic() ? length() : length() + 1;
+ }
+ bool empty() {
+ return size()==0;
+ }
+
+ virtual knot& cell(Int) = 0;
+ virtual knot& operator[] (Int i) {
+#if 0
+ assert(cyclic() || (0 <= i && i <= length())); // Bounds check.
+#endif
+ return cell(i);
+ }
+ knot& front() {
+ return (*this)[0];
+ }
+ knot& back() {
+ return (*this)[length()];
+ }
+};
+
+
+// Defines a knotlist as a piece of another knotlist.
+class subknotlist : public knotlist {
+ knotlist& l;
+ Int a,b;
+public:
+ subknotlist(knotlist& l, Int a, Int b)
+ : l(l), a(a), b(b) {}
+
+ Int length() { return b-a; }
+ bool cyclic() { return false; }
+ knot& cell(Int i) { return l[a+i]; }
+};
+
+struct simpleknotlist : public knotlist {
+ cvector<knot> nodes;
+ bool cycles;
+
+ simpleknotlist(cvector<knot> nodes, bool cycles=false)
+ : nodes(nodes), cycles(cycles) {}
+
+ Int length() { return cycles ? (Int) nodes.size() : (Int) nodes.size() - 1; }
+ bool cyclic() { return cycles; }
+ knot& cell(Int j) { return nodes[j]; }
+};
+
+// A protopath is a path being made.
+struct protopath {
+ bool cycles;
+ Int n;
+ mem::vector<solvedKnot> nodes;
+
+ protopath(Int n, bool cycles)
+ : cycles(cycles), n(n), nodes(n) {}
+
+ solvedKnot& operator[](Int j) {
+ return nodes[imod(j,n)];
+ }
+
+ bool& straight(Int j) {
+ return (*this)[j].straight;
+ }
+ pair& pre(Int j) {
+ return (*this)[j].pre;
+ }
+ pair& point(Int j) {
+ return (*this)[j].point;
+ }
+ pair& post(Int j) {
+ return (*this)[j].post;
+ }
+
+ void controlEnds() {
+ if (!cycles) {
+ solvedKnot& start=(*this)[0];
+ solvedKnot& end=(*this)[n-1];
+
+ start.pre=start.point;
+ end.post=end.point;
+ }
+ }
+ // Once all the controls are set, return the final (constant) path.
+ path fix() {
+ return path(nodes,n,cycles);
+ }
+};
+
+
+// Represents properties that can be computed along a knotlist.
+// Examples include distances (d), turning angles (psi), and the linear
+// equations used to solve for the thetas.
+template <typename T>
+class knotprop {
+protected:
+ knotlist& l;
+
+ // Calculate the property for the usual case in the iteration (and for a
+ // cyclic knot, the only case), at the index given.
+ virtual T mid(Int) = 0;
+
+ // The special cases, these default to the usual case: mid.
+ virtual T solo(Int j) // Calculates the property for a list of length 0.
+ {
+ return mid(j);
+ }
+ virtual T start(Int j) // Calculates it at the start of the list.
+ {
+ return mid(j);
+ }
+ virtual T end(Int j) // Calculate it at the end.
+ {
+ return mid(j);
+ }
+
+ virtual cvector<T> linearCompute()
+ {
+ Int n=l.length();
+ cvector<T> v;
+ if (n==0)
+ v.push_back(solo(0));
+ else {
+ v.push_back(start(0));
+ for (Int j=1; j<n; ++j)
+ v.push_back(mid(j));
+ v.push_back(end(n));
+ }
+ return v;
+ }
+
+ virtual cvector<T> cyclicCompute()
+ {
+ Int n=l.length();
+ cvector<T> v;
+ for (Int j=0; j<n; ++j)
+ v.push_back(mid(j));
+ return v;
+ }
+
+ virtual cvector<T> linearBackCompute()
+ {
+ Int n=l.length();
+ cvector<T> v;
+ if (n==0)
+ v.push_back(solo(0));
+ else {
+ v.push_back(end(n));
+ for (Int j=1; j<n; ++j)
+ v.push_back(mid(n-j));
+ v.push_back(start(0));
+ }
+ return v;
+ }
+
+ virtual cvector<T> cyclicBackCompute()
+ {
+ Int n=l.length();
+ cvector<T> v;
+ for (Int j=1; j<=n; ++j)
+ v.push_back(mid(n-j));
+ return v;
+ }
+
+public:
+ virtual ~knotprop() {}
+
+ virtual cvector<T> compute() {
+ return l.cyclic() ? cyclicCompute() : linearCompute();
+ }
+
+ // Compute the values in the opposite order. This is needed for instance if
+ // the i-th calculation needed a result computed in the i+1-th, such as in the
+ // back substitution for solving thetas.
+ virtual cvector<T> backCompute() {
+ cvector<T> v=l.cyclic() ? cyclicBackCompute() : linearBackCompute();
+
+ // Even though they are computed in the backwards order, return them in the
+ // standard order.
+ reverse(v.begin(),v.end());
+ return v;
+ }
+
+ knotprop(knotlist& l)
+ : l(l) {}
+};
+
+// A knot transforms, it takes in one knotlist and transforms it knot for knot
+// into a new one.
+class knottrans : public knotprop<knot> {
+protected:
+ virtual knot mid(Int j) {
+ /* By default, just copy the knot. */
+ return l[j];
+ }
+
+public:
+ virtual ~knottrans() {}
+
+ knottrans(knotlist& l)
+ : knotprop<knot>(l) {}
+
+ virtual simpleknotlist trans() {
+ return simpleknotlist(compute(),l.cyclic());
+ }
+};
+
+// Like a knotprop, but it doesn't compute a vector of values for the knot. It
+// iterates over the knotlist calling method for side-effect. For instance,
+// this is used to plug control points into protopaths.
+class knoteffect {
+protected:
+ knotlist& l;
+
+ virtual void mid(Int) = 0;
+
+ // The special cases, these default to the usual case: mid.
+ virtual void solo(Int j) {
+ mid(j);
+ }
+ virtual void start(Int j) {
+ mid(j);
+ }
+ virtual void end(Int j) {
+ mid(j);
+ }
+
+ virtual void linearExec()
+ {
+ Int n=l.length();
+ if (n==0)
+ solo(0);
+ else {
+ start(0);
+ for (Int j=1; j<n; ++j)
+ mid(j);
+ end(n);
+ }
+ }
+
+ virtual void cyclicExec()
+ {
+ Int n=l.length();
+ for (Int j=0; j<n; ++j)
+ mid(j);
+ }
+
+ virtual void linearBackExec()
+ {
+ Int n=l.length();
+ if (n==0)
+ solo(0);
+ else {
+ end(n);
+ for (Int j=1; j<n; ++j)
+ mid(n-j);
+ start(0);
+ }
+ }
+
+ virtual void cyclicBackExec()
+ {
+ Int n=l.length();
+ for (Int j=1; j<=n; ++j)
+ mid(n-j);
+ }
+
+public:
+ virtual ~knoteffect() {}
+
+ virtual void exec() {
+ if (l.cyclic())
+ cyclicExec();
+ else
+ linearExec();
+ }
+
+ virtual void backCompute() {
+ if (l.cyclic())
+ cyclicBackExec();
+ else
+ linearBackExec();
+ }
+
+ knoteffect(knotlist& l)
+ : l(l) {}
+};
+
+path solve(knotlist& l);
+
+path solveSimple(cvector<pair>& z);
+
+double velocity(double theta, double phi, tension t);
+
+} // namespace camp
+
+GC_DECLARE_PTRFREE(camp::eqn);
+GC_DECLARE_PTRFREE(camp::tension);
+
+#endif // KNOT_H