JCylinder3D.hh

Go to the documentation of this file.

#ifndef __JCYLINDER3D__
#define __JCYLINDER3D__

#include <istream>
#include <ostream>
#include <limits>
#include <utility>
#include <cmath>

#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<double, double> 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<class T>
    JCylinder3D(T __begin,
                T __end,
                const double precision = std::numeric_limits<double>::epsilon()) :
      JCircle2D(__begin, __end, precision),
      zmin(0.0),
      zmax(0.0)
    {
      if (__begin != __end) {

        zmin = std::numeric_limits<double>::max();
        zmax = std::numeric_limits<double>::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<JPosition2D&>(*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<JCircle2D&>(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<JCylinder3D>(JFormat_t(9, 3, std::ios::fixed | std::ios::showpos)));

      out << static_cast<const JCircle2D&>(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<JCircle2D&>(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<const JCircle2D&>(cylinder);
      out << cylinder.zmin << cylinder.zmax;

      return out;
    }
    

  protected:
    double zmin;
    double zmax;
  };
}

#endif