FDH14.cc

Go to the documentation of this file.

#include <esg/geometry/FDH14.h>
#include <esg/mesh/FDH14Mesh.h>

using namespace esg;

const double   FDH14::SCALE         = .5773502691896; // 1/sqrt(3)
const unsigned FDH14::DIRS          =  7;             // number of directions
const unsigned FDH14::SIZE          =  14;            // 2 * DIRS
const float    FDH14::FDHMat[14][3] =                 // matrix of FDH
{
    {      1,      0,      0 },
    {      0,      1,      0 },
    {      0,      0,      1 },
    {  SCALE,  SCALE,  SCALE },
    {  SCALE, -SCALE,  SCALE },
    { -SCALE,  SCALE,  SCALE },
    {  SCALE,  SCALE, -SCALE },
    {     -1,      0,      0 },
    {      0,     -1,      0 },
    {      0,      0,     -1 },
    { -SCALE, -SCALE, -SCALE },
    { -SCALE,  SCALE, -SCALE },
    {  SCALE, -SCALE, -SCALE },
    { -SCALE, -SCALE,  SCALE }
};


void FDH14::_rotateX(float a)
{
    Matrix4 trMat; trMat.rotX(a);
    RotInfo r(trMat);
    float*  newVal = new float [2*_dirs];
        
    for (register unsigned i = 0; i < 2*_dirs; i++)
        newVal[i] = RotInfo::fast_lp(r.get(i), _values);

    delete [] _values;
    _values = newVal;
}

void FDH14::_rotateY(float a)
{
    Matrix4 trMat; trMat.rotY(a);
    RotInfo r(trMat);
    float*  newVal = new float [2*_dirs];
        
    for (register unsigned i = 0; i < 2*_dirs; i++)
        newVal[i] = RotInfo::fast_lp(r.get(i), _values);

    delete [] _values;
    _values = newVal;
}

void FDH14::_rotateZ(float a)
{
    Matrix4 trMat; trMat.rotZ(a);
    RotInfo r(trMat);
    float*  newVal = new float [2*_dirs];
        
    for (register unsigned i = 0; i < 2*_dirs; i++)
        newVal[i] = RotInfo::fast_lp(r.get(i), _values);

    delete [] _values;
    _values = newVal;
}

void FDH14::_rotate(float a, const Vector3& axis)
{
    Matrix4 trMat; trMat.rotationGL(a, axis);
    RotInfo r(trMat);
    float*  newVal = new float [2*_dirs];
        
    for (register unsigned i = 0; i < 2*_dirs; i++)
        newVal[i] = RotInfo::fast_lp(r.get(i), _values);

    delete [] _values;
    _values = newVal;
}

void FDH14::_rotate(const Matrix3& rotMat)
{
    Matrix4 trMat(rotMat);
    RotInfo r(trMat);
    float*  newVal = new float [2*_dirs];
        
    for (register unsigned i = 0; i < 2*_dirs; i++)
        newVal[i] = RotInfo::fast_lp(r.get(i), _values);

    delete [] _values;
    _values = newVal;
}

void FDH14::_transform(const Matrix4& trMat)
{
    Vector3 trVec, scVec;
    Matrix3 rotMat;

    trMat.get(rotMat, trVec, scVec);

    if ((scVec.x == scVec.y) && (scVec.y == scVec.z)) {
        RotInfo  r(trMat);
        float*   newVal = new float [2*_dirs];
        
        for (register unsigned i = 0; i < 2*_dirs; i++)
            newVal[i] = RotInfo::fast_lp(r.get(i), _values);

        delete [] _values;
        _values = newVal;
    } else 
        FDH::_transform(trMat);
}

Mesh* FDH14::_mesh(int d) const
{
    float ar, br, cr;

    ar = _values[0] + _values[7];
    br = _values[1] + _values[8];
    cr = _values[2] + _values[9];
    
    return (Mesh*) new FDH14Mesh(_mat, _values, ar, br, cr);
}

void FDH14::_duplicate_attributes(const Geometry& src)
{
    FDH::_duplicate_attributes(src);
}



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

#define BVH_INIT_CACHE  \
   if (!pPE->pCache) { \
       pPE->pCache = new FDHCache(4); \
       co = &(((FDHCache*)pPE->pCache)->oval[0]); \
       cd = &(((FDHCache*)pPE->pCache)->dval[0]); \
       *(co+0) = SCALE * ( origin.x + origin.y + origin.z); \
       *(co+1) = SCALE * ( origin.x - origin.y + origin.z); \
       *(co+2) = SCALE * (-origin.x + origin.y + origin.z); \
       *(co+3) = SCALE * ( origin.x + origin.y - origin.z); \
       *(cd+0) = SCALE * ( direction.x + direction.y + direction.z); \
       *(cd+1) = SCALE * ( direction.x - direction.y + direction.z); \
       *(cd+2) = SCALE * (-direction.x + direction.y + direction.z); \
       *(cd+3) = SCALE * ( direction.x + direction.y - direction.z); \
   } else { \
       co = &(((FDHCache*)pPE->pCache)->oval[0]); \
       cd = &(((FDHCache*)pPE->pCache)->dval[0]); \
   }

void FDH14::rayIntersection(PointEnv*      pPE,
                            int            mask,
                            const Vector3& origin,
                            const Vector3& direction,
                            float          maxDist)
{
    if (!pPE) {
        std::cerr << "FDH14::rayIntersection(): ";
        std::cerr << "No point environment structure defined" << std::endl;
        return;
    }

    pPE->mask = ENV_HAVE_NOTHING;


    register double t_in, t_out, t1, t2;

    if (mask&ENV_WANT_UV_COORD)
        fprintf(stderr,"FDH14::rayIntersection(): UV mapping is not implemented");

    if (mask & ENV_WANT_NORMAL) {
        /*
         * Use a little slower macros
         */
        int plane1 = 0, plane2 = 0;
        FDH_CLIP_LINE_INDEX_INIT(origin.x, direction.x, 0, 7);
        FDH_CLIP_LINE_INDEX_CONT(origin.y, direction.y, 1, 8);
        FDH_CLIP_LINE_INDEX_CONT(origin.z, direction.z, 2, 9);
        if (mask & ENV_USE_CACHE) {
            register float * co, * cd;
            BVH_INIT_CACHE;
            FDH_CLIP_LINE_INDEX_CONT(*co, *cd, 3, 10); co++; cd++;
            FDH_CLIP_LINE_INDEX_CONT(*co, *cd, 4, 11); co++; cd++;
            FDH_CLIP_LINE_INDEX_CONT(*co, *cd, 5, 12); co++; cd++;
            FDH_CLIP_LINE_INDEX_CONT(*co, *cd, 6, 13);
        } else {
            register double np0 = SCALE*(origin.x   +origin.y   +origin.z);
            register double np1 = SCALE*(direction.x+direction.y+direction.z);
            FDH_CLIP_LINE_INDEX_CONT(np0, np1, 3, 10);
            np0 = SCALE * (origin.x    - origin.y    + origin.z);
            np1 = SCALE * (direction.x - direction.y + direction.z);
            FDH_CLIP_LINE_INDEX_CONT(np0, np1, 4, 11);
            np0 = SCALE * (-origin.x    + origin.y    + origin.z);
            np1 = SCALE * (-direction.x + direction.y + direction.z);
            FDH_CLIP_LINE_INDEX_CONT(np0, np1, 5, 12);
            np0 = SCALE * (origin.x    + origin.y    - origin.z);
            np1 = SCALE * (direction.x + direction.y - direction.z);
            FDH_CLIP_LINE_INDEX_CONT(np0, np1, 6, 13);
        }
        if (t_in < Geometry::EPS)    // ray starts inside FDH
            pPE->normal.set(_mat[plane2][0],_mat[plane2][1],_mat[plane2][2]);
        else
            pPE->normal.set(_mat[plane1][0],_mat[plane1][1],_mat[plane1][2]);
        pPE->normalOrientation = PointEnv::OUTWARDS_NORMAL;
        pPE->mask |= ENV_HAVE_NORMAL;
    } else {
        /*
         * Use a little bit more efficient macros.
         */
        FDH_CLIP_LINE_INIT(origin.x, direction.x, 0, 7);
        FDH_CLIP_LINE_CONT(origin.y, direction.y, 1, 8);
        FDH_CLIP_LINE_CONT(origin.z, direction.z, 2, 9);
        if (mask & ENV_USE_CACHE) {
            register float * co, * cd;
            BVH_INIT_CACHE;
            FDH_CLIP_LINE_CONT(*co, *cd, 3, 10); co++; cd++;
            FDH_CLIP_LINE_CONT(*co, *cd, 4, 11); co++; cd++;
            FDH_CLIP_LINE_CONT(*co, *cd, 5, 12); co++; cd++;
            FDH_CLIP_LINE_CONT(*co, *cd, 6, 13);
        } else {
            register double np0 = SCALE*(origin.x   +origin.y   +origin.z);
            register double np1 = SCALE*(direction.x+direction.y+direction.z);
            FDH_CLIP_LINE_CONT(np0, np1, 3, 10);
            np0 = SCALE * (origin.x    - origin.y    + origin.z);
            np1 = SCALE * (direction.x - direction.y + direction.z);
            FDH_CLIP_LINE_CONT(np0, np1, 4, 11);
            np0 = SCALE * (-origin.x    + origin.y    + origin.z);
            np1 = SCALE * (-direction.x + direction.y + direction.z);
            FDH_CLIP_LINE_CONT(np0, np1, 5, 12);
            np0 = SCALE * (origin.x    + origin.y    - origin.z);
            np1 = SCALE * (direction.x + direction.y - direction.z);
            FDH_CLIP_LINE_CONT(np0, np1, 6, 13);
        } 
    }
    
    pPE->mask |= ENV_HAVE_INTERFERENCE;
    
    if (mask & ENV_WANT_INTERSECTION) {
        pPE->intersection.set(direction);
        if (t_in < Geometry::EPS)    // ray starts inside FDH
            pPE->intersection.scaleAdd(t_out, origin);
        else
            pPE->intersection.scaleAdd(t_in, origin);
        pPE->mask |= ENV_HAVE_INTERSECTION;
    }
    
    if (mask & ENV_WANT_DISTANCE) {
        if (t_in < Geometry::EPS) pPE->distance = t_out;
        else                      pPE->distance = t_in;
        pPE->mask |= ENV_HAVE_DISTANCE;
    }
}

#undef BVH_INIT_CACHE

Interval FDH14::extent(const Vector3& direction) const
{
    Interval interval;
    int      counter = 0;
    Vector3  dir(direction);
    
    if (direction.x < 0.) counter++;
    if (direction.y < 0.) counter++;
    if (direction.z < 0.) counter++;
    if (counter > 1) dir.negate(); // lower direction
    
    RotInfo r(dir);

    interval.max =   RotInfo::fast_lp(r.get(0), _values);
    interval.min = - RotInfo::fast_lp(r.get(1), _values);
    
    return interval;
}

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





//                        subclass FDH_14::RotInfo                           //

const double FDH14::RotInfo::ROT_EPS        = .0;
const double FDH14::RotInfo::ROT_SCALE      = .8660254;
const int    FDH14::RotInfo::_SgnCtrl[8][4] =
{
    /* --- */ {10, 12, 11, 13},
    /* +-- */ {12, 10,  6,  4},
    /* -+- */ {11,  6, 10,  5},
    /* ++- */ { 6, 11, 12,  3},
    /* --+ */ {13,  4,  5, 10},
    /* +-+ */ { 4, 13,  3, 12},
    /* -++ */ { 5,  3, 13, 11},
    /* +++ */ { 3,  5,  4,  6}
};

#define SWAPF(x,y) { float ftmp; ftmp=x; x=y; y=ftmp; }
void FDH14::RotInfo::fdh_precomp(FDH_Precomp *r, const float c[])
{
    float d[3];
    int  pt[3];
    register int i, j;
    
    j = 0;
    for (i = 0; i < 3; i++) {
        if (c[i] < 0.) {
            pt[i] = i + 7;
            d[i] = -c[i];
        } else {
            j = j | (1 << i);
            pt[i] = i;
            d[i] = c[i];
        }
    }
    
    r->p[2] = _SgnCtrl[j][0];
    
    if (d[0] <= d[1]) {
        if (d[1] <= d[2]) {
            /* d[0]<=d[1]<=d[2] */
            r->p[0] = pt[1];
            r->p[1] = pt[2];
            r->p[3] = _SgnCtrl[j][0 + 1];
        } else if (d[0] <= d[2]) {
            /* d[0]<=d[2]< d[1] */
            r->p[0] = pt[2];
            r->p[1] = pt[1];
            r->p[3] = _SgnCtrl[j][0 + 1];
            SWAPF(d[1], d[2]);
        } else {
            /* d[2]< d[0]<=d[1] */
            r->p[0] = pt[0];
            r->p[1] = pt[1];
            r->p[3] = _SgnCtrl[j][2 + 1];
            SWAPF(d[0], d[2]);
            SWAPF(d[1], d[2]);
        }
    } else {                      /* d[1]< d[0] */
        if (d[0] <= d[2]) {
            /* d[1]< d[0]<=d[2] */
            r->p[0] = pt[0];
            r->p[1] = pt[2];
            r->p[3] = _SgnCtrl[j][1 + 1];
            SWAPF(d[0], d[1]);
        } else if (d[1] <= d[2]) {
            /* d[1]<=d[2]< d[0] */
            r->p[0] = pt[2];
            r->p[1] = pt[0];
            r->p[3] = _SgnCtrl[j][1 + 1];
            SWAPF(d[0], d[1]);
            SWAPF(d[1], d[2]);
        } else {
            /* d[2]< d[1]< d[0] */
            r->p[0] = pt[1];
            r->p[1] = pt[0];
            r->p[3] = _SgnCtrl[j][2 + 1];
            SWAPF(d[0], d[2]);
        }
    }
    
    r->v1[0] = d[1] - d[0];
    r->v1[1] = d[2] - d[0];
    r->v1[2] = d[0] / SCALE;
    
    r->v2[0] = d[2] - d[1];
    r->v2[1] = (d[0] + d[1]) * ROT_SCALE;
    r->v2[2] = (d[1] - d[0]) * ROT_SCALE;
}
#undef SWAPF



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

FDH14::RotInfo::RotInfo()
{
    float trMat [4][4] = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0}, {0,0,0,1}};
  
    rotation = 0;
    r = NULL;
    recompute(trMat);
}

FDH14::RotInfo::RotInfo(const float trMat[4][4])
{
    rotation = 0;
    r = NULL;
    recompute(trMat);
}

FDH14::RotInfo::RotInfo(const Matrix4f& trMat)
{
    rotation = 0;
    r = NULL;
    recompute(trMat);
}

FDH14::RotInfo::RotInfo(const Vector3& direction)
{
    rotation = 1;
    r = NULL;
    recompute(direction);
}

void FDH14::RotInfo::recompute(const float trMat[4][4])
{
    float c[3];
    
    if (rotation == 1) {
        for (register int i = 0; i < 4; i++) // Check if values are actual
            for (register int j = 0; j < 4; j++)
                if (mat[i][j] != trMat[i][j]) goto NOT_ACTUAL;
        return;
    }
    
 NOT_ACTUAL:
    if (r) delete [] r;
    if ((r = new FDH_Precomp [SIZE]) == NULL) return;
    
    /* Precomupte array r */
    for (register unsigned i = 0; i < SIZE; i++) {
        for (register int j = 0; j < 3; j++) {
#if 0 // probably bug - reverse rotation
            c[j] = FDHMat[i][0] * trMat[j][0] +
                FDHMat[i][1] * trMat[j][1] +
                FDHMat[i][2] * trMat[j][2];
#endif
            c[j] = FDHMat[i][0] * trMat[0][j] +
                FDHMat[i][1] * trMat[1][j] +
                FDHMat[i][2] * trMat[2][j];
        }
        fdh_precomp(&r[i], c);
        r[i].tau = FDHMat[i][0] * trMat[0][3] +
            FDHMat[i][1] * trMat[1][3] +
            FDHMat[i][2] * trMat[2][3];
    }
    
    for (register int i = 0; i < 4; i++)
        for (register int j = 0; j < 4; j++)
            mat[i][j] = trMat[i][j];
    rotation = 1;
}

void FDH14::RotInfo::recompute(const Matrix4f& trMat)
{
    float c[3];
    
    if (rotation == 1) {
        for (register int i = 0; i < 4; i++) // Check if values are actual
            for (register int j = 0; j < 4; j++)
                if (mat[i][j] != trMat(i,j)) goto NOT_ACTUAL;
        return;
    }
    
 NOT_ACTUAL:
    if (r) delete [] r;
    if ((r = new FDH_Precomp [SIZE]) == NULL) return;
    
    /* Precomupte array r */
    for (register unsigned i = 0; i < SIZE; i++) {
        for (register int j = 0; j < 3; j++) {
#if 0 // probably bug - reverse rotation
            c[j] = FDHMat[i][0] * trMat(j,0) +
                FDHMat[i][1] * trMat(j,1) +
                FDHMat[i][2] * trMat(j,2);
#endif
            c[j] = (FDHMat[i][0] * trMat(0,j) +
                    FDHMat[i][1] * trMat(1,j) +
                    FDHMat[i][2] * trMat(2,j));
        }
        fdh_precomp(&r[i], c);
        r[i].tau = (FDHMat[i][0] * trMat(0,3) +
                    FDHMat[i][1] * trMat(1,3) +
                    FDHMat[i][2] * trMat(2,3));
    }
    
    for (register int i = 0; i < 4; i++)
        for (register int j = 0; j < 4; j++)
            mat[i][j] = trMat(i,j);
    rotation = 1;
}


void FDH14::RotInfo::recompute(const Vector3& direction)
{
    float backdir[3], dir[3];
    
    dir[0] = (float) direction.x;
    dir[1] = (float) direction.y;
    dir[2] = (float) direction.z;
    
    for (register int i = 0; i < 3; i++) backdir[i] = -dir[i];
    
    if (rotation == 0) { 
        // It is not necessary to compute values if there was not any rotation.
        // We can use stored values instead.
        if (mat[0][3] == direction.x &&
            mat[1][3] == direction.y &&
            mat[2][3] == direction.z) return;
        if (mat[0][3] == -direction.x &&
            mat[1][3] == -direction.y &&
            mat[2][3] == -direction.z) return;
    }
    
    delete [] r;
    if ((r = new FDH_Precomp [2]) == NULL) return;
    
    fdh_precomp(&r[0], dir);
    fdh_precomp(&r[1], backdir);
    r[0].tau = 0;
    r[1].tau = 0;
    
    mat[0][3] = direction.x;
    mat[1][3] = direction.y;
    mat[2][3] = direction.z;
    rotation = 0;
}

FDH14::RotInfo* FDH14::RotInfo::clone() const
{
    FDH14::RotInfo* pRet = new RotInfo();
    unsigned        rSize;
    
    for (register unsigned i = 0; i < 4; i++)
        for (register unsigned j = 0; j < 4; j++)
            pRet->mat[i][j] = mat[i][j];
    
    pRet->rotation = rotation;
    
    if (rotation) rSize = SIZE;
    else rSize = 2;
    
    if ((pRet->r = new FDH_Precomp [ rSize ]) == NULL) {
        delete pRet;
        return NULL;
    }
    
    for (register unsigned i = 0; i < rSize; i++) pRet->r[i] = r[i];
    
    return pRet;
}

Matrix4 FDH14::RotInfo::getMat() const
{
    Matrix4 retMat;
    for (register unsigned i = 0; i < 4; i++)
        for (register unsigned j = 0; j < 4; j++)
            retMat.setElement(i,j, mat[i][j]);
    return retMat;
}

---