Hemisphere.cc

Go to the documentation of this file.

#include <esg/geometry/Hemisphere.h>

using namespace esg;

Mesh* Hemisphere::_mesh(int density) const
{
    fprintf(stderr,"Hemisphere::_mesh(): Unimplemented\n");
    return NULL;
}

void Hemisphere::_duplicate_attributes(const Geometry& src)
{
    Sphere::_duplicate_attributes(src);
    _xAxis.set(((Hemisphere&)src)._xAxis);
    _yAxis.set(((Hemisphere&)src)._yAxis);
    _zAxis.set(((Hemisphere&)src)._zAxis);
    _backfaceCulling = ((Hemisphere&)src)._backfaceCulling;
}

void Hemisphere::_rotateX(float a)
{
    Matrix3 trMat;
    trMat.rotX(a);
    trMat.transform(_centre);
    trMat.transform(_xAxis);
    trMat.transform(_yAxis);
    trMat.transform(_zAxis);
}

void Hemisphere::_rotateY(float a)
{
    Matrix3 trMat;
    trMat.rotY(a);
    trMat.transform(_centre);
    trMat.transform(_xAxis);
    trMat.transform(_yAxis);
    trMat.transform(_zAxis);
}

void Hemisphere::_rotateZ(float a)
{
    Matrix3 trMat;
    trMat.rotZ(a);
    trMat.transform(_centre);
    trMat.transform(_xAxis);
    trMat.transform(_yAxis);
    trMat.transform(_zAxis);
}

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

void Hemisphere::_rotate(const Matrix3& rotMat)
{
    rotMat.transform(_centre);
    rotMat.transform(_xAxis);
    rotMat.transform(_yAxis);
    rotMat.transform(_zAxis);
}

void Hemisphere::_transform(const Matrix4& trMat)
{
    Matrix3 rMat;
    trMat.get(rMat);
    trMat.transform(_centre);
    rMat.transform(_xAxis);
    rMat.transform(_yAxis);
    rMat.transform(_zAxis);
}


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

Hemisphere::Hemisphere(const Vector3& c, float r, const Vector3& z, bool bfc)
    : Sphere(c,r), _zAxis(z), _backfaceCulling(bfc)
{
    (fabs(z.y) > fabs(z.x)) ? _xAxis.set(1., 0., 0.) : _xAxis.set(0., 1., 0.);
    _yAxis.cross(_zAxis, _xAxis);
    _yAxis.normalize();
    _xAxis.cross(_zAxis, _yAxis);
    _xAxis.normalize();
}

Hemisphere::Hemisphere(const Vector3& z, bool bfc)
    : _zAxis(z), _backfaceCulling(bfc)
{
    (fabs(z.y) > fabs(z.x)) ? _xAxis.set(1., 0., 0.) : _xAxis.set(0., 1., 0.);
    _yAxis.cross(_zAxis, _xAxis);
    _yAxis.normalize();
    _xAxis.cross(_zAxis, _yAxis);
    _xAxis.normalize();
}

Hemisphere::Hemisphere(bool bfc)
    : _xAxis(Vector3(1., 0., 0.)),
      _yAxis(Vector3(1., 0.,-1.)),
      _zAxis(Vector3(0., 1., 0.)),
      _backfaceCulling(bfc)
{
}


void Hemisphere::rayIntersection(PointEnv*      pPE,
                                 int            mask,
                                 const Vector3& origin,
                                 const Vector3& direction,
                                 float          maxDist)
{
    pPE->mask = ENV_HAVE_NOTHING;

    Vector3          v(_centre - origin);
    register double  v2 = v.dot(v);
    register double& t0 = pPE->distance;

    fprintf(stderr,"Hemisphere::_rayIntersection(): Unimplemented\n");
    return; // TO DO

    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 Hemisphere::mapToUV(const Vector3& v, Vector2& uv)
{
    fprintf(stderr,"Hemisphere::mapToUV(): Unimplemented\n");
    return false; // TO DO 
}

void Hemisphere::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)) {
        pe.intersection.set(sampleUniformly(ESG_DBL_RAND, ESG_DBL_RAND, pdf));
        pe.intersection.scale(_radius/pe.intersection.length());
        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 *= _radius2;
}

bool Hemisphere::randomDirection(const Vector3& pov, Vector3&, double* pdf) 
{
    fprintf(stderr,"Hemisphere::randomDirection(): Not implemented\n");
    return false;
}

Interval Hemisphere::extent(const Vector3& direction) const
{
    fprintf(stderr,"Hemisphere::extent(): Unimplemented\n");
    return Interval(FLT_MAX,FLT_MIN);   //MAXFLOAT,MINFLOAT);
}

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

double Hemisphere::radius(const Vector3& centroid) const
{
    fprintf(stderr,"Hemisphere::radius(): Unimplemented\n");
    return - MAXDOUBLE; // TO DO
}

bool Hemisphere::separation(Geometry& geom, Vector3* pDir)
{
    fprintf(stderr,"Hemisphere::separation(): Unimplemented\n");
    return false; // TO DO
}

double Hemisphere::distance(const Geometry& geom, Vector3* pDir)
{
    fprintf(stderr,"Hemisphere::distance(): Unimplemented\n");
    return - MAXDOUBLE; // TO DO
}

void Hemisphere::dump(const char* intent, const char* tab)
{
    printf("%s geometry Hemisphere {\n", 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 %s backfaceCulling %s\n", intent, tab, ((_backfaceCulling) ? "yes" : "no"));
    printf("%s %s orientation [ # x axis, y axis, z axis\n", intent, tab);
    printf("%s %s %s %f %f %f\n", intent, tab, tab, _xAxis.x, _xAxis.y, _xAxis.z);
    printf("%s %s %s %f %f %f\n", intent, tab, tab, _yAxis.x, _yAxis.y, _yAxis.z);
    printf("%s %s %s %f %f %f\n", intent, tab, tab, _zAxis.x, _zAxis.y, _zAxis.z);
    printf("%s %s ]\n", intent, tab);
    printf("%s }\n", intent);
}

void Hemisphere::setZenith(const Vector3& zenith)
{
    _zAxis.set(zenith);
    (fabs(zenith.y) > fabs(zenith.x))
        ? _xAxis.set(1., 0., 0.)
        : _xAxis.set(0., 1., 0.);
    _yAxis.cross(_zAxis, _xAxis);
    _yAxis.normalize();
    _xAxis.cross(_zAxis, _yAxis);
    _xAxis.normalize();
}


Vector2 Hemisphere::world2local2D(const Vector3& dir) const
{
    Vector3 aux(dir.x * _xAxis.x + dir.y * _xAxis.y + dir.z * _xAxis.z,  
                dir.x * _yAxis.x + dir.y * _yAxis.y + dir.z * _yAxis.z,  
                dir.x * _zAxis.x + dir.y * _zAxis.y + dir.z * _zAxis.z); 
    double len = sqrt(aux.x * aux.x + aux.y * aux.y);                 
    double p   = 0.0;                                                 
    if (len != 0.0) {                                                 
        p = acos(aux.x/len);                                          
        if (aux.y < 0.) p = 2 * M_PI - p;                             
    }
    return Vector2(acos(aux.z), p);
}

#define ESG_ANGLES_TO_VECTOR(a, cosE, sinE, d) { \
    double cos_a = cos(a); \
    Vector3 aux(_xAxis); \
    aux.scale(cos_a); \
    (d).set(_yAxis); \
    (d).scaleAdd((a<M_PI) ? sqrt(1-cos_a*cos_a) : -sqrt(1-cos_a*cos_a), aux); \
    aux.set(_zAxis); /* ^^^^^^^^ fast sin(a) ^^^^^^^^^ */  \
    aux.scale(cosE); \
    (d).scaleAdd(sinE, aux); \
}

Vector3 Hemisphere::sampleUniformly(double  elevationRand,
                                    double  azimuthRand,
                                    double* pPDFVal) const
{
    double azimuth = 2. * PI * azimuthRand;
    double cosElev = sqrt(1. - elevationRand); 
    double sinElev = sqrt(elevationRand);

    Vector3 dir;
    
    ESG_ANGLES_TO_VECTOR(azimuth, cosElev, sinElev, dir);

    if (pPDFVal) *pPDFVal = .5 / PI;

    return dir;
}

Vector2 Hemisphere::sampleUniformly2D(double  elevationRand,
                                      double  azimuthRand,
                                      double* pPDFVal) const
{
    if (pPDFVal) *pPDFVal = .5 / PI;
    return Vector2(acos(sqrt(1. - elevationRand)), 2. * PI * azimuthRand);
}

Vector2 Hemisphere::dir2uvUniformly(const Vector3& dir) const
{
    Vector3 aux(dir.x * _xAxis.x + dir.y * _xAxis.y + dir.z * _xAxis.z,  
                dir.x * _yAxis.x + dir.y * _yAxis.y + dir.z * _yAxis.z,  
                dir.x * _zAxis.x + dir.y * _zAxis.y + dir.z * _zAxis.z);

    double len = sqrt(aux.x * aux.x + aux.y * aux.y);                 
    double phi = 0.0;
    if (len != 0.0) {                                                 
        phi = acos(aux.x/len) / (2 * M_PI);
        if (aux.y < 0.) phi = 1. - phi;
    }

    return Vector2(1. - aux.z * aux.z, phi);
    //                  ^^^^^  => cos(acos(aux.z)) = aux.z
}

Vector3 Hemisphere::samplePriorToZenith(double  elevationRand,
                                        double  azimuthRand,
                                        double* pPDFVal) const
{
    double azimuth = 2. * PI * azimuthRand;
    double sinElev = elevationRand;
    double cosElev = sqrt(1. - elevationRand * elevationRand);

    Vector3 dir;
    
    ESG_ANGLES_TO_VECTOR(azimuth, cosElev, sinElev, dir);

    if (pPDFVal) *pPDFVal = cosElev / PI;

    return dir;
}

Vector2 Hemisphere::samplePriorToZenith2D(double  elevationRand,
                                          double  azimuthRand,
                                          double* pPDFVal) const
{
    double cosElev = sqrt(1. - elevationRand * elevationRand);
    if (pPDFVal) *pPDFVal = cosElev / PI;
    return Vector2(asin(elevationRand), 2. * PI * azimuthRand);
}

Vector2 Hemisphere::dir2uvPriorToZenith(const Vector3& dir) const 
{
    Vector3 aux(dir.x * _xAxis.x + dir.y * _xAxis.y + dir.z * _xAxis.z,  
                dir.x * _yAxis.x + dir.y * _yAxis.y + dir.z * _yAxis.z,  
                dir.x * _zAxis.x + dir.y * _zAxis.y + dir.z * _zAxis.z);

    double len = sqrt(aux.x * aux.x + aux.y * aux.y);                 
    double phi = 0.0;                                                 
    if (len != 0.0) {                                                 
        phi = acos(aux.x/len) / (2 * M_PI);
        if (aux.y < 0.) phi = 1. - phi;
    }                                                                 
    
    return Vector2(sqrt(1. - aux.z * aux.z), phi);
    //             ^^^^^^^^^^^^^^^^^^^^^^^^ = sin(acos(aux.z))
}

Vector3 Hemisphere::samplePriorToPoweredZenith(double  elevationRand,
                                               double  azimuthRand,
                                               int     power,
                                               double* pPDFVal) const
{
    // here the elevation angle is the angle between direction and base plane!
    double azimuth = 2. * PI * azimuthRand;
    double sinElev = 1. - pow(elevationRand,  1./(power+1.));
    double cosElev = sqrt(1. - sinElev * sinElev);

    Vector3 dir;

    ESG_ANGLES_TO_VECTOR(azimuth, cosElev, sinElev, dir);

    if (pPDFVal) *pPDFVal = (power + 1.) * pow(cosElev, power) / (2. * PI);

    return dir;
}

Vector2 Hemisphere::samplePriorToPoweredZenith2D(double  elevationRand,
                                                 double  azimuthRand,
                                                 int     power,
                                                 double* pPDFVal) const
{
    double sinElev = 1. - pow(elevationRand,  1./(power+1.));
    double cosElev = sqrt(1. - sinElev * sinElev);
    if (pPDFVal) *pPDFVal = (power + 1.) * pow(cosElev, power) / (2. * PI);
    return Vector2(asin(sinElev), 2. * PI * azimuthRand);
}

Vector2 Hemisphere::dir2uvPriorToPoweredZenith(const Vector3& dir,
                                               int            power) const
{
    Vector3 aux(dir.x * _xAxis.x + dir.y * _xAxis.y + dir.z * _xAxis.z,  
                dir.x * _yAxis.x + dir.y * _yAxis.y + dir.z * _yAxis.z,  
                dir.x * _zAxis.x + dir.y * _zAxis.y + dir.z * _zAxis.z);

    double len = sqrt(aux.x * aux.x + aux.y * aux.y);                 
    double phi = 0.0;                                                 
    if (len != 0.0) {                                                 
        phi = acos(aux.x/len) / (2 * M_PI);
        if (aux.y < 0.) phi = 1. - phi;
    }                                                                 
    
    return Vector2(1. - pow(sqrt(1. - aux.z*aux.z), ((double)power+1.)), phi);
    //                      ^^^^^^^^^^^^^^^^^^^^^ = sin(acos(aux.z))
}

#undef ESG_ANGLES_TO_VECTOR

---