From 4d5b8732c97c73ce7cf04e6641a239ceb0447d8d Mon Sep 17 00:00:00 2001 From: Franklin Wei Date: Sun, 3 Feb 2019 17:51:06 -0500 Subject: Add circular arc and spiral curves --- curve.cpp | 89 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ curve.h | 36 +++++++++++++++++++++++--- quat.cpp | 38 +++++++++++++++++++++++++++ quat.h | 25 ++++++++++++++++++ vec3.h | 3 +++ 5 files changed, 188 insertions(+), 3 deletions(-) create mode 100644 quat.cpp create mode 100644 quat.h diff --git a/curve.cpp b/curve.cpp index 9b63e8d..73f40c1 100644 --- a/curve.cpp +++ b/curve.cpp @@ -1,3 +1,4 @@ +#include #include #include "curve.h" @@ -21,3 +22,91 @@ vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl) return sum; } + +vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl) +{ + //cout << "Integrating loop of length " << this->angle << " radians" << endl; + + double r = this->radius.magnitude(), dtheta = dl / r; + + //cout << "R = " << r << ", dtheta = " << dtheta << endl; + + quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate(); + + //cout << "rot = " << rot << ", crot = " << crot << endl; + + double len = this->angle * r, l; + + vec3 sum = 0; + + quat radius = this->radius; + + for(l = 0; l < len; l += dl) + { + vec3 ds = this->normal.cross(radius).normalize() * dl; + sum += integrand(this->center + radius, ds); + + radius = rot * radius * crot; + } + + if(l < len) + { + dl = len - l; + dtheta = dl / r; + rot = quat::from_angleaxis(dtheta, normal); + crot = rot.conjugate(); + + vec3 ds = this->normal.cross(radius).normalize() * dl; + sum += integrand(this->center + radius, ds); + } + + return sum; +} + +vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl) +{ + //cout << "Integrating loop of length " << this->angle << " radians" << endl; + + double r = this->radius.magnitude(), dtheta = dl / r; + + //cout << "R = " << r << ", dtheta = " << dtheta << endl; + + quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate(); + + //cout << "rot = " << rot << ", crot = " << crot << endl; + + double len = this->angle * r, l; + + vec3 sum = 0; + + quat radius = this->radius; + + /* how far along the axis we've moved */ + vec3 axis = 0; + + /* we add this axial component each iteration */ + vec3 dp = this->normal * dtheta * pitch / (2 * M_PI); + + for(l = 0; l < len; l += dl, axis += dp) + { + vec3 ds = this->normal.cross(radius).normalize() * dl + dp; + sum += integrand(this->origin + axis + radius, ds); + + /* rotate by dtheta about normal */ + radius = rot * radius * crot; + } + + if(l < len) + { + dl = len - l; + dtheta = dl / r; + rot = quat::from_angleaxis(dtheta, normal); + crot = rot.conjugate(); + dp = this->normal * dtheta * pitch / (2 * M_PI); + + vec3 ds = this->normal.cross(radius).normalize() * dl + dp; + sum += integrand(this->origin + axis + radius, ds); + } + + return sum; +} diff --git a/curve.h b/curve.h index c896a23..e9e9642 100644 --- a/curve.h +++ b/curve.h @@ -1,5 +1,10 @@ +#ifndef CURVE_H +#define CURVE_H + #include #include "vec3.h" +#include "quat.h" + using namespace std; class Curve { @@ -19,10 +24,35 @@ public: class Arc : Curve { private: vec3 center; - double radius; - double angle[2]; /* start and end angles */ + + /* these are relative to the center (direction will be determined + * by RHR of normal), and should be orthonormal */ + vec3 radius, normal; + + /* how many radians the arc extends for (can be greater than 2pi) */ + double angle; public: - LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {}; + Arc(vec3 c_, vec3 r_, vec3 n_, double th) : center(c_), radius(r_), normal(n_), angle(th) {}; + + vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta); +}; + +class Spiral : Curve { +private: + vec3 origin; + + /* these are relative to the center (direction will be determined + * by RHR of normal), and should be orthonormal */ + vec3 radius, normal; + + /* how many radians the arc extends for (can be greater than 2pi) */ + double angle; + + /* space between turns (2pi) */ + double pitch; +public: + Spiral(vec3 c_, vec3 r_, vec3 n_, double th, double p) : origin(c_), radius(r_), normal(n_), angle(th), pitch(p) {}; vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta); }; +#endif diff --git a/quat.cpp b/quat.cpp new file mode 100644 index 0000000..aa53af6 --- /dev/null +++ b/quat.cpp @@ -0,0 +1,38 @@ +#include "quat.h" +#include + +quat::quat(double w_, double x_, double y_, double z_) : w(w_), x(x_), y(y_), z(z_) { } +quat::quat(double x_, double y_, double z_) : w(0), x(x_), y(y_), z(z_) { } +quat::quat(double w_, vec3 vec) : w(w_), x(vec[0]), y(vec[1]), z(vec[2]) { } +quat::quat(vec3 vec) : w(0), x(vec[0]), y(vec[1]), z(vec[2]) { } +quat::quat() : w(0), x(0), y(0), z(0) { } + +quat::operator vec3() +{ + return vec3(this->x, this->y, this->z); +} + +quat operator*(const quat &lhs, const quat &rhs) +{ + return quat(lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z, + lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y, + lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x, + lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w); +} + +quat quat::conjugate() const +{ + return quat(this->w, -this->x, -this->y, -this->z); +} + +quat quat::from_angleaxis(double angle, vec3 axis) +{ + double si = std::sin(angle / 2); + double co = std::cos(angle / 2); + return quat(co, si * axis[0], si * axis[1], si * axis[2]); +} + +std::ostream &operator<<(std::ostream &os, const quat &q) +{ + return os << "(" << q.w << ", " << q.x << ", " << q.y << ", " << q.z << ")"; +} diff --git a/quat.h b/quat.h new file mode 100644 index 0000000..65a828a --- /dev/null +++ b/quat.h @@ -0,0 +1,25 @@ +#ifndef QUAT_H +#define QUAT_H +#include "vec3.h" +#include + +class quat { +public: + double w, x, y, z; +public: + quat(double w, double x, double y, double z); + quat(double x, double y, double z); + quat(double w, vec3 vec); + quat(vec3 vec); + quat(); + + operator vec3(); + + quat conjugate() const; + + static quat from_angleaxis(double angle, vec3 axis); +}; + +quat operator*(const quat &, const quat &); +std::ostream &operator<<(std::ostream &os, const quat &); +#endif diff --git a/vec3.h b/vec3.h index 1f01098..abfc6c3 100644 --- a/vec3.h +++ b/vec3.h @@ -1,3 +1,5 @@ +#ifndef VEC3_H +#define VEC3_H #include class vec3 { public: @@ -25,3 +27,4 @@ class vec3 { }; std::ostream &operator<<(std::ostream &output, const vec3 &v); +#endif -- cgit v1.1