Login Form

OpenCL Marching Cubes

OpenCL Marching Cubes


(tradução breve) 


Get sample source code.


CMSoft brings a very versatile and useful tool, Marching Cubes adapted to GPU acceleration using OpenCL. OpenCLTemplate Marching Cubes is another resource made available to users who want to have easy access to GPU accelerated mathematical tools. The video below shows an example of the implementation:




1. Introduction


Marching Cubes is an algorithm used in a very wide range of applications, ranging from:

- Medical visualizations such as CT and MRI scan images;
- Special 3D effects and 3D modelling of metaballs or metasurfaces:
- Analysis of oil reservoirs in the oil and gas industry;
- Reconstitution of surfaces whose data has been acquired through seismic methods.

This OpenCL implementation allows the user to run and display objects using OpenGL.


2. Overview


CLMarchingCubes is implemented in a quite straightforward way to use: 

1. Create the function data you want to analyze and create isosurfaces from;
2. Create a new CLCalc.Programs.MarchingCubes object;
3. Set initial point and stepsizes using InitValues and Increment properties;
4. Calculate the isosurface using the CalcIsoSurface method;
5. Retrieve the results using the GetEdgeInfo method.

Using my ATI Radeon 4890 I have been able to draw isosurfaces in real time using a 51x51x51 grid without any problems. The calculations can be done much faster than that but preparing the data for OpenGL visualization takes some time.

The OpenCL/OpenGL proved not to be very effective because the clAcquireGLObjects seems to be a rather slow function. Optimized drivers in the future may make the interop as the best choice but right now using OpenCL to calculate the data and manually passing to OpenGL is faster.


3. OpenCLTemplate Marching Cubes commands description


3.1 Setting up


In order to create a new Marching Cubes stance, you need to create a 3-dimensional array of floats which may or may not contain the data you want to analyze. Either way, this vector will define the dimensions of the problem that will be instantiated by the Marching Cubes implementation. 

For example:


float[, ,] vals = new float[40, 50, 60];

CLCalc.CLPrograms.MarchingCubes MC = new CLCalc.CLPrograms.MarchingCubes(vals)


Will create a new MarchingCubes algorithm that will handle a grid with 40 X divisions, 50 Y divisions and 60 Z divisions. You may populate vals with the actual values you want to analyze from the beginning or choose to recreate the data to be analyzed later on.

The vertex coordinates are generated assuming that the step in each direction is 1.0 and that the initial point is [0,0,0]. You may change that by setting the InitValues and Increments properties.

The InitValues property takes a 3-dimensional float array which has the X, Y and Z initial coordinates.

The Increments property takes a 3-dimensional float array which has the X, Y and Z step sizes. For example:


MC.InitValues = new float[] { -2.0f, -2.0f, -2.0f };

MC.Increments = new float[] { 4.0f / (float)(maxX-1), 4.0f / (float)(maxY-1), 4.0f / (float)(maxZ-1) };


Will set the initial point to [-2, -2, -2] and set step values so that the final value is [2, 2, 2].

If you don't want to compute normal vectors you can set MC.ComputeNormals = false. From my experience, that doesn't improve performance very much.


3.2 Modifying point values


There are two ways to modify the grid point values. You can either pass the new values directly from the Host code, using a float[,,] variable that has the correct dimensions with which you initialized the Marching Cubes algorithm or manipulate the OpenCL memory buffer directly.

To pass new values from the host memory, use the SetFuncVals method:


float[,,] FuncValues;
//Populate FuncValues and create correct dimentions


To manipulate OpenCL memory buffer directly, manipulate the public MC.CLFuncVals variable in your own kernel.


3.3 Calculating the isosurface


In order to calculate the isosurface simply to invoke the CalcIsoSurface method:




The argument to this function is the desired value of the isosurface, i.e., the region where the function values would assume 1.0f value in this case.


If you pass the OpenCL variables using the custom constructor, you can use them to interoperate with OpenGL. They are the 3 vertex buffer objects to draw the OpenGl scene: vertex coordinates, vertex normals and element array index. Since this method has proven to be much slower, I will not go into details. The source code presented shows how to use the OpenCL OpenGL interoperation for this purpose.


3.4 Retrieving results


The method GetEdgeInfo is used to retrieve data. MarchingCubes returns the vertex coordinates of the edges, their normal vectors and the element array index containing the triangles:


List<int> ElemArray = new List<int>();
List<float> Coords = new List<float
List<float> Normals = new List<float

MC.GetEdgeInfo(out Coords, out Normals, out ElemArray);


The n-th vertex is so that its coordinates are: X - Coords[3*n], Y - Coords[3*n+1], Z - Coords[3*n+2]. Its normal vector's coordinates are X - Normals[3*n], Y - Normals[3*n+1], Z - Normals[3*n+2].

The element array contains information to build the triangles: the m-th triangle should link the vertexes ElemArray[3*m], ElemArray[3*m+1] and ElemArray[3*m+2] . This is exactly the way to build OpenGL vertex buffer objects and the code below shows this construction (assuming that you have already used glGenBuffers and stored the results in the bufs[] array):


List<int> ElemArray = new List<int>();
List<float> Coords = new List<float>();
List<float> Normals = new List<float>(); MC.GetEdgeInfo(out Coords, out Normals, out ElemArray);

//Copy buffer objs

EdgeCoords = Coords.ToArray();
EdgeNormals = Normals.ToArray();
ElementData = ElemArray.ToArray();

GL.BindBuffer(BufferTarget.ArrayBuffer, bufs[1]);
GL.BufferData(BufferTarget.ArrayBuffer, (IntPtr)(EdgeCoords.Length * sizeof(float)), EdgeCoords, BufferUsageHint.StreamDraw);

GL.BindBuffer(BufferTarget.ArrayBuffer, bufs[2]);
GL.BufferData(BufferTarget.ArrayBuffer, (IntPtr)(EdgeNormals.Length * sizeof(float)), EdgeNormals, BufferUsageHint.StreamDraw);

GL.BindBuffer(BufferTarget.ElementArrayBuffer, bufs[3]);
GL.BufferData(BufferTarget.ElementArrayBuffer, (IntPtr)(ElementData.Length * sizeof(int)), ElementData, BufferUsageHint.StreamDraw);


4. Conclusion


The Marching Cubes algorithm is a highly parallelizable algorithm that can be implemented in OpenCL and allows for real time visualization without problems of a 51x51x51 grid using a Radeon 4970 GPU. OpenCL/OpenGL interop proved not to be very effective because of the time-consuming clAcquireGLObject. It is expected that new drivers will fix this issue in the future.

This particular algorithm greatly benefits from multiple GPU cores since the number of work-items directly depends on the XYZ dimensions of the problem.


5. References


Lorensen, W.E. and Cline, H.E.,
Marching Cubes: a high resolution 3D surface reconstruction algorithm,
Computer Graphics, Vol. 21, No. 4, pp 163-169 (Proc. of SIGGRAPH), 1987.

Watt, A., and Watt, M.,
Advanced Animation and Rendering Techniques,
Addison-Wesley, 1992.

Polygonising a scalar field
http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/ access 19 jun 2010
Written by Paul Bourke May 1994   


Get sample source code.

Copyright © 2014 CMSoft. All Rights Reserved.
Joomla! is Free Software released under the GNU/GPL License.
Design by handy online shop & windows 7 forum