Source Code for Module examples.voxelization

  1  #!/usr/bin/env python 
  2   
  3  #! This file is a literate Python program. You can compile the documentation 
  4  #! using mylit (http://pypi.python.org/pypi/mylit/). 
  5  ## title = "glitter Example: Voxelization" 
  6  ## stylesheet = "pygments_style.css" 
  7   
  8  # <h1><i>glitter</i> Example: Voxelization</h1> 
  9   
 10  # <h2>Summary</h2> 
 11   
 12  # This program will create an invisible OpenGL context to voxelize a mesh given 
 13  # as an <a href="www.hdfgroup.org/HDF5/">HDF5</a> data file. The resulting 
 14  # volume is written to a file. 
 15   
 16  #! TODO describe voxelization algorithm 
 17  #! TODO include spreadbits.py so the data can be viewed in the volume viewer 
 18   
 19  # <h2>Front matter</h2> 
 20   
 21  # <h3>Module docstring</h3> 
 22   
 23  # The module docstring is used as a description of this example in the 
 24  # generated documentation: 
 25  """Voxelize a mesh. 
 26   
 27  @author: Stephan Wenger 
 28  @date: 2012-02-29 
 29  """ 
 30   
 31  # We can usually import classes and functions contained in <i>glitter</i> 
 32  # submodules directly from glitter: 
 33  from glitter import VertexArray, TextureArray2D, uint32, Framebuffer, State, ShaderProgram 
 34   
 35  # <h2>Shaders</h2> 
 36   
 37  # The vertex shader transforms the position as usual. 
 38  vertex_code = """ 
 39  #version 400 core 
 40   
 41  layout(location=0) in vec4 in_position; 
 42  out vec4 ex_position; 
 43   
 44  void main() { 
 45      gl_Position = ex_position = in_position; 
 46  } 
 47  """ 
 48   
 49  # The fragment code will perform the actual work. It will quantize each 
 50  # fragment into a depth layer, compute a bitmask corresponding to that depth, 
 51  # and write the bit mask into the color attachments. The number of color 
 52  # attachments is only known at runtime, which is why a <code>%d</code> 
 53  # placeholder will have to be filled in just before the shader is compiled. 
 54  fragment_code = """ 
 55  #version 400 core 
 56   
 57  in vec4 ex_position; 
 58  uniform bool solid; 
 59  layout(location=0) out uvec4 fragmentColor[%d]; 
 60   
 61  uvec4 voxelize(uint z, uint offset, uint fill, uint base) { 
 62      if (z >= (offset + 1u) * 128u) return uvec4(0u); 
 63      if (z < offset * 128u) return uvec4(fill); 
 64      z -= offset * 128u; 
 65      return uvec4( 
 66          z < 128u ? (z >= 96u ? base << (z - 96u) : fill) : 0u, 
 67          z <  96u ? (z >= 64u ? base << (z - 64u) : fill) : 0u, 
 68          z <  64u ? (z >= 32u ? base << (z - 32u) : fill) : 0u, 
 69          z <  32u ?             base << (z -  0u)         : 0u 
 70      ); 
 71  } 
 72   
 73  void main() { // if solid, switch on all bits >= z; else, switch on bit == z 
 74      uint z = uint(round(128.0 * float(fragmentColor.length()) * 0.5 * (ex_position.z + 1.0))); 
 75      for (int i = 0; i < fragmentColor.length(); ++i) 
 76          fragmentColor[i] = voxelize(z, i, solid ? 0xffffffffu : 0u, solid ? 0xffffffffu : 1u); 
 77  } 
 78  """ 
 79   
 80  # <h2>Voxelization</h2> 
 81   
 82  # The function <code>voxelize()</code> performs solid or boundary voxelization 
 83  # of a mesh into a volume of given size. 
 84 -def voxelize(mesh, size, solid=True): 
 85      """Voxelize a mesh into a cube of length C{size}. 
 86   
 87      @param mesh: The mesh to voxelize. 
 88      @type mesh: L{VertexArray} 
 89      @param size: The length of the resulting cube. Should be a multiple of 128. 
 90      @type size: C{int} 
 91      @param solid: Whether to perform solid voxelization instead of boundary voxelization. 
 92      @type solid: C{bool} 
 93   
 94      @attention: The result of boundary voxelization is likely to have holes. To 
 95      obtain a watertight volume, you can voxelize the volume from different 
 96      directions and combine the results, or use the gradient of a solid volume. 
 97      """ 
 98   
 99      # An array of 2D integer textures will hold the output volume encoded in 
100      # its bits: 
101      volume = TextureArray2D(shape=(size // 128, size, size, 4), dtype=uint32) 
102   
103      # We create a framebuffer with all layers of the volume attached and 
104      # directly bind it: 
105      with Framebuffer(*[volume[i] for i in range(len(volume))]) as fbo: 
106          fbo.clear() 
107          # Depending on the voxelization mode, we set the logic op mode. We then 
108          # create and activate a shader program and draw the mesh. 
109          with State(logic_op_mode="XOR" if solid else "OR", color_logic_op=True): 
110              with ShaderProgram(vertex=vertex_code, fragment=fragment_code % len(volume), solid=solid): 
111                  mesh.draw() 
112   
113      # Finally, the texture array is returned. The client code can decide to 
114      # download the data or to process it further on the GPU. 
115      return volume 
116   
117  # <h2>Initialization and main loop</h2> 
118   
119  # Finally, if this program is being run from the command line, we set up all 
120  # the previously mentioned objects and start rendering. The program will expect 
121  # two command line parameters: the name of an HDF5 input file, and the name of 
122  # an HDF5 output file. An optional third parameter specifies the size of the 
123  # resulting volume. 
124  if __name__ == "__main__": 
125      # We need to read a mesh filename from <code>sys.argv</code>, so import 
126      # <code>sys</code>. 
127      import sys 
128   
129      # We assume the mesh is stored in a <a 
130      # href="http://www.hdfgroup.org/HDF5/">HDF5</a> file, so import <a 
131      # href="h5py.alfven.org"><code>h5py</code></a>. 
132      import h5py 
133   
134      # We open the HDF5 file specified on the command line for reading, extract 
135      # the vertices and indices, and store them in a vertex array: 
136      with h5py.File(sys.argv[1]) as f: 
137          mesh = VertexArray(f["vertices"], elements=f["indices"]) 
138   
139      # Now the voxelization is performed and the result is returned as an array 
140      # of 2D integer textures. 
141      volume = voxelize(mesh, int(sys.argv[3]) if len(sys.argv) > 3 else 128) 
142   
143      # Finally, we open the output file, download the data, and store it in a 
144      # dataset. The dataset can typically be compressed a lot, which is why we 
145      # enable LZF compression. However, this may take a long time, so you can 
146      # disable it if you are not short on disk space. 
147      with h5py.File(sys.argv[2], "w") as f: 
148          f.create_dataset("data", data=volume.data, compression="lzf") 
149   
150      # Note how we did not need to create an OpenGL context explicitly. The 
151      # first object that needs an OpenGL context (in this case the 
152      # <code>VertexArray</code>) will create an invisible context as soon as it 
153      # is initialized. This is very convenient; on the other hand, it means that 
154      # if you create a context manually (e.g., because you need a visible 
155      # window), you have to initialize it before any other objects. Otherwise, 
156      # these objects may end up in the wrong context, and sharing data across 
157      # context is a problem of its own. 
158