Sphere.cc

Go to the documentation of this file.

#include <esg/geometry/Sphere.h>
#include <esg/mesh/SphereMesh.h>
#include <esg/explorer/POGExplorer.h>
#include <esg/explorer/RadiusExplorer.h>

using namespace esg;

Mesh* Sphere::_mesh(int density) const
{
    Mesh*  pMesh = (SphereMesh*) new SphereMesh(density);
    if (!pMesh) return NULL;
    if (_radius != 1.0)
        pMesh->scale(_radius, _radius, _radius);
    if (_centre.x || _centre.y || _centre.z)
        pMesh->translate(_centre);
    return pMesh;
}

void Sphere::_duplicate_attributes(const Geometry& src)
{
    _radius   = ((Sphere&)src)._radius;
    _radius2  = _radius * _radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
    _centre.set(((Sphere&)src)._centre);
}

void Sphere::_rotateX(float a)
{
    Matrix3 trMat;
    trMat.rotX(a);
    trMat.transform(_centre);
}

void Sphere::_rotateY(float a)
{
    Matrix3 trMat;
    trMat.rotY(a);
    trMat.transform(_centre);
}

void Sphere::_rotateZ(float a)
{
    Matrix3 trMat;
    trMat.rotZ(a);
    trMat.transform(_centre);
}

void Sphere::_rotate(float a, const Vector3& axis)
{
    Matrix4 trMat;
    trMat.rotationGL(a, axis);
    trMat.transform(_centre);
}

void Sphere::_rotate(const Matrix3& rotMat)
{
    rotMat.transform(_centre);
}

void Sphere::_translate(float x, float y, float z)
{
    _centre.add(Vector3(x, y, z));
}

void Sphere::_transform(const Matrix4& trMat)
{
    trMat.transform(_centre);
}

void Sphere::_scale(float s)
{
    _centre.scale(s);
    _radius *= s;
    _radius2  = _radius * _radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}


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

Sphere::Sphere()
{
    _radius   = 1.;
    _radius2  = 1.;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
    _centre.set(0., 0., 0.);
}

Sphere::Sphere(const Vector3& c, float r)
{
    _radius   = r;
    _radius2  = r * r;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
    _centre.set(c);
}


Sphere::Sphere(float x, float y, float z, float r)
{
    _radius   = r;
    _radius2  = r * r;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
    _centre.set(x,y,z);
}

Sphere::Sphere(const Geometry& g)
{
    _centre.set(g.centroid());
    _radius   = g.radius(_centre);
    _radius2  = _radius * _radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}

Sphere::Sphere(const Geometry& g1, const Geometry& g2)
{
    Vector3 c1(g1.centroid());
    Vector3 c2(g2.centroid());
    Vector3 v(c2 - c1);
    double  r1 = g1.radius(c1);
    double  r2 = g2.radius(c2);
    double  d  = v.length();

    if ((d + r2) <= r1) { // second inside first
        _centre.set(c1);
        _radius   = r1;
        _radius2  = r1 * r1;
    } else {
        if ((d + r1) <= r2) { // first inside second
            _centre.set(c2);
            _radius    = r2;
            _radius2   = r2 * r2;
        } else {
            c2.scale((1.0 + (r2 - r1) / d) / 2.0);
            _centre.set(c1);
            _centre.scaleAdd((1.0 + (r1 - r2) / d) / 2.0, c2);
            _radius    = (r1 + r2 + d) / 2.0;
            _radius2   = _radius * _radius;
        }
    }

    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}

Sphere::Sphere(SceneGraphObject& obj)
{
    POGExplorer explorer;
    explorer.explore(obj);
    
    _centre.set(explorer.result());
    
    RadiusExplorer explorer2(_centre);
    explorer2.explore(obj);
    
    _radius   = explorer2.result();
    _radius2  = _radius * _radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}

Sphere::Sphere(List<SceneGraphObject>& list)
{
    _centre.set(0,0,0);
    
    unsigned          num = 0;
    SceneGraphObject* pObject = list.firstItem();
    while (pObject) {
        POGExplorer pogExplorer;
        pogExplorer.explore(*pObject);
        _centre.add(pogExplorer.result());
        num++;
        pObject = list.nextItem();
    }
    
    if (num) _centre.scale(1.0/num);

    _radius = -MAXFLOAT;

    float r;
    pObject = list.firstItem();
    while (pObject) {
        RadiusExplorer radExplorer(_centre);
        radExplorer.explore(*pObject);
        r = radExplorer.result();
        if (r > _radius) _radius = r;
        pObject = list.nextItem();
    }
    _radius2  = _radius * _radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}

void Sphere::rayIntersection(PointEnv*      pPE,
                             int            mask,
                             const Vector3& origin,
                             const Vector3& direction,
                             float          maxDist)
{
    if (!pPE) {
        fprintf(stderr,"Sphere error: No point environment structure defined\n");
        return;
    }
 
    pPE->mask = ENV_HAVE_NOTHING;

    Vector3 v(_centre - origin);
    double  v2 = v.dot(v);
    double& t0 = pPE->distance;
    
    t0 = v.dot(direction);
    if (v2 > _epsMinus) { // outside of or at the sphere
        if (t0 < Geometry::EPS) return; // sphere is 'behind' the ray
        register double td2 = _radius2 - v2 + (t0*t0);
        if (td2 < 0.0) return; // ray doesn't hit the sphere
        // if ray starts on the surface then get farter (inner) intersection
        t0 = (v2 < _epsPlus) ? t0 + sqrt(td2) : t0 - sqrt(td2);
    } else  // ray starts inside the sphere
        t0 += sqrt(_radius2 - v2 + (t0*t0));

    if (t0 > maxDist) return;

    if (mask & (ENV_WANT_INTERSECTION|ENV_WANT_UV_COORD|ENV_WANT_NORMAL)) {
        pPE->intersection.set(direction);
        pPE->intersection.scaleAdd(t0, origin);
        if (mask & ENV_WANT_UV_COORD) {
            Vector3 xyz(pPE->intersection);
            xyz.sub(_centre);
            double auxV = acos(xyz.z / _radius);
            if (xyz.y >= 0.) 
                pPE->uvCoord.set(acos(xyz.x/(_radius*sin(auxV))) / (2*PI),
                                 auxV / PI);
            else
                pPE->uvCoord.set((PI+acos(xyz.x/(_radius*sin(auxV)))) / (2*PI),
                                 auxV / PI);
            pPE->mask |=
                ENV_HAVE_INTERFERENCE|ENV_HAVE_INTERSECTION|ENV_HAVE_DISTANCE|ENV_HAVE_UV_COORD;
        } else 
            pPE->mask |=
                ENV_HAVE_INTERFERENCE|ENV_HAVE_INTERSECTION|ENV_HAVE_DISTANCE;
        if (mask & ENV_WANT_NORMAL) {
            pPE->normal.set(pPE->intersection - _centre);
            //pPE->normal.scale(1./_radius); // unsharp normalization 
            pPE->normal.normalize();
            pPE->normalOrientation = PointEnv::OUTWARDS_NORMAL;
            pPE->mask |= ENV_HAVE_NORMAL;
        }
    } else 
        pPE->mask |= ENV_HAVE_INTERFERENCE|ENV_HAVE_DISTANCE;
}

bool Sphere::mapToUV(const Vector3& v, Vector2& uv)
{
    Vector3 xyz(v);
    xyz.sub(_centre);
    xyz.normalize();

    // re-arrange for pole distortion
    double aux = xyz.x; xyz.x = xyz.z; xyz.z = aux; 

    double r   = sqrt((xyz.x * xyz.x) + (xyz.y * xyz.y));
    double rho = sqrt((r * r) + (xyz.z * xyz.z));
    double theta, phi;

    if(r == 0.0) {
        theta = 0.0;
        phi   = 0.0;
    } else {
        phi   = (xyz.z == 0.0) ? PI / 2.0 : acos(xyz.z / rho);
        theta = (xyz.y == 0.0) ? PI / 2.0 : asin(xyz.y / r) + (PI / 2.0);
    }

    uv.set(1 - (phi / PI), 1 - (theta / PI));

    return true;
    
#if 0 // funkci mapovani na dve polokoule
    Vector3 xyz(v);
    xyz.sub(_centre);
    xyz.normalize();

    // re-arrange for pole distortion
    double aux = xyz.x; xyz.x = xyz.z; xyz.z = aux; 

    double r   = sqrt((xyz.x * xyz.x) + (xyz.y * xyz.y));
    double rho = sqrt((r * r) + (xyz.z * xyz.z));
    double theta, phi;

    if(r == 0.0) {
        theta = 0.0;
        phi   = 0.0;
    } else {
        phi   = (xyz.z == 0.0) ? PI / 2.0 : acos(xyz.z / rho);
        theta = (xyz.y == 0.0) ? PI / 2.0 : asin(xyz.y / r) + (PI / 2.0);
    }

    uv.set(1 - (phi / PI), 1 - (theta / PI));
    return true;
#endif
}

void Sphere::randomSample(int mask, PointEnv& pe, double* pdf) 
{
    pe.mask = ENV_HAVE_NOTHING;
    
    if (mask & (ENV_WANT_SURFACE_POINT|ENV_WANT_NORMAL|ENV_WANT_UV_COORD)) {
        double d;
        do { // rejection sampling
            pe.intersection.set(2. * ESG_DBL_RAND - 1.,
                                2. * ESG_DBL_RAND - 1.,
                                2. * ESG_DBL_RAND - 1.);
            d = pe.intersection.lengthSquared();
        } while (d > 1.);
        pe.intersection.scale(_radius/sqrt(d));
        pe.intersection.add(_centre);
        pe.mask |= ENV_HAVE_SURFACE_POINT;

        if ((mask&ENV_WANT_UV_COORD) && mapToUV(pe.intersection, pe.uvCoord))
            pe.mask |= ENV_HAVE_UV_COORD;
        
        if (mask & ENV_WANT_NORMAL) {
            pe.normal.set(pe.intersection);
            pe.normal.normalize();
            pe.mask |= ENV_HAVE_NORMAL;
            pe.normalOrientation = PointEnv::OUTWARDS_NORMAL;
        }
    }

    if (pdf) *pdf = .25 / (M_PI * _radius2);
}

bool Sphere::randomDirection(const Vector3& pov, Vector3& dir, double* pdf) 
{
    /*
     * Reference: 
     *    Arvo J. et al.:
     *    State of the Art in Monte Carlo Ray Tracing
     *    for Relistic Image Sythesis, 
     *    SIGGRAPH2001 Course 29, Aug 13, 2001. Page 33.
     */

    double  wLength;
    Vector3 u(0.0, 1.0, 0.0);
    Vector3 v;
    Vector3 w(_centre);

    w.sub(pov);
    wLength = w.length();
    w.scale(1./wLength);

    (fabs(w.x) > fabs(w.y)) ? u.set(0.0, 1.0, 0.0) : u.set(1.0, 0.0, 0.0);
    
    v.cross(w, u); v.normalize();
    u.cross(w, v); u.normalize();

    double rd       = _radius / wLength;
    double cosAlpha = ESG_DBL_RAND * (sqrt(1. - rd * rd) - 1.) + 1.;
    double phi      = ESG_DBL_RAND * 2 * M_PI;
    double alpha    = acos(cosAlpha);
    double sinAlpha = sin(alpha);

    Vector3 aux(cos(phi) * sinAlpha, sin(phi) * sinAlpha, cosAlpha);
    aux.normalize();
    
    dir.set(aux.x * u.x + aux.y * v.x + aux.z * w.x,
            aux.x * u.y + aux.y * v.y + aux.z * w.y,
            aux.x * u.z + aux.y * v.z + aux.z * w.z);

    // PDF: to do
    if (pdf)
        cerr << "Sphere::randomDirection(): PDF is not implemented" << endl;
    
    return true;
}

Interval Sphere::extent(const Vector3& direction) const
{
  Interval retInt;
  float    hlp;

  hlp = _centre.dot(direction);
  retInt.min = hlp - _radius;
  retInt.max = hlp + _radius;
  return retInt;
}

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

Vector3 Sphere::centroid() const
{
    return _centre;
}

double Sphere::radius() const
{
    return _radius;
}

double Sphere::radius(const Vector3& centroid) const
{
    return (((Vector3)(centroid-_centre)).length() + _radius);
}

bool Sphere::separation(Geometry& geom, Vector3* pDir) 
{
    register double rads = geom.radius() + _radius;
    Vector3 c(geom.centroid());
    c.sub(_centre);
    register double ns = c.normSquared();
    if (ns > rads * rads) {
        ns = sqrt(ns);
        if (pDir) pDir->set(c.x/ns, c.y/ns, c.z/ns); // normalize
        return true;
    }
    return false;
}

double Sphere::distance(const Geometry& geom, Vector3* pDir)
{
    Vector3 c = geom.centroid();
    register double d = _radius + geom.radius(c);
    c.sub(_centre);
    register double d2 = c.length();
    if (pDir) pDir->set(c.x / d2, c.y / d2, c.z / d2);
    return (d2 - d);
}

void Sphere::dump(const char* intent, const char* tab)
{
    printf("%s geometry Sphere {", intent);
    printf("%s %s centre %f %f %f\n", intent, tab, _centre.x, _centre.y, _centre.z);
    printf("%s %s radius %f\n", intent, tab, _radius);
    printf("%s }", intent);
}


void Sphere::setCentre(const Vector3& centre)
{
    _centre.set(centre);
}


void Sphere::setCentre(float x, float y, float z)
{
    _centre.set(x,y,z);
}

void Sphere::setRadius(float radius)
{
    _radius   = radius;
    _radius2  = radius * radius;
    _epsPlus  = (_radius + Geometry::EPS) * (_radius + Geometry::EPS);
    _epsMinus = (_radius - Geometry::EPS) * (_radius - Geometry::EPS);
}

---