Menger Sponge in C++ using OpenGL

This past summer I was going through some old projects and came across a Menger Sponge visualizer that I wrote back in college. A Menger Sponge is simple fractal that has infinite surface area and encloses zero volume. The sponge is constructed in successive iterations and the first four iterations are rendered in the video below.

The sponge starts as a single cube that is segmented into twenty-seven equally sized cubes. The center cubes of each face and that of the parent cube are then discarded and the process is applied again to each of the remaining cubes. Visually, the process looks like the following:

The geometry of the process is straight forward. Starting with a cube’s origin, , and edge length, , each of the children’s attributes can be calculated. Each child’s edge length is given by . Each child’s origin given by . The constant represents a child’s relative position (e.g., ) to its parent.

The following implementation isn’t particularly well written, but it accomplishes the desired end result. The point and Cube classes achieve the logic that I’ve outlined above. Cube can be thought of as a tree structure that is generated upon instantiation. The visualize() method pumps out the desired OpenGL commands to produce the faces of the cubes.

#include <GL\glut.h>

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

//=================================================================================
//=================================================================================
class point
{
public:
	point(GLfloat x, GLfloat y, GLfloat z, point* ref = NULL);
	void visualize();

	GLfloat x,y,z;
};

point::point(GLfloat x, GLfloat y, GLfloat z, point* ref)
{
	this->x = x;
	this->y = y;
	this->z = z;

	if(ref != NULL)
	{
		this->x += ref->x;
		this->y += ref->y;
		this->z += ref->z;
	}
}

//=================================================================================
//=================================================================================

class Cube
{
public:
	Cube(point* origin, GLfloat edgelength, GLfloat depth);
	~Cube();

	void visualize();

private:
	void MakeFace(int i, int j, int k, int l);
	void ActAsContainer(point* o, GLfloat e, GLfloat d);
	void ActAsCube(point* o, GLfloat e);

	point** PointCloud;
	Cube** SubCubes;
};

Cube::Cube(point* origin, GLfloat edgelength, GLfloat depth)
{
	if(depth <= 1.0)
	{
		ActAsCube(origin, edgelength);
	} else {
		ActAsContainer(origin, edgelength, depth);
	}
}

Cube::~Cube()
{
	int i;

	if(PointCloud != NULL)
	{
		for(i = 0; i < 8; i++)
			delete PointCloud[i];
		delete[] PointCloud;
	}

	if(SubCubes != NULL)
	{
		for(i = 0; i < 20; i++)
			delete SubCubes[i];
		delete[] SubCubes;
	}
}

void Cube::ActAsCube(point* o, GLfloat e)
{
	GLfloat ne = e / 2.0;

	PointCloud = new point*[8];		// This is the actual physical cube coordinates;
	SubCubes = NULL;

	PointCloud[0] = new point( ne,  ne,  ne, o);	// net
	PointCloud[1] = new point( ne, -ne,  ne, o);	// set
	PointCloud[2] = new point(-ne,  ne,  ne, o);	// nwt
	PointCloud[3] = new point(-ne, -ne,  ne, o);	// swt
	PointCloud[4] = new point( ne,  ne, -ne, o);	// neb
	PointCloud[5] = new point( ne, -ne, -ne, o);	// seb
	PointCloud[6] = new point(-ne,  ne, -ne, o);	// nwb
	PointCloud[7] = new point(-ne, -ne, -ne, o);	// swb
}

void Cube::ActAsContainer(point* o, GLfloat e, GLfloat d)
{
	GLfloat ne = e / 3.0;

	SubCubes = new Cube*[20];	// These are the centers of each sub cube structure
	PointCloud = NULL;

	SubCubes[0] = new Cube(new point(-ne,  ne,  ne, o), ne, d-1.0);
	SubCubes[1] = new Cube(new point(0.0,  ne,  ne, o), ne, d-1.0);
	SubCubes[2] = new Cube(new point( ne,  ne,  ne, o), ne, d-1.0);
	SubCubes[3] = new Cube(new point( ne, 0.0,  ne, o), ne, d-1.0);
	SubCubes[4] = new Cube(new point( ne, -ne,  ne, o), ne, d-1.0);
	SubCubes[5] = new Cube(new point(0.0, -ne,  ne, o), ne, d-1.0);
	SubCubes[6] = new Cube(new point(-ne, -ne,  ne, o), ne, d-1.0);
	SubCubes[7] = new Cube(new point(-ne, 0.0,  ne, o), ne, d-1.0);
	
	SubCubes[8] = new Cube(new point( ne,  ne,  0.0, o), ne, d-1.0);
	SubCubes[9] = new Cube(new point( ne, -ne,  0.0, o), ne, d-1.0);
	SubCubes[10] = new Cube(new point(-ne, ne,  0.0, o), ne, d-1.0);
	SubCubes[11] = new Cube(new point(-ne, -ne,  0.0, o), ne, d-1.0);
	
	SubCubes[12] = new Cube(new point(-ne,  ne, -ne, o), ne, d-1.0);
	SubCubes[13] = new Cube(new point(0.0,  ne, -ne, o), ne, d-1.0);
	SubCubes[14] = new Cube(new point( ne,  ne, -ne, o), ne, d-1.0);
	SubCubes[15] = new Cube(new point( ne, 0.0, -ne, o), ne, d-1.0);
	SubCubes[16] = new Cube(new point( ne, -ne, -ne, o), ne, d-1.0);
	SubCubes[17] = new Cube(new point(0.0, -ne, -ne, o), ne, d-1.0);
	SubCubes[18] = new Cube(new point(-ne, -ne, -ne, o), ne, d-1.0);
	SubCubes[19] = new Cube(new point(-ne, 0.0, -ne, o), ne, d-1.0);
}

void Cube::MakeFace(int i, int j, int k, int l)
{
		glVertex3f(PointCloud[i]->x, PointCloud[i]->y, PointCloud[i]->z);
		glVertex3f(PointCloud[j]->x, PointCloud[j]->y, PointCloud[j]->z);
		glVertex3f(PointCloud[k]->x, PointCloud[k]->y, PointCloud[k]->z);
		glVertex3f(PointCloud[l]->x, PointCloud[l]->y, PointCloud[l]->z);

}

void Cube::visualize()
{
	int i;

	if(PointCloud != NULL)
	{
		glBegin(GL_QUADS);
			glColor3f(1.0,0.0,0.0);// top
			MakeFace(0,2,3,1);
			glColor3f(0.0,1.0,1.0);//bottom
			MakeFace(4,6,7,5);
			
			glColor3f(0.0,1.0,0.0);// north
			MakeFace(0,2,6,4);
			glColor3f(1.0,0.0,1.0);// south
			MakeFace(1,3,7,5);
			
			glColor3f(0.0,0.0,1.0);//east
			MakeFace(0,4,5,1);
			glColor3f(1.0,1.0,0.0);// west
			MakeFace(2,6,7,3);
		glEnd();
	}

	if(SubCubes != NULL)
	{
		for(i = 0; i < 20; i++)
		{
			SubCubes[i]->visualize();
		}
	}
}

The implementation of the program is your run-of-the-mill OpenGL boilerplate. The application takes in an argument dictating what order of sponge it should produce. It sets up the camera and positions the sponge at the origin. The sponge is left stationary, while the camera is made to orbit upon each display(). On idle(), a redisplay message is sent back to the OpenGL system in order to achieve the effect that the sponge is spinning.

//=================================================================================
//=================================================================================
Cube* MengerCube;

void idle()
{
	glutPostRedisplay();
}

void display()
{
	static GLfloat rtri = 0.0;
	
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	gluLookAt(1.0,1.0,1.0, 0.0,0.0,0.0,0.0,1.0,0.0);
	glRotatef((rtri+=0.932), 1.0, 0.5, -1.0);

	MengerCube->visualize();

	glutSwapBuffers();
}

void reshape(int w, int h)
{
	glViewport(0,0,w,h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-8.0, 8.0,-8.0, 8.0,-8.0, 8.0);
}

void init()
{
	glShadeModel(GL_SMOOTH);
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glClearDepth(1.0f);
	glEnable(GL_DEPTH_TEST);
	glColor3f(1.0, 1.0, 1.0);
}

GLfloat getDepth(char* depth)
{
	int k = atoi(depth);
	if(k <= 1) return 1.0;
	else if (k>= 5) return 5.0;
	else return (GLfloat) k;
}

int main(int argc, char* argv[])
{
	GLfloat depth;
	bool viewAsPointCloud = false;
	point origin(0.0, 0.0, 0.0);

	printf("%d\n",argc);

	switch(argc)
	{
	case 2:
		depth = getDepth(argv[1]);
		break;
	default:
		depth = 2.0;
		break;
	}

	MengerCube = new Cube(&origin, 8.0, depth);

	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glEnable(GL_DEPTH_TEST);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(0,0);
	glutCreateWindow(*argv);
	glutReshapeFunc(reshape);
	glutDisplayFunc(display);
	glutIdleFunc(idle);
	init();
	glutMainLoop();

	delete MengerCube;
}