How the Marching Cubes Algorithm Works
Leslie Aine
Full-Stack Developer | JavaScript | React | Ruby on Rails | Redux | GitHub| Postgresql | HTML5 | Tailwind | CSS | Creative | Swimmer | Tennis
Remember your Maths Geometry class where you were asked to ‘shade the unwanted region?’ How about when you were taught the equation of a circle? And an ellipse? The Algebra class about Linear Interpolation? Lookup tables? Vectors?
If you took any basic Maths classes, these are familiar topics… or words, if the depth of the concepts is a thing of the past for you. And that is ok. Unless of course, you intend to develop an algorithm that can generate a smooth 3D surface from volumetric data, like the one CT and MRI scans use. In this case, you might want to track back and hammer these concepts in your mind.
Or… you are about 38 years late as someone has already developed this ground-breaking algorithm that proved very influential in coming up with the current 3D visualization techniques.
So, before I explain how all these classes can help you come down to this important algorithm, let us talk about life before 1987(when these geniuses, William Lorensen and Harvey Cline, came up with the Marching Cubes Algorithm). Were there no CT scans? MRIs? 3D models from data?
Actually, there were 3D models, but these were very raw, rigid, and unsophisticated as visualizing data was often in 2D. CT scans, MRIs, and other imaging equipment also existed but used only 2D visuals. The thing is, you had to be really, really good at visualizing, like an architect, or like in your Technical Drawing class because these scans came up with 2D visualizations. With these, the doctors, nurses, or whoever had to look at these pictures would have to piece the 2D visualizations together to come up with how the 3D image would look like… IN THEIR MINDS!??
Or, in older video games, terrains, and surfaces are so blocky and ugly! A classic example is Minecraft. Now, Minecraft has maintained this blocky image structure for aesthetic reasons among others, but its initial blocky structure was not deliberate. We are all familiar with pixels, but are we familiar with a voxel?? A voxel is the 3D equivalent of a pixel(easiest description I can think of) and there were some voxel-based video games, eg Outcast, or Delta Force. These games were blocky, had jagged terrain, the environment was not as pleasing or realistic as the smooth surfaces we see in today’s games. Voxel-based games like Astroneer rely on voxel terrains, but they use algorithms like Marching Cubes to create smooth, rounded surfaces for planets, caves, etc which result in visually appealing and seamless environments.
Ok, so we can already see how influential and important developing this algorithm has been to our world(why didn't someone just come up with this sooner????). Let us get into the ideas behind this very smart and cool piece of code.
In programming, you are taught that if you want to solve a big problem, you want to break it down to the simplest and basic form, solve that, and keep building on these simpler form solutions. So, for a 3D visualization problem, where each space can be treated as a cube, the simplest form of the cube is its 2D partner, the square.
Therefore before developing the Marching Cube Algorithm, the Marching Square Algorithm was created first. Building on its basics, everything else was a matter of working with a cube instead of a square.
And what about this square?
So, if I was drawing a circle, it would be easy as shown below. The equation of a circle, or ellipse is implicit. You can draw the circle and either shade within or outside it to make it stand out. Meaning, there is a boundary for this shading.
If I had a circle of radius 5cm at a coordinate position(0,0) on my grid, I would have the circle equation(or for our boundary in this case) as ;
From (x-h)^2 + (y-k)^2 = r^2
where (h,k) = (0,0) and r = 5
We now have our equation as ;
x^2 + y^2 = 25
Meaning any point where x2 + y2 > 25 is outside our circle and x2 + y2 < 25 is inside our circle. So, we can just shade this region within our circle since we have a boundary.
If you were coding, you could simply assign points within the circle a different color, say red, from the ones outside this boundary, say, blue, all this based on the boundary from the equation. So now you have a simple circle shaded red in a sea of blue.
Simple, right?! But if you were given a large dataset, would it come with an equation to define a boundary? For example, when handed a CT scan image, without anything shaded yet since we have no absolute boundary to make markings, with a grid of pixel values, with each pixel representing the density of tissue at that point, does it come with a function to help map out a bone, or a muscle, or an organ? Besides, do everybody’s organs conform to the same pattern, both in size and shape, to develop a function that can outline them from this CT scan?
The answer is no. And this is where our marching squares algorithm comes in. Even with a large abstract dataset, this algorithm is capable of making this distinction for you and creating these boundaries for you. The data does not directly give you a boundary but that boundary can be implied by your data and this algorithm will map it out for you.
You could be asking, how? Remember how we said our CT scan is made up of a grid of pixels? So, each pixel would represent tissue density at that particular point(density is the measure for a scan, if it were a weather map, we would be working with temperature). So since each pixel represents density at a point, we expect the pixel intensity value to keep changing at different points depending on the density of that particular point. If we knew a specific value that corresponds to the density of a bone, our problem is already half solved.
Lucky for us, these values are known and a bone density threshold value is between 700HU to 3000HU, depending on the bone tissue. This value is fed to our algorithm and is referred to as a threshold value.
For context, if we were working again with a weather map, it would be up to us to look for what to map out, maybe I want to map out regions with a temperature of 85?C, we would set 85 as our threshold value, and so on.
The algorithm then processes each square of the grid of pixel values in the scan. Pause. Yes, each square, hence the name, marching square.
What does it do to each square to determine a boundary? It checks if the pixel value(or in our case, density value) at each corner of a square is greater than or less than our threshold value.
This brings us to the next important part of our algorithm. I would even dare to say, the part that makes it special. Now, some of our squares will have all corners above or below our threshold value. These squares are obvious about where they lie even without a boundary. Those whose corners are all greater than our threshold value can be considered to be outside our not-yet-outlined ‘boundary’ and vice versa.
But what about those that have one, two, or three corners but not all four corners adhering to the same above or below the threshold value? Something like this below;
So, we would assume that if point A is above our threshold but points B, C and D are not, then our boundary must be FG as shown below;
But how can we determine where points F or G accurately lie on their respective lines? Simple, linear interpolation!
领英推荐
Linear interpolation combined with binary encoding and using Lookup Tables can help find the accurate locations of points F and G quite fast.
The need to use lookup tables and binary encoding before performing the actual linear interpolation is because unlike you and I, the computer does not automatically know to draw the line FG. So it has to first determine how the boundary will pass through a square. With us, we would make a fair assumption like we did assuming that line FG is how our boundary would pass through the square but the computer needs a lookup table to come up with this assumption.
Here’s how;
So, the binary code for this square is 1001 (from top-left to bottom-left, moving clockwise).
Using the lookup table
Each combination has a predefined way for the boundary to pass through the square. This information is stored in a lookup table.
This lookup table allows the algorithm to quickly decide how the boundary should pass through the square, without needing to recalculate the boundary for each square from scratch.
Once the algorithm determines how the boundary passes through a square using the lookup table, the next step is to find the exact points where the boundary crosses the edges of the square.
This is where your basic linear interpolation comes in. I know most of us will remember this interpolation, but I will still explain it simply below;
Let’s say the algorithm determines that the boundary crosses the top edge of a square between the top-left and top-right corners.
The boundary line (representing 1000HU) must cross the top edge somewhere between these two corners. Linear interpolation helps find the exact position.
Here’s how it works:
2. Calculate the ratio of how far the boundary is between the two points:
3. Now, apply this ratio to the horizontal distance between the corners (for example, if the corners are 1 unit apart):
By performing this interpolation for an edge where the boundary crosses, the algorithm finds the exact position of the boundary line along the edge of the square.
Like we suggested that our squares are marching, the algorithm repeats this process for each square until it finds the exact positions of our boundary line along the edges of the squares, allowing it to draw a smooth boundary line across the grid.
As you can already probably see, the more pixels, the more the squares, the more accurate the boundary lines can be and hence the smoother your edges will be(Probably why images with higher pixels are smooth and not so blocky and ugly like the ones with low pixels???).
Ok, so we have solved the problem in a 2D space, what about 3D? Well, in 3D visualization, instead of a square, we have a cube. So instead of a grid of squares, it is a grid of cubes, each cube with 8 corner points(instead of four for a square).
In this case, we are determining how a surface, not a boundary line, crosses a cube. Because in a 3D space, one line is not sufficient to create a boundary, but a plane, or surface could as illustrated below.
Each corner is binary encoded as the squares, and in this case, we use a 3D lookup table to determine how a surface(not a line) crosses a cube so instead of finding boundary lines on a 2D plane, Marching Cubes extracts triangles that form the surface in 3D space. These triangles are then pieced together to form a smooth, continuous 3D surface.
Like in Marching Squares, the algorithm uses linear interpolation along the edges of the cube to find the exact points where the surface crosses. This allows the extracted surface to be smooth and accurate, even in coarse datasets.
And voila! This is the logic behind the algorithm that was a significant milestone in the evolution of 3D visualization!
Before Marching Cubes, there were 3D models, but they were limited to wireframes, rough polygons, or blocky voxel-based models. Surface details were often missing or approximated, and creating smooth 3D surfaces from real-world data (e.g. medical scans) was extremely difficult. Marching Cubes revolutionized 3D modeling by enabling the automated extraction of smooth surfaces from volumetric data, ushering in a new era of medical imaging, scientific visualization, and 3D graphics.
If you made it this far, you are obviously very inquisitive, curious, and enjoy learning so please feel free to leave comments or suggestions on improvements or on what you would like me to write about next.