Menger Sponge Intersection

#include <dgl/dgl.hpp> //my own library for window management etc etc etc - not important#include <dgl/dgli.hpp>#include "menger.hpp"void cb_main();void cb_mdown(int but);/*************************************************************************************************************************/Menger::Sponge sponge;int main(){    //init window    dgl::WINDOW_STYLE style;    style.width = 800;    style.height = 600;    style.xPos = 50;    style.yPos = 50;    style.name = "Menger Cube";    style.fullscreen = false;    style.maximise = false;    style.minimise = true;    style.resize = false;    dgl::init(style);    dgl::mainFunc(cb_main);    dgl::mouseDownFunc(cb_mdown);    /*********************/    //opengl settings    glClearColor(0,0,0,0);    dgl::proj2D(style.width,style.height,true); //sets up a 2d orthographic projection, true sets it to have y=0 at bottom of screen    /*********************/    //init menger cube instance    sponge.iterations = 7; //number of iterations    sponge.p = vec3(0,0,0); //position    sponge.A = vec33::scale(200,200,200); //transformation    sponge.iA = sponge.A.inverse();    sponge.itA = sponge.iA.transpose();    /*********************/    //enter main loop of program    dgl::mainloop();    return 0;}//captures right clicks and saves a screenshot to sshot.jpgvoid cb_mdown(int but){    if(but==dgl::MOUSE_RIGHT)    {        dgli::writeJpeg("sshot.jpg",90,dgli::capture());    }}/*************************************************************************************************************************///return pixel colour from ray cast through screen coordinate (sx,sy)vec3 cast(real sx, real sy){    //build ray    Menger::Ray cray;    cray.p = vec3(0,0,0); //camera is at origin    cray.v = vec3(sx-400.0,sy-300.0,400.0).normal(); //camera is facing foward along z-axis, no rotation to this    //transform into menger cube coordinate system    Menger::Ray tray = sponge.transform(cray);    //find intersection with time value + normal should it intersect    real t;    vec3 n;    if(Menger::intersectSponge(sponge,tray,t,n))    {        real c = 1-(t/400.0-0.25);        //transform normal back into world space, shift to (0,0,0) to (1,1,1) ranges, scale by time value)        return (sponge.normal(n).normal()*0.5 + vec3(0.5,0.5,0.5))*c;    }    else //background colour        return vec3(0,0,0);}/*************************************************************************************************************************///used in multisamplingconst real f13 = 1.0/3.0;//used in rotation of cubereal r = 0.0;void cb_main(){    /***************************************** render ********************************************/    //rendered as points    glBegin(GL_POINTS);    //skip over 5x5 blocks of pixels    for(int ay = 0; ay<600; ay+=5)    {        for(int ax = 0; ax<800; ax+=5)        {            //iterate through block's pixels            const int y2 = ay+5;            for(int y = ay; y<y2; y++)            {                const int x2 = ax+5;                for(int x = ax; x<x2; x++)                {                    real rx = x;                    real ry = y;                    /*16x multisampling*/                    //vec3 col = cast(rx,ry)*0.25 + (cast(rx-f13,ry)+cast(rx,ry-f13)+cast(rx+f13,ry)+cast(rx,ry+f13))*0.125                    //         +(cast(rx-f13,ry-f13)+cast(rx-f13,ry+f13)+cast(rx+f13,ry+f13)+cast(rx+f13,ry-f13))*0.0625;                    /*4x multisampling*/                    //vec3 col = (cast(rx-f13,ry-f13)+cast(rx-f13,ry+f13)+cast(rx+f13,ry+f13)+cast(rx+f13,ry-f13))*0.25;                    /*single ray/pixel*/                    vec3 col = cast(rx,ry);                    glColor3f(col.x,col.y,col.z);                    glVertex2i(x,y);                }            }        }    }    glEnd();    dgl::swapBuffers();    /*********************************** rotate menger sponge  ****************************************/    r += 0.06;    sponge.A*=vec33::rot_y(0.12*cos(r)*sin(3*r))*vec33::rot_x(0.12*sin(r)*sin(r*r))*vec33::rot_z(0.12*cos(2*r)*sin(r));    sponge.iA = sponge.A.inverse();    sponge.itA = sponge.iA.transpose();}

#ifndef menger_hpp#define menger_hpp#include <floats.hpp> //my own library for vectors+matrices of floats - not important//typedef incase i want to swap and use doubles instead of floatstypedef float real;typedef float3 vec3;typedef float33 vec33;/*************************************************************************************/namespace Menger{    //Ray class    class Ray    {        public:            vec3 p; //origin            vec3 v; //direction    };    /*************************************************************************************/    //Menger Sponge class    class Sponge    {        public:            int iterations; //number of subdivisions            vec3 p; //position            vec33 A, iA, itA; //transformation + inverse and inverse transpose                              //note that with a unit transformation, the sponge will be contained in a unit-cube            Ray transform(const Ray & r) const; //transforms ray to sponge's coordinate system            vec3 normal(const vec3 & n) const; //transforms normal to world coordinate system    };    /*************************************************************************************/    //used in ray traversal, 3 way min-max    inline real Min(real x, real y, real z)    {        real t;        return ((t=((x<y)?x:y))<z)?t:z;    }    inline real Max(real x, real y, real z)    {        real t;        return ((t=((x>y)?x:y))>z)?t:z;    }    /*************************************************************************************/    //given t values, find the index of the first node in the cubic-octree (my name for it)    int find_firstNode(real tx0, real ty0, real tz0,                       real tx1, real tx2,                       real ty1, real ty2,                       real tz1, real tz2);    //given t values, find the index of the next node in the tree (either x,y or z)    int find_nextNode(real tx3, real ty3, real tz3,                      int x, int y, int z);    //process a node of the cubic-octree    bool proc_subtree(real tx0, real ty0, real tz0,                      real tx3, real ty3, real tz3,                      int cite);    /*************************************************************************************/    //intersect a ray with a menger sponge (ray in sponge coordinates)    //out: t, n    //time value for ray, and normal at intersection (in sponge coordinates)    bool intersectSponge(const Sponge & o, Ray r, real & t, vec3 & n);};#endif

#include "menger.hpp"/*************************************************************************************/Menger::Ray Menger::Sponge::transform(const Menger::Ray & r) const{    Ray t;    t.p = iA*(r.p-p);    t.v = iA*r.v;    return t;}vec3 Menger::Sponge::normal(const vec3 & n) const{    return itA*n;}/*************************************************************************************///used so that i don't have to pass them into all of the methods that make use of themstatic real tval;static vec3 norm;static bool flipx;static bool flipy;static bool flipz;/*************************************************************************************//*note on cell indices: they are ordered according the formulaindex = x+(y*3)+(z*9)*/int Menger::find_firstNode(real tx0, real ty0, real tz0,                           real tx1, real tx2,                           real ty1, real ty2,                           real tz1, real tz2){    real max = Menger::Max(tx0,ty0,tz0);    int index = 0;    if(tx0==max)    {//enters cell in yz plane        if(ty2<tx0)      index += 6;        else if(ty1<tx0) index += 3;        if(tz2<tx0)      index += 18;        else if(tz1<tx0) index += 9;    }else if(ty0==max)    {//xz        if(tx2<ty0)      index += 2;        else if(tx1<ty0) index ++;        if(tz2<ty0)      index += 18;        else if(tz1<ty0) index += 9;    }else    {//xy        if(tx2<tz0)      index += 2;        else if(tx1<tz0) index ++;        if(ty2<tz0)      index += 6;        else if(ty1<tz0) index += 3;    }    return index;}/******************************************/int Menger::find_nextNode(real tx3, real ty3, real tz3,                          int x, int y, int z){    real min = Menger::Min(tx3,ty3,tz3);    if(tx3==min)      return x; //exits through yz plane    else if(ty3==min) return y; //xz    else              return z; //xy}/******************************************///used to determine time values for the nodes in the sub-treestatic const real f13 = 1.0/3.0;static const real f23 = 2.0/3.0;bool Menger::proc_subtree(real tx0, real ty0, real tz0,                          real tx3, real ty3, real tz3,                          int cite){    //does not intersect    if(tx3<=0||ty3<=0||tz3<=0)        return false;    //at lowest subdivision, traversal can end    if(cite==0)    {        //tvalue        tval = Menger::Max(tx0,ty0,tz0);        //normal - since ray traverses positively along all axi, normal component is always -1, it is negated to 1        //         the ray was reflected along that direction.        if(tx0==tval)      norm = vec3(flipx?1:-1,0,0);        else if(ty0==tval) norm = vec3(0,flipy?1:-1,0);        else               norm = vec3(0,0,flipz?1:-1);        return true;    }    //t values for mid-planes    real tx1 = f23*tx0+f13*tx3;    real tx2 = f13*tx0+f23*tx3;    real ty1 = f23*ty0+f13*ty3;    real ty2 = f13*ty0+f23*ty3;    real tz1 = f23*tz0+f13*tz3;    real tz2 = f13*tz0+f23*tz3;    //determine the first node in the sub-tree that is intersected    int currNode = find_firstNode(tx0,ty0,tz0,tx1,tx2,ty1,ty2,tz1,tz2);    do    {        //switch through the current node, if it is a cube that is not one of the wholes, process sub-tree        // proc_subtree being called only for certain cases is hard-coding the subdivion pattern of the menger cube,        // if you were to swap it round and have proc_subtree called elsewhere, you'd end up with a different pattern.        //advance to next node in the sub-tree, the index to step into next is hardcoded in the parameter for each exit plane        switch(currNode)        {            case 0: if(proc_subtree(tx0,ty0,tz0, tx1,ty1,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty1,tz1,1,3,9);                    break;            case 1: if(proc_subtree(tx1,ty0,tz0, tx2,ty1,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx2,ty1,tz1,2,4,10);                    break;            case 2: if(proc_subtree(tx2,ty0,tz0, tx3,ty1,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty1,tz1,27,5,11);                    break;            case 3: if(proc_subtree(tx0,ty1,tz0, tx1,ty2,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty2,tz1,4,6,12);                    break;            case 4:                    currNode = find_nextNode(    tx2,ty2,tz1,5,7,13);                    break;            case 5: if(proc_subtree(tx2,ty1,tz0, tx3,ty2,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty2,tz1,27,8,14);                    break;            case 6: if(proc_subtree(tx0,ty2,tz0, tx1,ty3,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty3,tz1,7,27,15);                    break;            case 7: if(proc_subtree(tx1,ty2,tz0, tx2,ty3,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx2,ty3,tz1,8,27,16);                    break;            case 8: if(proc_subtree(tx2,ty2,tz0, tx3,ty3,tz1, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty3,tz1,27,27,17);                    break;            case 9: if(proc_subtree(tx0,ty0,tz1, tx1,ty1,tz2, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty1,tz2,10,12,18);                    break;            case 10:                    currNode = find_nextNode(    tx2,ty1,tz2,11,13,19);                    break;            case 11:if(proc_subtree(tx2,ty0,tz1, tx3,ty1,tz2, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty1,tz2,27,14,20);                    break;            case 12:                    currNode = find_nextNode(    tx1,ty2,tz2,13,15,21);                    break;            case 13:                    currNode = find_nextNode(    tx2,ty2,tz2,14,16,22);                    break;            case 14:                    currNode = find_nextNode(    tx3,ty2,tz2,27,17,23);                    break;            case 15:if(proc_subtree(tx0,ty2,tz1, tx1,ty3,tz2, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty3,tz2,16,27,24);                    break;            case 16:                    currNode = find_nextNode(    tx2,ty3,tz2,17,27,25);                    break;            case 17:if(proc_subtree(tx2,ty2,tz1, tx3,ty3,tz2, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty3,tz2,27,27,26);                    break;            case 18:if(proc_subtree(tx0,ty0,tz2, tx1,ty1,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty1,tz3,19,21,27);                    break;            case 19:if(proc_subtree(tx1,ty0,tz2, tx2,ty1,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx2,ty1,tz3,20,22,27);                    break;            case 20:if(proc_subtree(tx2,ty0,tz2, tx3,ty1,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty1,tz3,27,23,27);                    break;            case 21:if(proc_subtree(tx0,ty1,tz2, tx1,ty2,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty2,tz3,22,24,27);                    break;            case 22:                    currNode = find_nextNode(    tx2,ty2,tz3,23,25,27);                    break;            case 23:if(proc_subtree(tx2,ty1,tz2, tx3,ty2,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx3,ty2,tz3,27,26,27);                    break;            case 24:if(proc_subtree(tx0,ty2,tz2, tx1,ty3,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx1,ty3,tz3,25,27,27);                    break;            case 25:if(proc_subtree(tx1,ty2,tz2, tx2,ty3,tz3, cite-1)) return true;                    currNode = find_nextNode(    tx2,ty3,tz3,26,27,27);                    break;            case 26:if(proc_subtree(tx2,ty2,tz2, tx3,ty3,tz3, cite-1)) return true;                    currNode = 27;                    break;        }    }while(currNode < 27); //currNode is set to 27 when leaving the sub-tree cell altogether,                           //aka doesn't intersect anything in this sub-tree    return false;}/*************************************************************************************/bool Menger::intersectSponge(const Menger::Sponge & o, Menger::Ray r, real & t, vec3 & n){    vec3 & p = r.p;    vec3 & v = r.v;    //the algorithm is much simpler if the ray is only ever traversed positively throug the axi    //exploiting the symmetry of the sponge, we can simply reflect in each axial plane with the same end-result    if(flipx = v.x<0)    {        p.x = -p.x;        v.x = -v.x;    }    if(flipy = v.y<0)    {        p.y = -p.y;        v.y = -v.y;    }    if(flipz = v.z<0)    {        p.z = -p.z;        v.z = -v.z;    }    //determine time values for intersection of each outer-cube plane    //to deal with axis-aligned rays easily, if the value is 0, it is very slightly increment    //so as not to effect outcome of rendering, but allow the algorithm to work properly    if(v.x==0) v.x=1e-10;    if(v.y==0) v.y=1e-10;    if(v.z==0) v.z=1e-10;    vec3 iv = v.inv(); //float3::inv() -> component wise multiplicative inverse aka (1/x,1/y,1/z)    real tx0 = (-0.5-p.x)*iv.x; real tx3 = (0.5-p.x)*iv.x;    real ty0 = (-0.5-p.y)*iv.y; real ty3 = (0.5-p.y)*iv.y;    real tz0 = (-0.5-p.z)*iv.z; real tz3 = (0.5-p.z)*iv.z;    //determine if ray intersects lowest cube    real t0 = Max(tx0,ty0,tz0);    real t1 = Min(tx3,ty3,tz3);    if(t0<t1 && t1>0)    {        //returns true if the ray intersects the sponge        if(proc_subtree(tx0,ty0,tz0,tx3,ty3,tz3, o.iterations))        {            //set n,t parameters to determined time value and normal            t = tval;            n = norm;            return true;        }    }    //ray has not intersected the sponge    return false;}