VoxelCube.cpp

Go to the documentation of this file.

#include "../common/standard.h"
#include "VoxelCube.h"
#include "GLBitmapBuffer.h"
#include "MarchingCubesTables.h"
#include "GLutils.h"

#include <GL/glu.h>

#include <map>


std::vector< Vertex > VoxelCube::calculateVoxelCorners (int xpos, int ypos, int zpos) const
{
    double halfVoxelSize = voxelSize / 2.;
    double halfCubeSize = voxelsPerSide * halfVoxelSize;

    Vertex origin = Vertex (-halfCubeSize + center[0], -halfCubeSize + center[1], -halfCubeSize + center[2]);

    origin[0] += voxelSize * xpos;
    origin[1] += voxelSize * ypos;
    origin[2] += voxelSize * zpos;

    std::vector< Vertex > rtn;
    for (int z = 0; z < 2; z++)
    {
        for (int y = 0; y < 2; y++)
        {
            for (int x = 0; x < 2; x++)
            {
                Vertex pt = origin;
                pt[0] += voxelSize*x;
                pt[1] += voxelSize*y;
                pt[2] += voxelSize*z;

                rtn.push_back (pt);
            }
        }
    }

    return rtn;
}


void VoxelCube::setTransparent(SurfaceVoxelList::iterator i)
{
    VoxelPos &p = *i;
    int index = PosToIndex (p[0], p[1], p[2]);
    solids[index] = false;
}


void VoxelCube::exposeSurroundingVoxels (SurfaceVoxelList::iterator i)
{
    VoxelPos &p = *i;

    unsigned int x = p[0];
    unsigned int y = p[1];
    unsigned int z = p[2];

/*
 * Enumerate the possible surrounding voxels
 */
    VoxelPos x1(x-1, y,   z);
    VoxelPos y1(x,   y-1, z);
    VoxelPos z1(x,   y,   z-1);
    VoxelPos x2(x+1, y,   z);
    VoxelPos y2(x,   y+1, z);
    VoxelPos z2(x,   y,   z+1);

/* 
 * Add surrounding voxels to the SVL. 
 */
    unsigned int maxVoxel = voxelsPerSide - 1;

    if (x > 0 && !isTransparent(x1)) svl.push_back (x1);
    if (y > 0 && !isTransparent(y1)) svl.push_back (y1);
    if (z > 0 && !isTransparent(z1)) svl.push_back (z1);

    if (x < maxVoxel && !isTransparent(x2)) svl.push_back (x2);
    if (y < maxVoxel && !isTransparent(y2)) svl.push_back (y2);
    if (z < maxVoxel && !isTransparent(z2)) svl.push_back (z2);
}



void VoxelCube::setFinalColor (SurfaceVoxelList::iterator i, RgbByte &color)
{
    VoxelPos &pos = *i;

    ColoredVoxelPos cvp(pos, color);

    svlColored.push_back (cvp);
}



void VoxelCube::fillItemBufferArray (Viewpoint &image)
{
    GLBitmapBuffer bitbuff;
    int width = image.getWidth();
    int height = image.getHeight();

/*
 * Make our bitmap buffer, and set up the image's perspective
 */
    bitbuff.create(width, height);
    bitbuff.select();

    // set the projection matrix
    glMatrixMode(GL_PROJECTION);
    glLoadMatrixd(image.getProjectionMatrix());

    // set the modelview matrix
    glMatrixMode(GL_MODELVIEW);
    glLoadMatrixd(image.getModelviewMatrix());

    // set the viewport
    glViewport (0, 0, width, height);
    glEnable (GL_DEPTH_TEST);        
    glEnable (GL_CULL_FACE);
/*
 * Draw to the bitmap buffer
 */
    glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    double halfVoxelSize = voxelSize / 2.;
    double halfCubeSize = voxelsPerSide * halfVoxelSize;

    glMatrixMode (GL_MODELVIEW);
    glPushMatrix();
        glTranslated(center[0], center[1], center[2]);
        glTranslated (-halfCubeSize, -halfCubeSize, -halfCubeSize);
        glTranslated (halfVoxelSize, halfVoxelSize, halfVoxelSize);
        
            SurfaceVoxelList::iterator i;
            unsigned int id;
            for (i = svl.begin(); i != svl.end(); i++)
            {
                VoxelPos &pos = *i;
                glPushMatrix();
                    glTranslated (pos[0]*voxelSize, pos[1]*voxelSize, pos[2]*voxelSize);

                    id = PosToIdentifier(pos[0], pos[1], pos[2]);
                    glColor3ubv ((GLubyte*)&id);
                    
                    drawCube (voxelSize, GL_QUADS);
                glPopMatrix();
            }

    glPopMatrix();

    glFinish();

    // copy the bitmap buffer to the image's item buffer
    glReadPixels (
        0, 0, 
        width, height,
        GL_RGBA, GL_UNSIGNED_BYTE,
        &image.imageBuffer[0]);

    bitbuff.unselect();
}



ColorVertex Interpolate(Vertex &v1, Vertex &v2, double tVal1, double tVal2, float isoLevel)
{
        ColorVertex interpolation;
        float mu;

        mu = float((isoLevel - tVal1))/(tVal2 - tVal1);
        interpolation.pos[0] = v1[0] + mu*(v2[0] - v1[0]);
        interpolation.pos[1] = v1[1] + mu*(v2[1] - v1[1]);
        interpolation.pos[2] = v1[2] + mu*(v2[2] - v1[2]);

        return interpolation;
}



ColorVertex VoxelCube::calculateIntersection(ColVector< 3, unsigned int> &curOrigin, unsigned int nEdgeNo, float isoLevel, ColorMap &colorMap)
{
        float x1, y1, z1, x2, y2, z2;
        unsigned int v1x = curOrigin[0], v1y = curOrigin[1], v1z = curOrigin[2];
        unsigned int v2x = curOrigin[0], v2y = curOrigin[1], v2z = curOrigin[2];
        
        switch (nEdgeNo)
        {
        case 0:
                v2y += 1;
                break;
        case 1:
                v1y += 1;
                v2x += 1;
                v2y += 1;
                break;
        case 2:
                v1x += 1;
                v1y += 1;
                v2x += 1;
                break;
        case 3:
                v1x += 1;
                break;
        case 4:
                v1z += 1;
                v2y += 1;
                v2z += 1;
                break;
        case 5:
                v1y += 1;
                v1z += 1;
                v2x += 1;
                v2y += 1;
                v2z += 1;
                break;
        case 6:
                v1x += 1;
                v1y += 1;
                v1z += 1;
                v2x += 1;
                v2z += 1;
                break;
        case 7:
                v1x += 1;
                v1z += 1;
                v2z += 1;
                break;
        case 8:
                v2z += 1;
                break;
        case 9:
                v1y += 1;
                v2y += 1;
                v2z += 1;
                break;
        case 10:
                v1x += 1;
                v1y += 1;
                v2x += 1;
                v2y += 1;
                v2z += 1;
                break;
        case 11:
                v1x += 1;
                v2x += 1;
                v2z += 1;
                break;
        }

        x1 = v1x*voxelSize;
        y1 = v1y*voxelSize;
        z1 = v1z*voxelSize;
        x2 = v2x*voxelSize;
        y2 = v2y*voxelSize;
        z2 = v2z*voxelSize;

        unsigned int nPointsInXDirection = (voxelsPerSide);
        unsigned int nPointsInSlice = nPointsInXDirection*(voxelsPerSide);
        int val1 = !isTransparent(v1x, v1y, v1z);
    int val2 = !isTransparent(v2x, v2y, v2z);

    ColorVertex intersection = Interpolate(
        Vertex(x1, y1, z1),
        Vertex(x2, y2, z2),
        val1, val2, isoLevel);

    VoxelPos va(v1x, v1y, v1z);
    VoxelPos vb(v2x, v2y, v2z);

    RgbByte color1 (0,0,0);
    RgbByte color2 (0,0,0);

    ColorMap::iterator i;
#if 0
    i = colorMap.find(va);
    if (i != colorMap.end()) color1 = (*i).second;

    i = colorMap.find(vb);
    if (i != colorMap.end()) color2 = (*i).second;

    intersection.color = LERP(0.5,
        color1,
        color2);
#else
    ColorMap::iterator i_first = colorMap.find(va);
    ColorMap::iterator i_second = colorMap.find(vb);

    if (i_first != colorMap.end())
    {
        if (i_second != colorMap.end())
            intersection.color = LERP(0.5, (*i_first).second, (*i_second).second);
        else
            intersection.color = (*i_first).second;
    }
    else
    {
        if (i_second != colorMap.end())
            intersection.color = (*i_second).second;
        else
            intersection.color = RgbByte(1,0,0);
    }
#endif


    
        return intersection;
}





unsigned int VoxelCube::GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
{
        switch (nEdgeNo) {
        case 0:
                return PosToIndex(nX, nY, nZ) + 1;
        case 1:
                return PosToIndex(nX, nY + 1, nZ);
        case 2:
                return PosToIndex(nX + 1, nY, nZ) + 1;
        case 3:
                return PosToIndex(nX, nY, nZ);
        case 4:
                return PosToIndex(nX, nY, nZ + 1) + 1;
        case 5:
                return PosToIndex(nX, nY + 1, nZ + 1);
        case 6:
                return PosToIndex(nX + 1, nY, nZ + 1) + 1;
        case 7:
                return PosToIndex(nX, nY, nZ + 1);
        case 8:
                return PosToIndex(nX, nY, nZ) + 2;
        case 9:
                return PosToIndex(nX, nY + 1, nZ) + 2;
        case 10:
                return PosToIndex(nX + 1, nY + 1, nZ) + 2;
        case 11:
                return PosToIndex(nX + 1, nY, nZ) + 2;
        default:
                // Invalid edge no.
                throw ArgError("GetEdgeID(): invalid edge number");
        }
}

struct TRIANGLE {
        unsigned int pointID[3];
};

typedef std::map< unsigned int, ColorVertex > IdToVertexMap;
typedef std::vector< TRIANGLE > TriangleArray;



void VoxelCube::constructColorMap (ColorMap &colorMap)
{
    if (!svlColored.size())
        return;

    ColoredSurfaceVoxelList::iterator i;
    for (i = svlColored.begin(); i != svlColored.end(); i++)
    {
        ColoredVoxelPos pos = (*i);
        
        VoxelPos p;
        p[0] = pos[0]; p[1] = pos[1]; p[2] = pos[2];
        colorMap[p] = pos.color;
    }
}



void VoxelCube::makeMesh(float isoLevel, QEM::ProgMesh &mesh)
{
    ColorMap colorMap;
    constructColorMap (colorMap);
    
    unsigned int nCells = voxelsPerSide - 1;

    IdToVertexMap idVertexMap;
    TriangleArray tris;

    for (int z = 0; z < nCells; z++)
    {
        for (int y = 0; y < nCells; y++)
        {
            for (int x = 0; x < nCells; x++)
            {
                unsigned int tableIndex = 0;
                ColVector<3, unsigned int> curOrigin (x,y,z);
                
                if (isTransparent (x,   y ,  z))   tableIndex |= 1;   // 1 << 0
                if (isTransparent (x,   y+1, z))   tableIndex |= 2;   // 1 << 1
                if (isTransparent (x+1, y+1, z))   tableIndex |= 4;   // 1 << 2
                if (isTransparent (x+1, y,   z))   tableIndex |= 8;   // 1 << 3
                if (isTransparent (x,   y,   z+1)) tableIndex |= 16;  // 1 << 4
                if (isTransparent (x,   y+1, z+1)) tableIndex |= 32;  // 1 << 5
                if (isTransparent (x+1, y+1, z+1)) tableIndex |= 64;  // 1 << 6
                if (isTransparent (x+1, y,   z+1)) tableIndex |= 128; // 1 << 7
                
                if (edgeTable[tableIndex] == 0)
                    continue;
                
#define EDGETABLE_HAS(n) edgeTable[tableIndex] & (1 << n)
                
                // Now create a triangulation of the isosurface in this
                // cell.
                if (EDGETABLE_HAS (3)) {
                    ColorVertex pt = calculateIntersection (curOrigin, 3, isoLevel, colorMap);
                    unsigned int id = GetEdgeID(x, y, z, 3);
                    idVertexMap[id] = pt;
                }
                if (EDGETABLE_HAS(0)) {
                    ColorVertex pt = calculateIntersection (curOrigin, 0, isoLevel, colorMap);
                    unsigned int id = GetEdgeID(x, y, z, 0);
                    idVertexMap[id] = pt;
                }
                if (EDGETABLE_HAS(8)) {
                    ColorVertex pt = calculateIntersection (curOrigin, 8, isoLevel, colorMap);
                    unsigned int id = GetEdgeID(x, y, z, 8);
                    idVertexMap[id] = pt;
                }
                
                if (x == nCells - 1) {
                    if (EDGETABLE_HAS(2)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 2, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 2);
                        idVertexMap[id] = pt;
                    }
                    if (EDGETABLE_HAS(11)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 11, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 11);
                        idVertexMap[id] = pt;
                    }
                }
                if (y == nCells - 1) {
                    if (EDGETABLE_HAS(1)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 1, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 1);
                        idVertexMap[id] = pt;
                    }
                    if (EDGETABLE_HAS(9)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 9, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 9);
                        idVertexMap[id] = pt;
                    }
                }
                if (z == nCells - 1) {
                    if (EDGETABLE_HAS(4)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 4, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 4);
                        idVertexMap[id] = pt;
                    }
                    if (EDGETABLE_HAS(7)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 7, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 7);
                        idVertexMap[id] = pt;
                    }
                }
                if ((x==nCells - 1) && (y==nCells - 1))
                {
                    if (EDGETABLE_HAS(10)) 
                    {
                        ColorVertex pt = calculateIntersection (curOrigin, 10, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 10);
                        idVertexMap[id] = pt;
                    }
                }
                if ((x==nCells - 1) && (z==nCells - 1))
                {
                    if (EDGETABLE_HAS(6)) 
                    {
                        ColorVertex pt = calculateIntersection (curOrigin, 6, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 6);
                        idVertexMap[id] = pt;
                    }
                }
                if ((y==nCells - 1) && (z==nCells - 1))
                {
                    if (EDGETABLE_HAS(5)) {
                        ColorVertex pt = calculateIntersection (curOrigin, 5, isoLevel, colorMap);
                        unsigned int id = GetEdgeID(x, y, z, 5);
                        idVertexMap[id] = pt;
                    }
                }
                
                for (unsigned int i = 0; triTable[tableIndex][i] != -1; i += 3) {

                    TRIANGLE triangle;
                    unsigned int pointID0, pointID1, pointID2;
                    pointID0 = GetEdgeID(x, y, z, triTable[tableIndex][i]);
                    pointID1 = GetEdgeID(x, y, z, triTable[tableIndex][i+1]);
                    pointID2 = GetEdgeID(x, y, z, triTable[tableIndex][i+2]);
                    
                    triangle.pointID[0] = pointID0;
                    triangle.pointID[1] = pointID1;
                    triangle.pointID[2] = pointID2;
                                        
                    tris.push_back(triangle);
                }                

            }
            
        }
        
    }

    IdToVertexMap::iterator vertex_iter;

    std::map< unsigned int, unsigned int> oursToTheirs;

    float halfWidth = (float)(voxelsPerSide * voxelSize) / 2;

    for (vertex_iter = idVertexMap.begin(); vertex_iter != idVertexMap.end(); vertex_iter++)
    {
        ColorVertex &v = (*vertex_iter).second;
        
        unsigned int id = mesh.addVertexWithColor (v.pos[0], v.pos[1], v.pos[2], v.color);
        oursToTheirs[(*vertex_iter).first] = id;
        MTRACE ("added vertex %.3d: %.2lf %.2lf %.2lf\n",
            id, v.pos[0], v.pos[1], v.pos[2]);
    }

    for (int i = 0; i < tris.size(); i++)
    {
        unsigned int a = tris[i].pointID[0];
        unsigned int b = tris[i].pointID[1];
        unsigned int c = tris[i].pointID[2];

        MTRACE("Trying to add face ids (%.3d -> %3d) (%.3d -> %3d) (%.3d -> %3d)\n",
            a, oursToTheirs[a],
            b, oursToTheirs[b],
            c, oursToTheirs[c]);

        if (a != b && b != c && a != c)
        {
            try {
                mesh.addFace(
                    oursToTheirs[a],
                    oursToTheirs[b],
                    oursToTheirs[c]);
            }
            catch (...)
            {
                MTRACE("caught\n");
            }
        }
    }

}


void VoxelCube::draw(bool bTranslateToCenter)
{
    double halfVoxelSize = voxelSize / 2.;
    double halfCubeSize = voxelsPerSide * halfVoxelSize;
    
    glPushAttrib (GL_ENABLE_BIT);

        glDisable (GL_TEXTURE_2D);

        glMatrixMode (GL_MODELVIEW);
        glPushMatrix();
            if (bTranslateToCenter)
                glTranslated(center[0], center[1], center[2]);

            glTranslated(-halfCubeSize, -halfCubeSize, -halfCubeSize);
            glTranslated(halfVoxelSize, halfVoxelSize, halfVoxelSize);

            // if there are colored voxel in our colored surface voxel list,
            // then let's draw that
            if (svlColored.size())
            {
                ColoredSurfaceVoxelList::iterator i;
                for (i = svlColored.begin(); i != svlColored.end(); i++)
                {
                    ColoredVoxelPos pos = (*i);
                    
                    glPushMatrix();
                    
                        glTranslated (
                            pos[0] * voxelSize,
                            pos[1] * voxelSize,
                            pos[2] * voxelSize);
                    
                        glColor3ub (pos.color[0], pos.color[1], pos.color[2]);
                        drawCube (voxelSize, GL_QUADS);

                    glPopMatrix();
                }
            }
            else // no colors are assigned yet, so we'll draw the svl with random colors
            {
                SurfaceVoxelList::iterator i;
                for (i = svl.begin(); i != svl.end(); i++)
                {
                    VoxelPos pos = (*i);
                    
                    glPushMatrix();
                    
                        glTranslated (
                            pos[0] * voxelSize,
                            pos[1] * voxelSize,
                            pos[2] * voxelSize);
                    
                        srand (pos[0]+pos[1]*voxelsPerSide+pos[2]*voxelsPerSide*voxelsPerSide);

                        glColor3f (
                            (rand()%100)/100.,
                            (rand()%100)/100.,
                            (rand()%100)/100.);
                        drawCube (voxelSize, GL_QUADS);

                    glPopMatrix();
                }
            }
        glPopMatrix();

    glPopAttrib (); // GL_ENABLE_BIT
}




double VoxelCube::getMinimumFarPlane (Vertex COP, Vertex lookatPoint) const
{
    ColVector3 viewDirection = lookatPoint - COP;
    viewDirection.normalize();

    double rtn = 0.0;

    for (int x = -1; x <= 1; x+= 2)
    {
        for (int y = -1; y <= 1; y+= 2)
        {
            for (int z = -1; z <= 1; z += 2)
            {
                Vertex corner = center;
                corner[0] += (x * voxelSize * voxelsPerSide) / 2.;
                corner[1] += (y * voxelSize * voxelsPerSide) / 2.;
                corner[2] += (z * voxelSize * voxelsPerSide) / 2.;

                ColVector3 diff = corner - COP;

                ColVector3 diffDirection = diff;
                diffDirection.normalize();

                double cos = diffDirection.dotProduct(viewDirection);

                if (cos > 0)
                {
                    double dist = cos * diff.getLength();
                    rtn = MAX (rtn, dist);
                }
            }
        }
    }

    return rtn;
}


bool VoxelCube::isTransparent (const VoxelPos &p) const
{
    return isTransparent (p[0], p[1], p[2]);
}

bool VoxelCube::isTransparent (int x, int y, int z) const
{
    return !solids[PosToIndex (x,y,z)];
}



void VoxelCube::setVoxelsPerSide (int vps)
{
    if (vps < 1 || vps > 255) throw ArgError("VoxelCube::setVoxelsPerSide(): voxels per side must be between 1 and 255");

    voxelsPerSide = vps;

    svl.clear();
    svlColored.clear();

   
    // add the 8 faces of the voxel cube to the SVL
    int max = voxelsPerSide - 1;
    for (int row = 0; row < voxelsPerSide; row++)
        for (int col = 0; col < voxelsPerSide; col++)
        {
            svl.push_back (VoxelPos(0,row,col));
            svl.push_back (VoxelPos(row,0,col));
            svl.push_back (VoxelPos(row,col,0));

            svl.push_back (VoxelPos(max,row,col));
            svl.push_back (VoxelPos(row,max,col));
            svl.push_back (VoxelPos(row,col,max));
        }


    // initialize our list of solid (vs. transparent) bools
    solids.resize(vps*vps*vps);

    // all of our voxels are initially solid;
    // initialize our "solids" array to all "true"
    std::fill(solids.begin(), solids.end(), true);
}

unsigned int VoxelCube::PosToIdentifier (int x, int y, int z) const
{
    unsigned int rtn = PosToIndex(x,y,z) + 1;
    return rtn;
}


unsigned int VoxelCube::PosToIndex (int x, int y, int z) const
{ 
    if (x >= voxelsPerSide ||
        y >= voxelsPerSide ||
        z >= voxelsPerSide)
    {
        throw ArgError ("VoxelCube::getIndex(): you've asked for voxel (%d,%d,%d), but there are only %d voxels per side!",
            x,y,z,voxelsPerSide);
    }

    return x + y*voxelsPerSide + z*voxelsPerSide*voxelsPerSide; 
}




Rect< double > VoxelCube::getBoundingRect (std::vector< Vertex > &vertices, const GLdouble modelMatrix[16], const GLdouble projMatrix[16], const GLint viewport[4]) const
{    
    Rect< double > bounds;
    
    double garbage;
    GLenum rtn = gluProject (
        vertices[0][0], vertices[0][1], vertices[0][2],
        modelMatrix,
        projMatrix,
        viewport,
        &bounds.left, &bounds.bottom,
        &garbage);

    if (rtn != GL_TRUE)
        throw Error("VoxelCube::getBoundingRect(): gluProject failure.");
    
    bounds.right = bounds.left;
    bounds.top = bounds.bottom;
    
   /*
    * Let's make a bounding rectangle for the projected voxel corners
    */

    for (int corner = 1; corner < 8; corner++)
    {
        double curx, cury;

        rtn = gluProject (
            vertices[corner][0], vertices[corner][1], vertices[corner][2],
            modelMatrix,
            projMatrix,
            viewport,
            &curx, &cury,
            &garbage);
        if (rtn != GL_TRUE)
            throw Error("VoxelCube::getBoundingRect(): gluProject failure.");
        
        bounds.expandToEnclosePoint (curx, cury);
    }

    // let's be conservative, and enlarge the bounding box by a pixel
    bounds.bottom -= 1;
    bounds.left   -= 1;
    bounds.right  +=1;
    bounds.top    +=1;
    
    Rect< double > screen(
        viewport[0], 
        viewport[2]-1,
        viewport[3]-1,
        viewport[1]);
    
    bounds = bounds.intersectRect (screen);
    return bounds;

}

---