#ifndef __JCYLINDER3D__ #define __JCYLINDER3D__ #include #include #include #include #include #include "JIO/JSerialisable.hh" #include "JLang/JManip.hh" #include "JMath/JConstants.hh" #include "JGeometry2D/JCircle2D.hh" #include "JGeometry3D/JVector3D.hh" #include "JGeometry3D/JAxis3D.hh" /** * \author mdejong */ namespace JGEOMETRY3D {} namespace JPP { using namespace JGEOMETRY3D; } namespace JGEOMETRY3D { using JIO::JReader; using JIO::JWriter; using JGEOMETRY2D::JVector2D; using JGEOMETRY2D::JCircle2D; /** * Cylinder object. * * The cylinder consists of a 2D circle in the (x,y) plane and a range in z * (i.e.\ axis of cylinder is parallel to the z-axis). */ class JCylinder3D : public JCircle2D { public: /** * Type definition of intersection. */ typedef std::pair intersection_type; /** * Default constructor. */ JCylinder3D() : JCircle2D(), zmin(0.0), zmax(0.0) {} /** * Constructor. * * \param circle 2D circle in (x,y) * \param zmin minimal z * \param zmax maximal z */ JCylinder3D(const JCircle2D& circle, const double zmin, const double zmax) : JCircle2D(circle) { this->zmin = zmin; this->zmax = zmax; } /** * Constructor. * * Determines smallest enclosing cylinder for any number of points. * * \param __begin begin of data * \param __end end of data * \param precision precision */ template JCylinder3D(T __begin, T __end, const double precision = std::numeric_limits::epsilon()) : JCircle2D(__begin, __end, precision), zmin(0.0), zmax(0.0) { if (__begin != __end) { zmin = std::numeric_limits::max(); zmax = std::numeric_limits::lowest(); for (T i = __begin; i != __end; ++i) { if (i->getZ() < zmin) zmin = i->getZ(); if (i->getZ() > zmax) zmax = i->getZ(); } } } /** * Get minimal z position. * * \return minimal z position */ double getZmin() const { return zmin; } /** * Get maximal z position. * * \return maximal z position */ double getZmax() const { return zmax; } /** * Set minimal z position. * * \param zmin minimal z position */ void setZmin(const double zmin) { this->zmin = zmin; } /** * Set maximal z position. * * \param zmax maximal z position */ void setZmax(const double zmax) { this->zmax = zmax; } /** * Add position. * * \param pos position * \return this cylinder */ JCylinder3D& add(const JVector3D& pos) { static_cast(*this).add(JPosition2D(pos.getX(), pos.getY())); zmin += pos.getZ(); zmax += pos.getZ(); return *this; } /** * Add (safety) margin. * * \param D margin */ void addMargin(const double D) { __r += D; zmin -= D; zmax += D; } /** * Get volume. * * \return volume */ inline double getVolume() const { return (getZmax() - getZmin()) * JMATH::PI * getRadius() * getRadius(); } /** * Get centre. * * \return centre */ JPosition3D getCenter() const { return JPosition3D(getPosition(), (getZmax() - getZmin())/2.0); } /** * Check whether given point is inside cylinder. * * \param pos position * \return true if inside; else false */ inline bool is_inside(const JVector3D& pos) const { return (pos.getZ() >= getZmin() && pos.getZ() <= getZmax() && JCircle2D::is_inside(JVector2D(pos.getX(), pos.getY()))); } /** * Get distance between cylinder wall and given position. * * \param pos position * \return distance */ inline double getDistance(const JVector3D& pos) const { JVector2D D(pos); D.sub(*this); double R = D.getLength(); if (R > this->getRadius()) { R -= this->getRadius(); double dz = 0.0; if (pos.getZ() > this->getZmax()) dz = pos.getZ() - this->getZmax(); else if (pos.getZ() < this->getZmin()) dz = this->getZmin() - pos.getZ(); else return R; return sqrt(R*R + dz*dz); } else { if (pos.getZ() > this->getZmax()) return pos.getZ() - this->getZmax(); else if (pos.getZ() < this->getZmin()) return this->getZmin() - pos.getZ(); else return 0.0; } } /** * Get square of distance between cylinder wall and given position. * * \param pos position * \return square of distance */ inline double getDistanceSquared(const JVector3D& pos) const { const double d = getDistance(pos); return d*d; } /** * Get intersection points of axis with cylinder. * * \param axis axis * \return up and down stream positions along axis */ inline intersection_type getIntersection(const JAxis3D& axis) const { double path[] = { 0.0, 0.0 }; if (fabs(axis.getDZ()) != 0.0) { // intersection with bottom or top surface const double Z[] = { axis.getDZ() > 0 ? this->getZmin() : this->getZmax(), axis.getDZ() > 0 ? this->getZmax() : this->getZmin() }; for (int i = 0; i != 2; ++i) { const double u = (Z[i] - axis.getZ()) / axis.getDZ(); const double x = axis.getX() + u * axis.getDX() - this->getX(); const double y = axis.getY() + u * axis.getDY() - this->getY(); if (x*x + y*y <= this->getRadius() * this->getRadius()) { path[i] = u; } } } if (fabs(axis.getDZ()) != 1.0) { // intersection with cylinder wall const double x = axis.getX() - this->getX(); const double y = axis.getY() - this->getY(); const double dx = axis.getDX(); const double dy = axis.getDY(); const double R = this->getRadius(); const double a = (dx * dx + dy * dy); const double b = 2*(dx * x + dy * y); const double c = (x * x + y * y) - R * R; const double q = b*b - 4*a*c; if (q >= 0) { const double u[] = { (-b - sqrt(q)) / (2*a), (-b + sqrt(q)) / (2*a) }; for (int i = 0; i != 2; ++i) { const double z = axis.getZ() + u[i] * axis.getDZ(); if (z >= this->getZmin() && z <= this->getZmax()) { path[i] = u[i]; } } } } return std::make_pair(path[0], path[1]); } /** * Read cylinder from input stream. * * \param in input stream * \param cylinder cylinder * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JCylinder3D& cylinder) { in >> static_cast(cylinder); in >> cylinder.zmin >> cylinder.zmax; return in; } /** * Write cylinder to output stream. * * \param out output stream * \param cylinder cylinder * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JCylinder3D& cylinder) { const JFormat format(out, getFormat(JFormat_t(9, 3, std::ios::fixed | std::ios::showpos))); out << static_cast(cylinder); out << ' '; out << format << cylinder.zmin; out << ' '; out << format << cylinder.zmax; return out; } /** * Read cylinder from input. * * \param in reader * \param cylinder cylinder * \return reader */ friend inline JReader& operator>>(JReader& in, JCylinder3D& cylinder) { in >> static_cast(cylinder); in >> cylinder.zmin >> cylinder.zmax; return in; } /** * Write cylinder to output. * * \param out writer * \param cylinder cylinder * \return writer */ friend inline JWriter& operator<<(JWriter& out, const JCylinder3D& cylinder) { out << static_cast(cylinder); out << cylinder.zmin << cylinder.zmax; return out; } protected: double zmin; double zmax; }; } #endif