summaryrefslogtreecommitdiff
path: root/graphics/asymptote/triple.h
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /graphics/asymptote/triple.h
Initial commit
Diffstat (limited to 'graphics/asymptote/triple.h')
-rw-r--r--graphics/asymptote/triple.h365
1 files changed, 365 insertions, 0 deletions
diff --git a/graphics/asymptote/triple.h b/graphics/asymptote/triple.h
new file mode 100644
index 0000000000..de6767a29f
--- /dev/null
+++ b/graphics/asymptote/triple.h
@@ -0,0 +1,365 @@
+/*****
+ * triple.h
+ * John Bowman
+ *
+ * Stores a three-dimensional point.
+ *
+ *****/
+
+#ifndef TRIPLE_H
+#define TRIPLE_H
+
+#include <cassert>
+#include <iostream>
+#include <cmath>
+#include <cstring>
+
+#include "common.h"
+#include "angle.h"
+#include "pair.h"
+
+namespace run {
+void transpose(double *a, size_t n);
+void inverse(double *a, size_t n);
+}
+
+namespace camp {
+
+typedef double Triple[3];
+
+class triple;
+
+bool isIdTransform3(const double* t);
+void copyTransform3(double*& d, const double* s, GCPlacement placement=NoGC);
+void multiplyTransform3(double*& t, const double* s, const double* r);
+
+void boundstriples(double& x, double& y, double& z, double& X, double& Y,
+ double& Z, size_t n, const triple* v);
+
+class triple : virtual public gc {
+ double x;
+ double y;
+ double z;
+
+public:
+ triple() : x(0.0), y(0.0), z(0.0) {}
+ triple(double x, double y=0.0, double z=0.0) : x(x), y(y), z(z) {}
+ triple(const Triple& v) : x(v[0]), y(v[1]), z(v[2]) {}
+
+ virtual ~triple() {}
+
+ void set(double X, double Y=0.0, double Z=0.0) { x=X; y=Y; z=Z; }
+
+ double getx() const { return x; }
+ double gety() const { return y; }
+ double getz() const { return z; }
+
+ // transform by row-major matrix
+ friend triple operator* (const double* t, const triple& v)
+ {
+ if(t == NULL)
+ return v;
+
+ double f=t[12]*v.x+t[13]*v.y+t[14]*v.z+t[15];
+ if(f == 0.0)
+ reportError("division by 0 in transform of a triple");
+
+ f=1.0/f;
+
+ return triple((t[0]*v.x+t[1]*v.y+t[2]*v.z+t[3])*f,
+ (t[4]*v.x+t[5]*v.y+t[6]*v.z+t[7])*f,
+ (t[8]*v.x+t[9]*v.y+t[10]*v.z+t[11])*f);
+ }
+
+ friend triple transformNormal(const double* t, const triple& v)
+ {
+ if(t == NULL)
+ return v;
+
+ double T[16];
+ memcpy(T,t,sizeof(double)*16);
+ T[3]=T[7]=T[11]=0.0;
+ run::inverse(T,4);
+ run::transpose(T,4);
+ triple V=T*v;
+ return unit(V);
+ }
+
+ friend void transformtriples(const double* t, size_t n, triple* d,
+ const triple* s)
+ {
+ if(n == 0 || d == NULL || s == NULL)
+ return;
+
+ for(size_t i=0; i < n; i++)
+ d[i]=t*s[i];
+ }
+
+ friend void copytriples(size_t n, triple* d, const triple* s)
+ {
+ if(d == NULL || s == NULL)
+ return;
+
+ for(size_t i=0; i < n; i++) d[i]=s[i];
+ }
+
+ friend void boundstriples(triple& Min, triple& Max, size_t n, const triple* v)
+ {
+ if(n==0 || v==NULL)
+ return;
+
+ double x,y,z;
+ double X,Y,Z;
+
+ X=x=v[0].getx();
+ Y=y=v[0].gety();
+ Z=z=v[0].getz();
+ for(size_t i=1; i < n; ++i) {
+ const double vx=v[i].getx();
+ x=fmin(x,vx); X=fmax(X,vx);
+ const double vy=v[i].gety();
+ y=fmin(y,vy); Y=fmax(Y,vy);
+ const double vz=v[i].getz();
+ z=fmin(z,vz); Z=fmax(Z,vz);
+ }
+
+ Min.set(x,y,z);
+ Max.set(X,Y,Z);
+ }
+
+ friend void ratiotriples(pair &b, double (*m)(double, double), bool &first,
+ size_t n, const triple* v)
+ {
+ if(n==0 || v==NULL)
+ return;
+
+ if(first) {
+ first=false;
+ const triple& v0=v[0];
+ b=pair(v0.x/v0.z,v0.y/v0.z);
+ }
+
+ double x=b.getx();
+ double y=b.gety();
+ for(size_t i=0; i < n; ++i) {
+ const triple& vi = v[i];
+ x=m(x,vi.x/vi.z);
+ y=m(y,vi.y/vi.z);
+ }
+ b=pair(x,y);
+ }
+
+ friend triple operator+ (const triple& z, const triple& w)
+ {
+ return triple(z.x + w.x, z.y + w.y, z.z + w.z);
+ }
+
+ friend triple operator- (const triple& z, const triple& w)
+ {
+ return triple(z.x - w.x, z.y - w.y, z.z - w.z);
+ }
+
+ friend triple operator- (const triple& z)
+ {
+ return triple(-z.x, -z.y, -z.z);
+ }
+
+ friend triple operator* (double s, const triple& z)
+ {
+ return triple(s*z.x, s*z.y, s*z.z);
+ }
+
+ friend triple operator* (const triple& z, double s)
+ {
+ return triple(z.x*s, z.y*s, z.z*s);
+ }
+
+ friend triple operator/ (const triple& z, double s)
+ {
+ if (s == 0.0)
+ reportError("division by 0");
+ s=1.0/s;
+ return triple(z.x*s, z.y*s, z.z*s);
+ }
+
+ const triple& operator+= (const triple& w)
+ {
+ x += w.x;
+ y += w.y;
+ z += w.z;
+ return *this;
+ }
+
+ const triple& operator-= (const triple& w)
+ {
+ x -= w.x;
+ y -= w.y;
+ z -= w.z;
+ return *this;
+ }
+
+ friend bool operator== (const triple& z, const triple& w)
+ {
+ return z.x == w.x && z.y == w.y && z.z == w.z;
+ }
+
+ friend bool operator!= (const triple& z, const triple& w)
+ {
+ return z.x != w.x || z.y != w.y || z.z != w.z;
+ }
+
+ double abs2() const
+ {
+ return x*x+y*y+z*z;
+ }
+
+ friend double abs2(const triple &v)
+ {
+ return v.abs2();
+ }
+
+ double length() const /* r */
+ {
+ return sqrt(abs2());
+ }
+
+ friend double length(const triple& v)
+ {
+ return v.length();
+ }
+
+ double polar() const /* theta */
+ {
+ double r=length();
+ if (r == 0.0)
+ reportError("taking polar angle of (0,0,0)");
+ return acos(z/r);
+ }
+
+ double azimuth() const /* phi */
+ {
+ return angle(x,y);
+ }
+
+ friend triple unit(const triple& v)
+ {
+ double scale=v.length();
+ if(scale == 0.0) return v;
+ scale=1.0/scale;
+ return triple(v.x*scale,v.y*scale,v.z*scale);
+ }
+
+ friend double dot(const triple& u, const triple& v)
+ {
+ return u.x*v.x+u.y*v.y+u.z*v.z;
+ }
+
+ friend triple cross(const triple& u, const triple& v)
+ {
+ return triple(u.y*v.z-u.z*v.y,
+ u.z*v.x-u.x*v.z,
+ u.x*v.y-u.y*v.x);
+ }
+
+ // Returns a unit triple in the direction (theta,phi), in radians.
+ friend triple expi(double theta, double phi)
+ {
+ double sintheta=sin(theta);
+ return triple(sintheta*cos(phi),sintheta*sin(phi),cos(theta));
+ }
+
+ friend istream& operator >> (istream& s, triple& z)
+ {
+ char c;
+ s >> std::ws;
+ bool paren=s.peek() == '('; // parenthesis are optional
+ if(paren) s >> c;
+ s >> z.x >> std::ws;
+ if(s.peek() == ',') s >> c >> z.y;
+ else z.y=0.0;
+ if(s.peek() == ',') s >> c >> z.z;
+ else z.z=0.0;
+ if(paren) {
+ s >> std::ws;
+ if(s.peek() == ')') s >> c;
+ }
+
+ return s;
+ }
+
+ friend ostream& operator << (ostream& out, const triple& z)
+ {
+ out << "(" << z.x << "," << z.y << "," << z.z << ")";
+ return out;
+ }
+
+};
+
+triple expi(double theta, double phi);
+
+// Return the component of vector v perpendicular to a unit vector u.
+inline triple perp(triple v, triple u)
+{
+ return v-dot(v,u)*u;
+}
+
+double xratio(const triple& v);
+double yratio(const triple& v);
+
+inline void bounds(double& x, double &X, double v)
+{
+ if(v < x) x=v;
+ else if(v > X) X=v;
+}
+
+inline void boundstriples(double& x, double& y, double& z,
+ double& X, double& Y, double& Z,
+ size_t n, const triple* v)
+{
+ X=x=v[0].getx();
+ Y=y=v[0].gety();
+ Z=z=v[0].getz();
+
+ for(size_t i=1; i < n; ++i) {
+ triple V=v[i];
+ bounds(x,X,V.getx());
+ bounds(y,Y,V.gety());
+ bounds(z,Z,V.getz());
+ }
+}
+
+extern const double third;
+
+// return the maximum distance squared of points c0 and c1 from
+// the respective internal control points of z0--z1.
+inline double Straightness(const triple& z0, const triple& c0,
+ const triple& c1, const triple& z1)
+{
+ triple v=third*(z1-z0);
+ return std::max(abs2(c0-v-z0),abs2(z1-v-c1));
+}
+
+// return the maximum perpendicular distance squared of points c0 and c1
+// from z0--z1.
+inline double Distance1(const triple& z0, const triple& c0,
+ const triple& c1, const triple& z1)
+{
+ triple Z0=c0-z0;
+ triple Q=unit(z1-z0);
+ triple Z1=c1-z0;
+ return std::max(abs2(Z0-dot(Z0,Q)*Q),abs2(Z1-dot(Z1,Q)*Q));
+}
+
+// return the perpendicular distance squared of a point z from the plane
+// through u with unit normal n.
+inline double Distance2(const triple& z, const triple& u, const triple& n)
+{
+ double d=dot(z-u,n);
+ return d*d;
+}
+
+} //namespace camp
+
+GC_DECLARE_PTRFREE(camp::triple);
+
+#endif