<div dir="ltr">FYI- vtkFlyingEdges2D and vtkFlyingEdges3D have been merged into VTK proper. Thanks to help from Rob, Sujin, Bill L., Utkarsh, Sankhesh and many others...<div><br></div><div>We will be presenting a paper at <a href="http://www.ldav.org/">LDAV 2015</a> in October. (Ask me for a preview copy of you are interested.)</div><div><br></div><div>Benchmarking indicates this is a very fast isocontouring algorithm, especially if you have multiple cores your machine (without generating seams or requiring point merging...build VTK with TBB using vtkSMPTools). So far we've run it on 2048^3 volumes with up to 160 threads, it runs great. Like any new algorithm I'm sure there are warts and implementation issues, so I'd enjoy hearing feedback. Maybe someday we'll slip it under the hood of vtkContourFilter, but I want to make sure it is rock solid first.</div><div><br></div><div>If you are really ambitious, there is room for further improvement in the implementation including parallelizing the prefix_sum (third pass), and in the fourth pass (decoupling the topology and geometry generation). I can tell you more if you are interested.</div><div><br></div><div>Best,</div><div>W</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Jan 10, 2015 at 10:04 AM, Will Schroeder <span dir="ltr"><<a href="mailto:will.schroeder@kitware.com" target="_blank">will.schroeder@kitware.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Over the holidays I gave myself the gift of designing and implementing a new contouring algorithm :-) To my utter delight and amazement I stumbled upon a novel approach that is readily parallelizable and yet wicked fast in serial mode as well. The algorithm (vtkFlyingEdges2D/vtkFlyingEdges3D) has been pushed into gerrit for review if anyone is interested (FlyingEdges branch). Consider this a work in progress, we are benchmarking and working through parallel implementation details now. We are planning a paper submission in the near future.<div><br></div><div>TLDR Details</div><div>+ As many of you know the VTK community is taking aim at seriously reworking the system to better support parallel computing. Much of this work is being led by Berk Geveci and his team. For example, Berk introduced the vtkSMPTools into VTK about a year ago, and we are actively working with some big name chip manufacturers, DOE Labs, research organizations, and other customers/collaborators to get this ambitious undertaking done. Hopefully individuals in this VTK community will also lend their significant expertise as well.</div><div>+ Under prodding by Berk, and inspired by Patrick O'Leary's "born parallel" slogan to rethinking systems and algorithms for parallel computing, the motivation was to take one of the most important visualization algorithms (isocontouring) and see if we could redesign it for emerging parallel approaches.</div><div>+ The venerable Marching Cubes algorithm introduced the key idea that the topology of an isosurface within a voxel can be captured with an index into a simple, finite case table. Yet the implementation of the algorithm can be slower than you might expect because, for example, on interior voxels, voxel edges can be processed up to four times, and vertices visited eight times (ignoring gradient calculations). Also the output mesh is dynamically created: a priori the number of output points and triangles is unknown, which means that memory resizing (realloc's) are required. Finally, the algorithm is often implemented with a "locator" or hash table to prevent the creation of duplicate points on shared voxel edges (the so-called point merging problem). Not good for parallel computing.</div><div>+ Ken Martin's extremely fast Synchronized Templates algorithm introduced the notion of a "voxel axes" or "edge triad" which consists of the three x-y-z voxel axes located at the voxel origin. This triad is moved through the volume in such a way as to intersect each voxel edge only one time (with special treatment on the +x+y+z volume boundaries). (An aside: Ken wrote this like 15 years ago, it's been the VTK standard ever since. Initially we filed a patent application [since abandoned] but that was before we realized how stupid it is for an open source company to hold patents. Live and learn :-))</div><div>+ Flying Edges builds on these concepts and introduces some novel ideas. It is a 3-pass algorithm: the first two passes gather information about the volume including "intersections" along x-edges, and the number of output primitives and points. The final pass generates output into preallocated arrays. Some of the key ideas include:</div><div>-- each pass traverses x-rows in the volume completely independently meaning that it is readily parallelizable.</div><div>-- the case table is based on the four voxel x-edges (which can be easily transformed to and from a vertex-based MC table). Because the cases are edge based, the computation of cases is parallel separable.</div><div>-- The novel concept of "computational trimming" is used to avoid massive amounts of work. Basically the fact that isosurfaces are topologically smooth is combined with the results of the first pass (building x-voxel-edge cases) to reason about where the isosurface can exist in the volume. This very simple topological reasoning enables the algorithm to skip portions of x-rows, entire x-rows, and even x-y slices very rapidly. </div><div>-- Some simple edge metadata is gathered for each x-row volume edge which can be used to allocate and partition memory for separate threads to write into without contention (i.e., no mutexes or locks). This metadata also enables the algorithm to process the voxels along each x-row in such a way as to synchronously increment the point ids and triangle ids without requiring point merging or hashing. In fact, topology generation (i.e., creating triangles) is completely independent from point/geometry generation (i.e., interpolating point coordinates and attributes along edges).</div><div>+ Early results show a 2-5x speedup over Synchronized Templates (>10x over vtkImageMarchingCubes) when running in serial. However, this is very data dependent and I suspect compiler dependent since there is a lot of C++ templating and in-lining going on (make sure to build optimized/release mode). Anyway I am very excited because once the parallel loops are functioning properly we should see much more performance improvement. The scalability may be hampered however because data movement over the memory bus may slow things up (computations are relatively modest compared to the total data processed). We'll be benchmarking performance over the next couple of months to see where we really are and of course we'll update this community at that time.</div><div>+ The memory overhead of the algorithm is roughly 2-bits per vertex. However as implemented now I'm avoiding bit packing and just use a byte (unsigned char) per vertex.</div><div>+ Other interesting tidbits:</div><div>-- The basic algorithm is extremely simple to implement. Most of the complications come from dealing with volume boundaries. It would be possible to significantly reduce code complexity by either padding the +x+y+z volume boundaries by a layer of voxels, or alternatively not processing the outermost layer of voxels.</div><div>-- Because I am lazy, the algorithm reuses the MC case table as enumerated in VTK (vtkMarchingCubesTriangleCases.h). However at instantiation the MC table is converted to a Flying Edges edge-based case table. Also note that the MC table as originally defined is for a left-handed coordinate system (see the original paper). So reordering of triangles is required to make sure that normals and triangle ordering is consistent.</div><div>-- On our radar is potential vectorization of operation like interpolation across voxel edges. As you know some of the big hardware manufacturers are really pushing this now and we'd like to learn how to best take advantage of this emerging resource.</div><div>-- Probably the biggest downside of this algorithm is that is can produce degenerate triangles (i.e., zero-area triangles in 3D, zero-length line segments in 2D). The major reason is that when preallocating output memory, we do so based on topological evaluation. However, during actual point generation in degenerate cases (a degenerate case occurs when one of more vertex scalar values == isocontour value) multiple intersection points can occur at the vertex producing degeneracies. Since discarding a degenerate triangle would mess up the output preallocation we can't just throw away degenerate primitives. However I have already thought through a solution to this problem but it is based on (spoiler alert) a much more complex case table that enumerates 3-values per vertex: (<,==,>) when compared to the isocontour value. This is a novel idea as well, mostly the academic literature completely ignores the reality of degeneracies. Maybe next holiday season we'll slip it into the algorithm.....</div><div><br></div><div>Okay I've got to get back to the pointy-haired management stuff. But please feel free to provide feedback, or offers of collaboration :-) In particular we'd love support, either in funding or brainpower, to do even more cool parallel processing stuff. Hoping to hear from you, and I hope your coming year is as exciting as mine is looking to be!</div><div><br></div><div>Best,</div><div>W<span class="HOEnZb"><font color="#888888"><br clear="all"><div><br></div>-- <br><div>William J. Schroeder, PhD<br>Kitware, Inc.<br>28 Corporate Drive<br>Clifton Park, NY 12065<br><a href="mailto:will.schroeder@kitware.com" target="_blank">will.schroeder@kitware.com</a><br><a href="http://www.kitware.com" target="_blank">http://www.kitware.com</a><br><a href="tel:%28518%29%20881-4902" value="+15188814902" target="_blank">(518) 881-4902</a></div>
</font></span></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">William J. Schroeder, PhD<br>Kitware, Inc.<br>28 Corporate Drive<br>Clifton Park, NY 12065<br><a href="mailto:will.schroeder@kitware.com" target="_blank">will.schroeder@kitware.com</a><br><a href="http://www.kitware.com" target="_blank">http://www.kitware.com</a><br>(518) 881-4902</div>
</div>