src/rcrt/primitives/DisplacedTriangle.cpp

Go to the documentation of this file.

#include "DisplacedTriangle.h"

using namespace std;

namespace rcrt
{

DisplacedTriangle::DisplacedTriangle(Vertex* av,Vertex* bv,Vertex* cv, SolidObject* parent) : Triangle(av,bv,cv,parent)
{
        da = Point3D(a->pos() + ((a->normal() / (a->normal() * faceNormal)) * parent->getMaxDispl()));
        db = Point3D(b->pos() + ((b->normal() / (b->normal() * faceNormal)) * parent->getMaxDispl()));
        dc = Point3D(c->pos() + ((c->normal() / (c->normal() * faceNormal)) * parent->getMaxDispl()));
                
        box.extend(da);
        box.extend(a->pos() + (faceNormal * (parent->getMaxDispl() + 100*numeric_limits<float>::epsilon())));
        box.extend(db);
        box.extend(b->pos() + (faceNormal * (parent->getMaxDispl() + 100*numeric_limits<float>::epsilon())));
        box.extend(dc);
        box.extend(c->pos() + (faceNormal * (parent->getMaxDispl() + 100*numeric_limits<float>::epsilon())));
        
        upper = new Triangle(new Vertex(da,a->normal()),new Vertex(db,b->normal()),new Vertex(dc,c->normal()),parent);
        
        ab = new ConvexQuad(parent,a->pos(),b->pos(),db,da);
        bc = new ConvexQuad(parent,b->pos(),c->pos(),dc,db);
        ca = new ConvexQuad(parent,c->pos(),a->pos(),da,dc);
        
        triangles = new vector<Triangle*>();
        
        int N = parent->getDisplDens();
        
        for (int i = 0; i < N; i++) {
                for (int j = 0; j < N; j++) {
                        for (int k = 0; k < N; k++) {
                                if (i+j+k == N-1) {
                                        std::cout << "lower" << std::endl;
                                        Point2D uv1 = Point2D(j/N,k/N);
                                        Point2D uv2 = Point2D((j+1)/N,k/N);
                                        Point2D uv3 = Point2D(j/N,(k+1)/N);
                                        
                                        Point3D po1 = getDispPoint(uv1);
                                        Point3D po2 = getDispPoint(uv2);
                                        Point3D po3 = getDispPoint(uv3);
                                        
                                        Vec3D sno = (po2 - po1).crossP(po3 - po1).normalize();
                                        Triangle* t = new Triangle(po1,po2,po3,sno,parent);
//                                      Triangle* t = new Triangle(new Vertex(po1,sno),new Vertex(po2,sno),new Vertex(po3,sno),parent, true);
                                        
                                        triangles->push_back(t);
                                }
                                
                                if (i+j+k == N-2) {
                                        std::cout << "upper" << std::endl;
                                        Point2D uv1 = Point2D(j/N,(k+1)/N);
                                        Point2D uv2 = Point2D((j+1)/N,k/N);
                                        Point2D uv3 = Point2D((j+1)/N,(k+1)/N);
                                        
                                        Point3D po1 = getDispPoint(uv1);
                                        Point3D po2 = getDispPoint(uv2);
                                        Point3D po3 = getDispPoint(uv3);
                                        
                                        Vec3D sno = (po2 - po1).crossP(po3 - po1).normalize();
                                        Triangle* t = new Triangle(po1,po2,po3,sno,parent);
//                                      Triangle* t = new Triangle(new Vertex(po1,sno),new Vertex(po2,sno),new Vertex(po3,sno),parent, true);
                                        
                                        triangles->push_back(t);
                                }
                        }               
                }
        }
        
        tree.setTraceables(triangles);
}

DisplacedTriangle::~DisplacedTriangle()
{
        if (ab) delete ab;
        if (bc) delete bc;
        if (ca) delete ca;
        
        if (upper) delete upper;
}

Intersection DisplacedTriangle::baseIntersect(Ray& ray) const
{
        ray.tris++;
        const Vec3D edge1 = b->pos()-a->pos();
    const Vec3D edge2 = c->pos()-a->pos();

    const Vec3D pvec = ray.dir().crossP(edge2);

    const float det = edge1 * pvec;
    if (fabs(det) < numeric_limits<float>::epsilon()) return Intersection();

    const float inv_det = 1.0f / det;

    const Vec3D tvec = ray.org()-a->pos();
    float lambda = tvec * pvec;
    lambda *= inv_det;

    if (lambda < 0.0f || lambda > 1.0f) return Intersection();

    const Vec3D qvec = tvec.crossP(edge1);
    float mue = ray.dir() * qvec;
    mue *= inv_det;

    if (mue < 0.0f || mue+lambda > 1.0f) return Intersection();

    float f = edge2 * qvec;
    f = f * inv_det - numeric_limits<float>::epsilon();
    if (ray.maxDist() <= f || f <  numeric_limits<float>::epsilon()  ) return Intersection();

    //Vec3D normal = lambda * b->normal() + mue * c->normal() + (1-lambda-mue) * a->normal();
    //normal.normalize();
    bool flip = ray.dir()*faceNormal >= numeric_limits<float>::epsilon();
    //std::cout << "[TRI-Intersect] t = " << f << " backside=" << flip << " distToOrg=" << (ray.org()-getCentroid()).norm() << " angle="<<ray.dir()*normal << std::endl;
    return Intersection(f, this, ray.atDistance(f), lambda, mue, flip);
}

Intersection DisplacedTriangle::intersect(Ray& ray) const
{
        
        //Intersection is = tree.intersect(ray);
        
        Intersection is;
        
        for(std::vector<Triangle*>::iterator i = triangles->begin(); i != triangles->end(); ++i) {
                is = (*i)->intersect(ray);
                if (is.isValid()) {
                        return is;
                }
        }
        
        if (is.isValid()) {
                is.setPrimitive(this);
                return is;
        }
        
        return Intersection();
        
        //intersect BBox
        if (!box.intersect(ray).isValid())
                return Intersection();
        
//      Intersection is = triangles.intersect(ray);
//      if (is.isValid()) {
//              is.setPrimitive(this);
//              return is;
//      }
        //std::cout << "intersecting volume" << std::endl;
        
        //const ImageTexture* const dispMap = parent->getDisplMap();
        float N = parent->getDisplDens();
        //const float maxDispl = parent->getMaxDispl();
        
        Intersection topIs, bottomIs, abIs, bcIs, caIs, nearIs, farIs;
        topIs = upper->intersect(ray);
        bottomIs = baseIntersect(ray);
//      bool nearFound = false;
//      bool farFound = false;
        
//      nearIs = topIs;
//      farIs = bottomIs;
        
        //std::cout << "topIS: " << topIs.isValid() << " bottomIS: " << bottomIs.isValid() << std::endl;
        
        if (false && (topIs.isValid() || bottomIs.isValid())) {
                //std::cout << "topIS: " << topIs.isValid() << " bottomIS: " << bottomIs.isValid() << std::endl;
                nearIs = topIs.isValid() ? topIs : bottomIs;
                Intersection is = nearIs;
                is.setPrimitive(this);
                is.setSNormalL(Vec3D(1,1,1));
                return is;
        }
        
        if (!(topIs.isValid() && bottomIs.isValid())) {
                if(topIs.isValid() || bottomIs.isValid()) {
                        nearIs = topIs.isValid() ? topIs : bottomIs;
                }
                
                abIs = ab->intersect(ray);
                bcIs = bc->intersect(ray);
                caIs = ca->intersect(ray);
                
                //std::cout << "abIS: " << abIs.isValid() << " bcIS: " << bcIs.isValid() << " caIS: " << caIs.isValid() << std::endl;
                
                if (false && (abIs.isValid() || bcIs.isValid() || caIs.isValid())) {
                        Intersection is = nearIs;
                        is.setPrimitive(this);
                        //is.setSNormal(Vec3D(1,1,1));
                        return is;
                }
                
                if (abIs.isValid() && abIs.getDistance() < nearIs.getDistance()) {
                        //std::cout << "setting ab to near" << std::endl;
                        farIs = nearIs;
                        nearIs = abIs;
                }
                else if (abIs.isValid() && abIs.getDistance() < farIs.getDistance()) {
                        //std::cout << "setting ab to far" << std::endl;
                        farIs = abIs;
                }
                if (bcIs.isValid() && bcIs.getDistance() < nearIs.getDistance()) {
                        //std::cout << "setting bc to near" << std::endl;
                        farIs = nearIs;
                        nearIs = bcIs;
                }
                else if (bcIs.isValid() && bcIs.getDistance() < farIs.getDistance()) {
                        //std::cout << "setting bc to far" << std::endl;
                        farIs = bcIs;
                }
                if (caIs.isValid() && caIs.getDistance() < nearIs.getDistance()) {
                        //std::cout << "setting ca to near" << std::endl;
                        farIs = nearIs;
                        nearIs = caIs;
                }
                else if (caIs.isValid() && caIs.getDistance() < farIs.getDistance()) {
                        //std::cout << "setting ca to far" << std::endl;
                        farIs = caIs;
                }
        }
        else {
                if(topIs.getDistance() < bottomIs.getDistance()) {
                        nearIs = topIs;
                        farIs = bottomIs;
                }
                else {
                        nearIs = bottomIs;
                        farIs = topIs;
                }
        }
        
        //std::cout << "nearIS: " << nearIs.isValid() << " farIS: " << farIs.isValid() << std::endl;
        
        if (!(nearIs.isValid() && farIs.isValid())) {
//              std::cout << "not both valid" << std::endl;
//              Intersection is = Intersection(1,this, Point3D(1,1,1),1,1,false);
//              //is.setPrimitive(this);
//              is.setSNormalL(Vec3D(1,1,1));
//              return is;
                return Intersection();
        }       
        
        if (false) {
                Intersection is = nearIs;
                is.setPrimitive(this);
                is.setSNormalL(Vec3D(1,1,1));
                return is;
        }
        
        //std::cout << "near and far valid" << std::endl;
        
        int i = 0,iStop = N,j = 0,jStop = N,k = 0,kStop = N;
        
        int lastChange = 0;
        
        Point3D subPA, subPB, subPC;
        Vec3D subNC;
        Point2D uva, uvb, uvc;
        bool rightOfC = false;
        
        // {iplus = 0, jminus = 1, kplus = 2, iminus = 3, jplus = 4,kminus = 5}
        
        // setting up start state
        
        // first intersection with bottom or top triangle
        
        //std::cout << "first is with top or bottom" << std::endl;
        
        if (nearIs.getPrimitive() == this || nearIs.getPrimitive() == upper) {
                
//              nearIs = topIs.isValid() ? topIs : bottomIs;
//              Intersection is = nearIs;
//              is.setPrimitive(this);
//              is.setSNormalL(Vec3D(1,1,1));
//              return is;
                
                const Point2D& bcoord = nearIs.getParams();
                float beta = bcoord[0];
                float gamma = bcoord[1];
                
                // get cell indices
                i = floor((1-beta-gamma) * N);
                j = floor(beta * N);
                k = floor(gamma * N);
                
                //std::cout << "i: " << i << " j: " << j << " k: " << k << std::endl;
                
                //first sub triangle is a lower triangle
                if (i+j+k == N-1) {
                        
                        //std::cout << "lower tri" << std::endl;
                        
//                      nearIs.setSNormalL(Vec3D(1,0,1));
//                      return nearIs;
                        
                        Point2D tuva = Point2D(j/N,k/N);
                        Point2D tuvb = Point2D((j+1)/N,k/N);
                        Point2D tuvc = Point2D(j/N,(k+1)/N);
                        
                        Point3D la = getPoint(tuva);
                        Point3D lb = getPoint(tuvb);
                        Point3D lc = getPoint(tuvc);
                        Point3D ua = upper->getPoint(tuva);
                        Point3D ub = upper->getPoint(tuvb);
                        Point3D uc = upper->getPoint(tuvc);
                        
                        //setting up points
                        
                        Point3D ta = getDispPoint(tuva);
                        Point3D tb = getDispPoint(tuvb);
                        Point3D tc = getDispPoint(tuvc);
                        
                        ConvexQuad testAB(parent,la,lb,ub,ua);
                        ConvexQuad testBC(parent,lb,lc,uc,ub);
                        ConvexQuad testCA(parent,lc,la,ua,uc);
                        
                        if (testAB.intersect(ray).isValid()) {
                                //std::cout << "testAB intersect" << std::endl;
                                //return nearIs;
                                lastChange = 0;
                                subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                        }
                        else if (testBC.intersect(ray).isValid()) {
                                //std::cout << "testBC intersect" << std::endl;
                                //return nearIs;
                                lastChange = 4;
                                subPA = tc; subPB = ta; subPC = tb; uva = tuvc; uvb = tuva; uvc = tuvb;
                        }
                        else if (testCA.intersect(ray).isValid()) {
                                //std::cout << "testCA intersect" << std::endl;
                                //return nearIs;
                                lastChange = 0;
                                subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                        }
                        else {
                                //std::cout << "no intersect" << std::endl;
                                lastChange = 0;
                                subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                                iStop = i; jStop = j; kStop = k;
                        }
                }
                
                //first sub triangle is a upper triangle
                if (i+k+j == N-2) {
                        
                        //std::cout << "upper tri" << std::endl;
                        
//                      nearIs.setSNormalL(Vec3D(0,1,0));
//                      return nearIs;
                        
                        Point2D tuva = Point2D(j/N,(k+1)/N);
                        Point2D tuvb = Point2D((j+1)/N,k/N);
                        Point2D tuvc = Point2D((j+1)/N,(k+1)/N);
                        
                        Point3D la = getPoint(tuva);
                        Point3D lb = getPoint(tuvb);
                        Point3D lc = getPoint(tuvc);
                        Point3D ua = upper->getPoint(tuva);
                        Point3D ub = upper->getPoint(tuvb);
                        Point3D uc = upper->getPoint(tuvc);
                        
                        Point3D ta = getDispPoint(tuva);
                        Point3D tb = getDispPoint(tuvb);
                        Point3D tc = getDispPoint(tuvc);
                        
                        //setting up points
                        
                        ConvexQuad testAB(parent,la,lb,ub,ua);
                        ConvexQuad testBC(parent,lb,lc,uc,ub);
                        ConvexQuad testCA(parent,lc,la,ua,uc);
                        
                        if (testAB.intersect(ray).isValid()) {
                                //std::cout << "testAB intersect" << std::endl;
                                //return nearIs;
                                lastChange = 1;
                                subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                        }
                        else if (testBC.intersect(ray).isValid()) {
                                //std::cout << "testBC intersect" << std::endl;
                                //return nearIs;
                                lastChange = 3;
                                subPA = ta; subPB = tb; subPC = tc; uva = tuva; uvb = tuvb; uvc = tuvc;
                        }
                        else if (testCA.intersect(ray).isValid()) {
                                //std::cout << "testCA intersect" << std::endl;
                                //return nearIs;
                                lastChange = 3;
                                subPA = ta; subPB = tb; subPC = tc; uva = tuva; uvb = tuvb; uvc = tuvc;
                        }
                        else {
                                //std::cout << "no intersect" << std::endl;
                                lastChange = 0;
                                subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                                iStop = i; jStop = j; kStop = k;
                        }
                }
        }
        // first Intersection with side quads
        else {
                
                //std::cout << "first IS with side" << std::endl;
                
                // determining cell to enter
                if(nearIs.getPrimitive() == ab) {
                        Point3D p = nearIs.getPosition();
                        Vec3D v = p - a->pos();
                        Vec3D orthInSide = (da - a->pos()).crossP(b->pos()-a->pos()).crossP(b->pos() - a->pos()).normalize();
                        orthInSide = orthInSide * (orthInSide * faceNormal);
                        Vec3D vp = orthInSide * (v * orthInSide);
                        //Vec3D vp = faceNormal * (v * faceNormal);
                        float f = fabs(vp.norm());
                        Point3D aAl = a->pos() + (da - a->pos()) * f;
                        Point3D bAl = b->pos() + (db - b->pos()) * f;
                        float x = (bAl - aAl).norm();
                        float y = (p - aAl).norm();
                        float ratio = y / x;
                        
                        //std::cout << "ratio: " << ratio << std::endl;
                        
                        // calculating indices - gamma=0
                        float alpha = 1 - ratio;
                        float beta = ratio;
                        float gamma = 0;
                        
                        // get cell indices
                        i = floor(alpha * N);
                        j = floor(beta * N);
                        k = floor(gamma * N);
                        
                        //setting up points
                        Point2D tuva = Point2D(j/N,k/N);
                        Point2D tuvb = Point2D((j+1)/N,k/N);
                        Point2D tuvc = Point2D(j/N,(k+1)/N);
                        
                        Point3D ta = getDispPoint(tuva);
                        Point3D tb = getDispPoint(tuvb);
                        Point3D tc = getDispPoint(tuvc);
                        
                        lastChange = 2;
                        subPA = ta; subPB = tb; subPC = tc; uva = tuva; uvb = tuvb; uvc = tuvc;
                }
                else if(nearIs.getPrimitive() == bc) {
                        Point3D p = nearIs.getPosition();
                        Vec3D v = p - b->pos();
                        Vec3D orthInSide = (db - b->pos()).crossP(c->pos()-b->pos()).crossP(c->pos() - b->pos()).normalize();
                        orthInSide = orthInSide * (orthInSide * faceNormal);
                        Vec3D vp = orthInSide * (v * orthInSide);
                        float f = fabs(vp.norm());
                        Point3D bAl = b->pos() + (db - b->pos()) * f;
                        Point3D cAl = c->pos() + (dc - c->pos()) * f;
                        float x = (cAl - bAl).norm();
                        float y = (p - bAl).norm();
                        float ratio = y / x;
                        
                        //std::cout << "ratio: " << ratio << std::endl;
                        
                        // calculating indices - gamma=0
                        float alpha = 0;
                        float beta = 1 - ratio;
                        float gamma = ratio;
                        
                        // get cell indices
                        i = floor(alpha * N);
                        j = floor(beta * N);
                        k = floor(gamma * N);
                        
                        //setting up points and uvs
                        Point2D tuva = Point2D(j/N,k/N);
                        Point2D tuvb = Point2D((j+1)/N,k/N);
                        Point2D tuvc = Point2D(j/N,(k+1)/N);
                        
                        Point3D ta = getDispPoint(tuva);
                        Point3D tb = getDispPoint(tuvb);
                        Point3D tc = getDispPoint(tuvc);
                        
                        lastChange = 0;
                        subPA = tb; subPB = tc; subPC = ta; uva = tuvb; uvb = tuvc; uvc = tuva;
                }
                else if(nearIs.getPrimitive() == ca) {
                        Point3D p = nearIs.getPosition();
                        Vec3D v = p - c->pos();
                        Vec3D orthInSide = (dc - c->pos()).crossP(a->pos()-c->pos()).crossP(a->pos() - c->pos()).normalize();
                        orthInSide = orthInSide * (orthInSide * faceNormal);
                        Vec3D vp = orthInSide * (v * orthInSide);
                        float f = fabs(vp.norm());
                        Point3D cAl = c->pos() + (dc - c->pos()) * f;
                        Point3D aAl = a->pos() + (da - a->pos()) * f;
                        float x = (aAl - cAl).norm();
                        float y = (p - cAl).norm();
                        float ratio = y / x;
                        
                        //std::cout << "ratio: " << ratio << std::endl;
                        
                        // calculating indices - gamma=0
                        float alpha = ratio;
                        float beta = 0;
                        float gamma = 1 - ratio;
                        
                        // get cell indices
                        i = floor(alpha * N);
                        j = floor(beta * N);
                        k = floor(gamma * N);
                        
                        //setting up points and uvs
                        Point2D tuva = Point2D(j/N,k/N);
                        Point2D tuvb = Point2D((j+1)/N,k/N);
                        Point2D tuvc = Point2D(j/N,(k+1)/N);
                        
                        Point3D ta = getDispPoint(tuva);
                        Point3D tb = getDispPoint(tuvb);
                        Point3D tc = getDispPoint(tuvc);
                        
                        lastChange = 2;
                        subPA = tc; subPB = ta; subPC = tb; uva = tuvc; uvb = tuva; uvc = tuvb;
                }
        }
        
        if (!(i+j+k == N-1 || i+j+k == N-2)) {
                std::cout << "omfg!!! i=" << i << " j=" << j << " k=" << k << " N=" << N << std::endl; 
        }
        
        if (false) {
                if (i+j+k == N-1) {
                        Intersection is = nearIs;
                        is.setPrimitive(this);
                        is.setSNormalL(Vec3D(1,0,0));
                        return is;
                }
                if (i+j+k == N-2) {
                        Intersection is = nearIs;
                        is.setPrimitive(this);
                        is.setSNormalL(Vec3D(0,1,0));
                        return is;
                }
        }
        
//      //determining end cells
//      
//      if (farIs.getPrimitive() == this || farIs.getPrimitive() == upper) {
//              // end intersection with top or bottom
//              const Point2D& bcoord = nearIs.getParams();
//              float beta = bcoord[0];
//              float gamma = bcoord[1];
//              
//              // get cell indices
//              iStop = floor((1-beta-gamma) * N);
//              jStop = floor(beta * N);
//              kStop = floor(gamma * N);
//      }
//      else {
//              // end intersection with side
//              if(farIs.getPrimitive() == ab) {
//                      Point3D p = farIs.getPosition();
//                      Vec3D v = p - a->pos();
//                      Vec3D orthInSide = (da - a->pos()).crossP(b->pos()-a->pos()).crossP(b->pos() - a->pos()).normalize();
//                      orthInSide = orthInSide * (orthInSide * faceNormal);
//                      Vec3D vp = orthInSide * (v * orthInSide);
//                      //Vec3D vp = faceNormal * (v * faceNormal);
//                      float f = fabs(vp.norm());
//                      Point3D aAl = a->pos() + (da - a->pos()) * f;
//                      Point3D bAl = b->pos() + (db - b->pos()) * f;
//                      float x = (bAl - aAl).norm();
//                      float y = (p - aAl).norm();
//                      float ratio = y / x;
//                      
//                      // calculating indices - gamma=0
//                      float alpha = 1 - ratio;
//                      float beta = ratio;
//                      float gamma = 0;
//                      
//                      // get cell indices
//                      iStop = floor(alpha * N);
//                      jStop = floor(beta * N);
//                      kStop = floor(gamma * N);
//              }
//              else if(farIs.getPrimitive() == bc) {
//                      Point3D p = farIs.getPosition();
//                      Vec3D v = p - b->pos();
//                      Vec3D orthInSide = (db - b->pos()).crossP(c->pos()-b->pos()).crossP(c->pos() - b->pos()).normalize();
//                      orthInSide = orthInSide * (orthInSide * faceNormal);
//                      Vec3D vp = orthInSide * (v * orthInSide);
//                      float f = fabs(vp.norm());
//                      Point3D bAl = b->pos() + (db - b->pos()) * f;
//                      Point3D cAl = c->pos() + (dc - c->pos()) * f;
//                      float x = (cAl - bAl).norm();
//                      float y = (p - bAl).norm();
//                      float ratio = y / x;
//                      
//                      // calculating indices - gamma=0
//                      float alpha = 0;
//                      float beta = 1 - ratio;
//                      float gamma = ratio;
//                      
//                      // get cell indices
//                      iStop = floor(alpha * N);
//                      jStop = floor(beta * N);
//                      kStop = floor(gamma * N);
//              }
//              else if(farIs.getPrimitive() == ca) {
//                      Point3D p = farIs.getPosition();
//                      Vec3D v = p - c->pos();
//                      Vec3D orthInSide = (dc - c->pos()).crossP(a->pos()-c->pos()).crossP(a->pos() - c->pos()).normalize();
//                      orthInSide = orthInSide * (orthInSide * faceNormal);
//                      Vec3D vp = orthInSide * (v * orthInSide);
//                      float f = fabs(vp.norm());
//                      Point3D cAl = c->pos() + (dc - c->pos()) * f;
//                      Point3D aAl = a->pos() + (da - a->pos()) * f;
//                      float x = (aAl - cAl).norm();
//                      float y = (p - cAl).norm();
//                      float ratio = y / x;
//                      
//                      // calculating indices - gamma=0
//                      float alpha = ratio;
//                      float beta = 0;
//                      float gamma = 1 - ratio;
//                      
//                      // get cell indices
//                      iStop = floor(alpha * N);
//                      jStop = floor(beta * N);
//                      kStop = floor(gamma * N);
//              }
//      }
        
        // brute Force
        if(true) {
                for (i = 0; i < N; i++) {
                        for (j = 0; j < N; j++) {
                                for (k = 0; k < N; k++) {
                                        if (i+j+k == N-1) {
                                                //std::cout << "lower" << std::endl;
                                                Point2D uv1 = Point2D(j/N,k/N);
                                                Point2D uv2 = Point2D((j+1)/N,k/N);
                                                Point2D uv3 = Point2D(j/N,(k+1)/N);
                                                
                                                Point3D po1 = getDispPoint(uv1);
                                                Point3D po2 = getDispPoint(uv2);
                                                Point3D po3 = getDispPoint(uv3);
                                                
                                                Vec3D sno = (po2 - po1).crossP(po3 - po1).normalize();
                                                
                                                Triangle t = Triangle(po1,po2,po3,sno,parent);
                                                
                                                Intersection is = t.intersect(ray);
                                                
                                                if(is.isValid()) {
                                                        is.setPrimitive(this);
                                                        is.setGNormalL(sno);
                                                        is.setSNormalL(sno);
                                                        return is;
                                                }
                                        }
                                        
                                        if (i+j+k == N-2) {
                                                //std::cout << "upper" << std::endl;
                                                Point2D uv1 = Point2D(j/N,(k+1)/N);
                                                Point2D uv2 = Point2D((j+1)/N,k/N);
                                                Point2D uv3 = Point2D((j+1)/N,(k+1)/N);
                                                
                                                Point3D po1 = getDispPoint(uv1);
                                                Point3D po2 = getDispPoint(uv2);
                                                Point3D po3 = getDispPoint(uv3);
                                                
                                                Vec3D sno = (po2 - po1).crossP(po3 - po1).normalize();
                                                
                                                Triangle t = Triangle(po1,po2,po3,sno,parent);

                                                Intersection is = t.intersect(ray);

                                                if(is.isValid()) {
                                                        is.setPrimitive(this);
                                                        is.setGNormalL(sno);
                                                        is.setSNormalL(sno);
                                                        return is;
                                                }
                                        }
                                }               
                        }
                }
                return Intersection();
        }
        
        //std::cout << "Idices: i=" << i << " j=" << j << " k=" << k << std::endl;
        
        // calculating C Normal
        
//      std::cout << "uvc: beta=" << uvc[0] << " gamma=" << uvc[1] << std::endl;
//      std::cout << "uva: beta=" << uva[0] << " gamma=" << uva[1] << std::endl;
//      std::cout << "uvb: beta=" << uvb[0] << " gamma=" << uvb[1] << std::endl;
        
        subNC = getDispNormal(uvc);
        
        // traversal
        
        float delta = 1 / N;
        
        //return Intersection();
        
        while(true) {
                //std::cout << "loop" << std::endl;
                //return Intersection();
                //std::cout << "Volume Traversal: change:" << lastChange << "  Idices: i=" << i << " j=" << j << " k=" << k << std::endl;
                Vec3D subNormal = (subPB - subPA).crossP(subPC - subPA).normalize();
                Triangle subTri(subPA, subPB, subPC, subNormal, parent);
                //std::cout << "Normal: " << subNormal << "  A: " << subPA << "  B: " << subPB << "  C: " << subPC << std::endl;
                Intersection subIs = subTri.intersect(ray);
                if (subIs.isValid()) {
//                      if (i+j+k == N - 1)
//                              subNormal = Vec3D(1,0,1);
//                      if (i+j+k == N - 2)
//                              subNormal = Vec3D(0,1,0);
                        //std::cout << "valid subIS" << std::endl;
                        subIs.setPrimitive(this);
                        subIs.setGNormalL(subNormal);
                        //std::cout << "setting sNormal " << subNormal << std::endl;
                        subIs.setSNormalL(subNormal);
                        //subIs.setUV(getUV(subIs.getParams()));
                        //std::cout << "normal set to " << subIs.getSNormal() << " and returning"<< std::endl;
                        return subIs;
                }
                //std::cout << "no valid subIS" << std::endl;
                if (i == iStop && j == jStop && k == kStop) {
                        //std::cout << "Stop Coordinates reached" << std::endl;
                        return Intersection();
                }
                rightOfC = ((subNC.crossP(ray.org() - subPC)) * ray.dir() < 0);
                if(rightOfC) {
                        subPA = subPC; uva = uvc;
                }
                else {
                        subPB = subPC; uvb = uvc;
                }
                // Take advantage of numbering. 5 = −1 mod 6
                lastChange = (lastChange + (rightOfC ? 1 : 5)) % 6;
                //std::cout << "change: " << lastChange << std::endl;
                if(lastChange == 3 /* change == iminus */) {
                        if (--i < 0){
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D((j+1)*delta, (k+1)*delta);
                }
                else if(lastChange == 0 /* change == iplus */) {
                        if (++i >= N ) {
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D( j*delta, k*delta);
                }
                else if(lastChange == 1 /* change == jminus */) {
                        if (--j < 0) {
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D( j*delta, (k+1)*delta);
                }
                else if(lastChange == 4 /* change == jplus */) {
                        if (++j >= N ) {
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D( (j+1)*delta, k*delta);
                }
                else if(lastChange == 5 /* change == kminus */) {
                        if (--k < 0) {
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D( (j+1)*delta, k*delta);
                }
                else if(lastChange == 2 /* change == kplus */) {
                        if (++k >= N ) {
                                //std::cout << "moving out of volume" << std::endl;
                                return Intersection();
                        }
                        uvc = Point2D( j*delta, (k+1)*delta);
                }
                subPC = getDispPoint(uvc);
                subNC = getDispNormal(uvc);
        }

}

Point3D DisplacedTriangle::getDispPoint(float beta, float gamma) const
{
        Point3D lp = getPoint(beta,gamma);
        return lp + ((upper->getPoint(beta,gamma)-lp).normalize()*parent->getMaxDispl()*parent->getDisplMap()->getColor(getUV(beta,gamma)).getLuminance());
}

Vec3D DisplacedTriangle::getDispNormal(float beta, float gamma) const
{
        Point3D lp = getPoint(beta,gamma);
        return (upper->getPoint(beta,gamma)-lp).normalize();
}

Point3D DisplacedTriangle::getDispPoint(const Point2D& p) const
{
        float beta = p[0];
        float gamma = p[1];
        return getDispPoint(beta,gamma);
}

Vec3D DisplacedTriangle::getDispNormal(const Point2D& p) const
{
        float beta = p[0];
        float gamma = p[1];
        return getDispNormal(beta,gamma);
}

}

---