Cylinder.cc

Go to the documentation of this file.

#include <esg/geometry/Cylinder.h>
#include <esg/geometry/Sphere.h>

using namespace esg;

void Cylinder::_duplicate_attributes(const Geometry& src)
{
    Geometry::_duplicate_attributes(src);
    _radius = ((Cylinder&)src)._radius;
    _radiusSqr = ((Cylinder&)src)._radiusSqr;
    _cap1Point.set(((Cylinder&)src)._cap1Point);
    _cap1Normal.set(((Cylinder&)src)._cap1Normal);
    _cap2Point.set(((Cylinder&)src)._cap2Point);
    _cap2Normal.set(((Cylinder&)src)._cap2Normal);
    _axis.set(((Cylinder&)src)._axis);
}

void Cylinder::_rotateX(float a)
{
    Matrix3 trMat;
    trMat.rotX(a);
    trMat.transform(_cap1Point);
    trMat.transform(_cap1Normal);
    trMat.transform(_cap2Point);
    trMat.transform(_cap2Normal);
    trMat.transform(_axis);
}

void Cylinder::_rotateY(float a)
{
    Matrix3 trMat;
    trMat.rotY(a);
    trMat.transform(_cap1Point);
    trMat.transform(_cap1Normal);
    trMat.transform(_cap2Point);
    trMat.transform(_cap2Normal);
    trMat.transform(_axis);
}

void Cylinder::_rotateZ(float a)
{
    Matrix3 trMat;
    trMat.rotZ(a);
    trMat.transform(_cap1Point);
    trMat.transform(_cap1Normal);
    trMat.transform(_cap2Point);
    trMat.transform(_cap2Normal);
    trMat.transform(_axis);
}

void Cylinder::_rotate(float a, const Vector3& axis)
{
    Matrix4 trMat;
    trMat.rotationGL(a, axis);
    trMat.transform(_cap1Point);
    trMat.transform(_cap1Normal);
    trMat.transform(_cap2Point);
    trMat.transform(_cap2Normal);
    trMat.transform(_axis);
}

void Cylinder::_rotate(const Matrix3& rMat)
{
    rMat.transform(_cap1Point);
    rMat.transform(_cap1Normal);
    rMat.transform(_cap2Point);
    rMat.transform(_cap2Normal);
    rMat.transform(_axis);
}

void Cylinder::_translate(float x, float y, float z)
{
    _cap1Point.add(Vector3(x, y, z));
    _cap2Point.add(Vector3(x, y, z));
}

void Cylinder::_transform(const Matrix4& trMat)
{
    Matrix3 rMat;
    Vector3 tVec;
    double sc = trMat.get(rMat, tVec);
    trMat.transform(_cap1Point);
    trMat.transform(_cap2Point);
    rMat.transform(_cap1Normal);
    rMat.transform(_cap2Normal);
    _axis.set(_cap2Point);
    _axis.sub(_cap1Point);
    _axis.normalize();
    if (sc != 1.0) {
        _radius *= sc;
        _radiusSqr = _radius * _radius;
    }
}

void Cylinder::_scale(float s)
{
    _cap1Point.scale(s);
    _cap2Point.scale(s);
    _radius *= s;
    _radiusSqr = _radius * _radius;
}


//----- public ------

Cylinder::Cylinder()
    : _radius(1.0),
      _radiusSqr(1.0),
      _cap1Point(0.0, 0.0, 0.0),
      _cap1Normal(0.0, -1.0, 0.0),
      _cap2Point(0.0, 1.0, 0.0),
      _cap2Normal(0.0, 1.0, 0.0),
      _axis(0.0, 1.0, 0.0)
{
}

Cylinder::Cylinder(const Vector3& point1,
                   const Vector3& normal1,
                   const Vector3& point2,
                   const Vector3& normal2,
                   double         radius)
    : _radius(radius),
      _radiusSqr(radius * radius),
      _cap1Point(point1),
      _cap1Normal(normal1),
      _cap2Point(point2),
      _cap2Normal(normal2)
{
    _axis.set(_cap2Point);
    _axis.sub(_cap1Point);
    _axis.normalize();
}

void Cylinder::rayIntersection(PointEnv*      pPE,
                               int            mask,
                               const Vector3& origin,
                               const Vector3& direction,
                               float          maxD) 
{
    pPE->mask = ENV_HAVE_NOTHING;
    
    /*
     * Algorithm reference:
     *    Heckbert P. S., Graphics Gems IV, Academic Press Professional,
     *    London, 1994, p. 356-365
     */
    
    static Vector3 RC, D, O, n;
    double  d, t, s;
    double  t_in =   MAXDOUBLE;
    double  t_out = -MAXDOUBLE;

    RC.set(origin);
    RC.sub(_cap1Point);
    n.cross(direction, _axis);
    double ln = n.length();

    // get intersection with cylinder
    if (ln == 0.0) { // ray is parallel to cylinder
        d = RC.dot(_axis);
        D.set(_axis);
        D.scaleAdd(-d, RC);
        d = D.length();
        if (d > _radius) return;
    } else {
        n.scale(1.0/ln); // normalize
        d = fabs(RC.dot(n)); // shortest distance
        if (d > _radius) return;
        O.cross(RC, _axis);
        t = -O.dot(n) / ln;
        O.cross(n, _axis);
        O.normalize();
        s = fabs(sqrt(_radiusSqr - d * d) / direction.dot(O));
        t_in  = t - s;
        t_out = t + s;
    }

    int surface = 0;

    // clip intersections by end-cap planes
    
    /*  Intersect the ray with the bottom end-cap plane.                */
    double dc = _cap1Normal.dot(direction);
    double dw = _cap1Normal.dot(origin) + _cap1Point.dot(_cap1Normal);

    if  (dc == 0.0) { // If parallel to bottom plane
        if (dw >= 0.0) return;
    } else {
        t  = -dw / dc;
        if (dc >= 0.0) {                            /* If far plane     */
            if (t < t_in) return;
            if (t > t_in && t < t_out) t_out = t;
        } else {                                    /* If near plane    */
            if (t > t_out) return;
            if (t > t_in && t < t_out) {
                t_in = t;
                surface = 1;
            }
        }
    }

    /* Intersect the ray with the top end-cap plane.                    */
    dc = _cap2Normal.dot(direction);
    dw = _cap2Normal.dot(origin) + _cap2Point.dot(_cap2Normal);

    if  (dc == 0.0) {           // If parallel to top plane
        if (dw >= 0.0) return;
    } else {
        t  = -dw / dc;
        if (dc >= 0.0) {                            /* If far plane     */
            if (t < t_in) return;
            if (t > t_in && t < t_out) t_out = t;
        } else {                                    /* If near plane    */
            if (t > t_out) return;
            if (t > t_in && t < t_out) {
                t_in = t;
                surface = 2;
            }
        }
    }

    if (t_in > t_out) return;

    if (mask & ENV_WANT_INTERSECTION) {
        pPE->intersection.set(direction);
        pPE->intersection.scaleAdd(t_in, origin);
        pPE->mask |= ENV_HAVE_INTERFERENCE | ENV_HAVE_INTERSECTION;
    } else
        pPE->mask |= ENV_HAVE_INTERFERENCE;

    if (mask & ENV_WANT_NORMAL) {
        switch (surface) {
            case 1: pPE->normal.set(_cap1Normal); break;
            case 2: pPE->normal.set(_cap2Normal); break;
            default:
                {
                    double d1 = pPE->intersection.dot(_axis);
                    double d2 = _cap1Point.dot(_axis);
                    Vector3 v(_axis);
                    v.scale(d2 - d1);
                    v.add(_cap1Point);
                    pPE->normal.set(pPE->intersection);
                    pPE->normal.sub(v);
                    pPE->normal.normalize();
                }
        }
        pPE->normalOrientation = PointEnv::OUTWARDS_NORMAL;
        pPE->mask |= ENV_HAVE_NORMAL;
    }
    
    if (mask & ENV_WANT_DISTANCE) {
        pPE->distance = t_in;
        pPE->mask |= ENV_HAVE_DISTANCE;
    }
}

bool Cylinder::mapToUV(const Vector3& v, Vector2& uv)
{
    cerr << "Cylinder::mapToUV(): Not implemented" << endl;
    return false;
}

void Cylinder::randomSample(int mask, PointEnv& pe, double* pdf)
{
    pe.mask = ENV_HAVE_NOTHING;
    cerr << "Cylinder::randomSample(): Not implemented" << endl;
}

Interval Cylinder::extent(const Vector3& direction) const
{
    switch ((int) direction.dot(_axis)) {
        case 0: // perpendicular to axis
            {
                double e = _cap1Point.dot(direction);
                return Interval(e - _radius, e + _radius);
            }
        case 1:  // parellel to axis
        case -1:
            {
                double e1 = _cap1Point.length();
                double e2 = _cap2Point.length();
                return ((e1 < e2)
                        ? Interval(e1, e2)
                        : Interval(e2, e1));
            }
        default:
            {
                double e1 = _cap1Point.length();
                double e2 = _cap2Point.length();
                return ((e1 < e2)
                        ? Interval(e1 - _radius, e2 + _radius)
                        : Interval(e2 - _radius, e1 + _radius));
            }
    }
}

Vector3 Cylinder::centroid(void) const
{
    Vector3 c(_cap1Point);
    c.add(_cap2Point);
    c.scale(0.5);
    return c;
}

double Cylinder::radius(const Vector3& centroid) const
{
    Vector3 v(_cap1Point);
    v.sub(centroid);
    double r1 = v.length();
    v.set(_cap2Point);
    v.sub(centroid);
    double r2 = v.length();
    return ((r1 > r2) ? r1 + _radius : r2 + _radius);
}

double Cylinder::radius(void) const
{
    Vector3 v(_cap1Point);
    v.sub(_cap2Point);
    double l = v.length()/2.0;
    return sqrt(l*l + _radiusSqr);
}

Geometry* Cylinder::clone(const Matrix4* pTrMat) const
{
    Cylinder* pRet = new Cylinder();
    pRet->_duplicate_attributes(*this);
    if (pTrMat) pRet->_transform(*pTrMat);
    return pRet;
}

bool Cylinder::separation(Geometry& geom, Vector3*  pDir)
{
    Interval myInt(_cap1Point.dot(_axis), _cap2Point.dot(_axis));
    Interval gInt(geom.extent(_axis));

    // separation by the end-cap planes
    if (myInt.min > myInt.max) {
        double a = myInt.min;
        myInt.min = myInt.max;
        myInt.max = a;
    }
    if (gInt.min > myInt.max || gInt.max < myInt.min) {
        if (pDir) pDir->set(_axis);
        return true;
    }

    // separation by the "pipe"
    Vector3 v(geom.centroid());
    v.sub(_cap1Point);
    v.normalize();
    Vector3 aux;
    aux.cross(_axis, v);
    aux.normalize();
    v.cross(aux, _axis);
    v.normalize();
    gInt = geom.extent(v);
    double d = _cap1Point.dot(v);
    if (gInt.min > d + _radius || gInt.max < d - _radius) {
        if (pDir) pDir->set(_axis);
        return true;
    }

    // no gap found
    return false;
}

double Cylinder::distance(const Geometry& geom, Vector3*  pDir)
{
    Interval myInt(_cap1Point.dot(_axis), _cap2Point.dot(_axis));
    Interval gInt(geom.extent(_axis));
    double   d;
    double   ret = -MAXDOUBLE;

    // distance of the end-cap planes
    if (myInt.min > myInt.max) {
        double a = myInt.min;
        myInt.min = myInt.max;
        myInt.max = a;
    }
    d = gInt.min - myInt.max;
    if (d > ret) {
        ret = d;
        if (pDir) pDir->set(_axis);
    }
    d = myInt.min - gInt.max;
    if (d > ret) {
        ret = d;
        if (pDir) pDir->set(_axis);
    }

    // distance of the "pipe"
    Vector3 v(geom.centroid());
    v.sub(_cap1Point);
    v.normalize();
    Vector3 aux;
    aux.cross(_axis, v);
    aux.normalize();
    v.cross(aux, _axis);
    v.normalize();
    gInt = geom.extent(v);
    double md = _cap1Point.dot(v);
    
    d = gInt.min - (md + _radius);
    if (d > ret) {
        ret = d;
        if (pDir) pDir->set(v);
    }
    d = (md - _radius) - gInt.max;
    if (d > ret) {
        ret = d;
        if (pDir) pDir->set(v);
    }


    // distance of bounding spheres
    v.set(_cap1Point);
    v.sub(_cap2Point);
    double l = v.length()/2.0;
    Sphere bs(centroid(), radius());
    d = bs.distance(geom, &v);
    if (d > ret) {
        ret = d;
        if (pDir) pDir->set(v);
    }

    return ret;
}

void Cylinder::dump(const char* intent, const char* tab)
{
}

---