This is nice. Could you add some documentation to your header file? And maybe add a header similar to other libigl files.

Also, I don't see why you need to implement the algorithm with an unordered map. I think it would be more efficient to generate duplicated vertices first, and then remap the duplicates at the end. An efficient way to achieve this would be to store the map (e1,e2),v1 in a std::vector<std::pair<std::pair<int,int>, int>, then sort the vector and count the duplicates.

Besides, it looks like you are only taking the edge midpoints as vertices positions, but shouldn't you be interpolating its positions as a weighted combination of the edge endpoints, based on the isovalues evaluated at the tet corners?

Sorting is silly when we can do the merge in linear time (and would probably be slower than what I have now). There's a function igl::remove_duplicates that dedups in linear time. One thing worth noting is that this approach will be much slower in debug mode (I tried it). If we do intend to parallelize, then it makes sense.

Currently igl::marching_cubes doesn't interpolate so I didn't want to have a different behavior. I do agree it's kind of strange. Maybe I'll make it a flag.

Sorting is not silly, in practice it is likely to be much faster than using a std::unordered_map (the asymptotic complexity is not everything, the hidden constant and the cache coherence also play a huge part).

If you look at the implementation of igl::remove_duplicates it has quadratic complexity, which is pretty bad. If you are talking about igl::remove_duplicate_vertices, what this function is doing is exactly a sort+unique on the rows of the matrix of vertices positions (so again, n log(n) complexity). But using one of those functions to remove duplicates is silly, this they both rely on the geometric position of the vertices, and might be prone to roundoff errors/depend on your choice of epsilon. It is much more sensible to do the sort+unique yourself (based on id of the edge endpoints), and it really is but 3 lines of code.

Finally, if you look at the code of igl::marching_cubes, you will find that it also interpolates the vertices positions, unless I have missed something else?

Just to give you some numbers, I implemented my version with the sort and tested it on the shared/cube.obj tetrahedralized with pq1.414a0.000001, which produced 1920674 tets. On this example my version with the sort was still 2x faster (30.328ms) than your version with the unordered_map (65.004ms), both compiled in release mode.

Ah, I didn't realize remove_duplicates was quadratic. I guess that makes sense. We can't quite use remove_duplicate_vertices here since we need to fix the face indexes as well and remove_duplicate_vertices doesn't do that. I guess I misremembered marching_cubes not interpolating. I'm about to send a fix now.

My original comment meant to discuss sorting vs. hashing to do deduplication after the initial pass (not hashing as you go vs. sorting at the end). The point was that you can deduplicate in linear time at the end.

To determine which was was best, I implemented 3 separate tests. The first has my original code (with a bug fix discussed below), The second does sorting to dedplicate the vertices and uses the sorted array to fix the face map. Finally, the third replaces the sort step with a step that uses hashing to do deduplication. They all use about the same amount of extra memory (sorting still requires we store a lookup to fix the faces).

While implementing sorting, I found a bug in the original code I submitted where this line: for(int e = 0; mt_cell_lookup[key][e] != -1; e++) loops past the end. The fix was changing it to: int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++). Sorting could not have worked with this so I'm glad I caught it. This bug was also causing my code to insert into the hash table about 2x as many times as it needed to.

I ran your benchmark above (using an isovalue of 0.5) over 10 runs of each implementation in Release mode and I get the following results:

MT with SORT took (average over 10 runs) 0.0975869s
MT with HASH took (average over 10 runs) 0.161265s
MT with HASH2 took (average over 10 runs) 0.0712481s

I concede that my original code is still slower than the sort-based implementation for this benchmark but with the bugfix, it's not quite 2x. Doing hash based deduplication does speed things and has better asymptotic performance so I think it's a better approach. If you have a way to do sorting without allocating vmap to fix the face topology then sorting could use less memory and that might be worth it.

In any case, here are the three functions I implemented: https://pastebin.com/WmbeTFuG. I'm going to send a PR with the third. Also here is the code to do the benchmarking: https://pastebin.com/X499PRUP.

EDIT: I refactored the code before putting it on pastebin and broke it. I updated this comment with the fixed version (old version here: https://pastebin.com/281CpDEK)

Right. I did fix the issue with e < 4 in my implementation with std::sort, but not in your code I was comparing against. This brought the total time of your original code down to 47.8476ms on this example (average over 10 runs). But you need to do the check before accessing mt_cell_lookup[key][e], otherwise you are still accessing a value out-of-range.

Your 3rd version took in average 30.3694ms on my computer, which is basically the same as my version with the sort (27.7268ms if I average over 10 runs, so the allocation of the output is amortized). In your benchmark code it seems you are counting the initial allocation of the output matrices only for the std::sort version, so it may bias things a bit.

Anyway, my case stands that both are comparable. The std::sort version a bit faster in my case, but it is certainly not slower. Your version with the hash also assumes that indices can be packed into a 64bit integer, and you are forcing the size of the types you are working with. Hash maps may also have performances that vary depending on the density of the output, how many collisions you get, etc. On the other hand the std::sort version will be more predictable.

My point in all this, is that there is no advantage in using a data structure that is complicated (hash map), for something that could be achieved with a std::vector and a std::sort. This is more of a general comment that you don't need to use a sledge-hammer to swat a fly. In regard to this particular implementation of marching test I honestly couldn't care less which version is used. I just wanted to carry my point across. In life you can often get away with a simple std::vector, and don't need hash maps, and it certainly doesn't hurt to think that way most of the time.

Regarding your specific commit you still need to change the following to adhere to the libigl coding style:

• Don't put helpers in the global namespace, use lambdas instead.
• Put braces on a new line.
• And I think the template instantiation should be a single (long) line.

I re-ran a new benchmark where each method gets a new vector (I didn't realize Eigen::resize() would no-op when the array was already the right size): https://pastebin.com/GS3XgKj8. Here it is with 100 runs:

MT with SORT took (average over 100 runs) 0.0840349s
MT with HASH took (average over 100 runs) 0.131921s
MT with HASH2 took (average over 100 runs) 0.0578004s

I think the point is that there is an advantage to using an unordered_map. This is a pretty small example (30ms is nothing). For the sake of argument, let's say both methods are equal for the small example. The hashing method has a lower asymptotic complexity and uses about the same amount of memory. For bigger examples, it will perform better.

If we never cared about big examples in the first placc, then we should go with the simplest piece of code (IMO that was the first thing I pushed). If we're separating deduplication out, we presumably want the option to parallelize so then we must care about large examples! In that case, we should go with a lower asymptotic complexity solution.

Furthermore, the hashing solution is less code and, (subjectively IMO) easier to read. I also don't think that a hash map of ints is a particularly complicated data structure (which is also my subjective opinion I guess).

Finally, it's also worth noting that there are trivial hash combines (see the one I added in igl::marching_cubes which is copypasta from boost::hash_combine) which can get rid of the bit-width requirement. But honesly, if your compiler doesn't have 2*sizeof(int32_t) != sizeof(int64_t), then you have much bigger problems (even 32 bit compilers support 64 bit int64_t).

I'll make sure I follow libIGL coding style. I'll fix my editor to code in that style to begin with.

Hashing is a complicated business, and there are a gazillion of papers about how to implement good hashing functions. You are free to think it is a magic black box with constant time complexity, but I'm just advocating that it doesn't hurt to think with std::vector from time to time.

Btw, since we are nitpicking, your implementation with the std::sort is allocating too many elements in vector<Eigen::RowVector3d> vertices, compared to the hash map version (which allocated only 1 row per unique edge). Three double is a lot of memory, and this is were the difference in timings come from.

If you fix this, then the std::sort version stays faster that your second hash map implementation. I took your last comparison code, and tried to see how it grows as I divide by 2 the volume in tetgen, and here the result that I got:

#tets:1920674
MT with SORT took (average over 100 runs) 0.0596599s
MT with HASH took (average over 100 runs) 0.133467s
MT with HASH2 took (average over 100 runs) 0.0660952s

#tets:3804863
MT with SORT took (average over 100 runs) 0.114474s
MT with HASH took (average over 100 runs) 0.278299s
MT with HASH2 took (average over 100 runs) 0.135596s

#tets:18793141
MT with SORT took (average over 100 runs) 0.782337s
MT with HASH took (average over 100 runs) 1.29649s
MT with HASH2 took (average over 100 runs) 0.933668s

To make sure we are on the same page, here is my implementation.

Regarding your PR, can you add a header to the files marching_tets.cpp|h (but not in the tutorial) like this one?

// This file is part of libigl, a simple c++ geometry processing library.
// 
// Copyright (C) 2018 Foo Bar <foobar@gmail.com>
// 
// This Source Code Form is subject to the terms of the Mozilla Public License 
// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
// obtain one at http://mozilla.org/MPL/2.0/.

And also fix the for loop to check for e < 4 && mt_cell_lookup[key][e] != -1, to avoid the out-of-bound read.

PS: A few minor tips:

• You can write emap.emplace(key, num_unique); instead of emap.insert(make_pair(key, num_unique));, and you can write edge_table.emplace_back(v1_idx, v2_idx); instead of edge_table.push_back(make_pair(v1_idx, v2_idx));
• You probably don't need to specify the return type -> int64_t in the lambda function.
• The static check that sizeof(std::int64_t) == 2*sizeof(std::int32_t) is probably unnecessary (according to the standard, int32_t is supposed to be exactly 32 bits, and int64_t is supposed to be exactly 64 bits).

Damn, I was trying to sneak in a fix before you noticed :P

I just pushed the change. Thanks for the feedback!

Ahah, too slow =)

Thanks for the changes. I'll merge once it's done building!

@alecjacobson can we merge?

@jdumas Since it appears we're in a battle of nits (and honestly it's fun to optimize this code to death). I updated my hashing code a bit. First off, I used the optimizations in your sort version (don't interpolate until deduplication, use min/max instead of branching, etc...). I also removed my std::vector of vertices.

I also realized I wasn't reserving memory for my hash table even though I know exactly how many elements I'm going to insert so that means I'm rehashing everything a bunch of times. Also for some stupid reason the default max load factor of std map is 1.0 (I know that stupid reason is memory but it should really be 0.5 for anything serious). Anyway here are my new results (original hash omitted because we all agree that it's slow) and a pastebin with benchmark code:

#tets 1920674
MT with SORT took (average over 100 runs) 0.0551038s
MT with HASH2 took (average over 100 runs) 0.0388276s

#tets 3804863
MT with SORT took (average over 100 runs) 0.11682s
MT with HASH2 took (average over 100 runs) 0.0909922s

#tets 7564920
MT with SORT took (average over 100 runs) 0.281586s
MT with HASH2 took (average over 100 runs) 0.23539s

#tets 15060886
MT with SORT took (average over 100 runs) 0.620236s
MT with HASH2 took (average over 100 runs) 0.523281s

I'll amend my original commit with these changes.

@alecjacobson which function are you referring to? I am pretty sure that we do not currently have an implementation of marching tetrahedra (which is a copyleft-free replacement of marching cubes!).

Oh ok. Do you mind if we merge this now, and eventually replace the implementation with yours (if there are advantages) when you will get your laptop back? The tutorial code will be useful in any case. We need this function for a couple of projects and I would like to merge it asap.

We can definitely wait a couple more days, we are not in a rush.

Finally got my laptop back. I have a duplicate function (actually it was in my dev branch all along).

I'll see which is better/faster and report back.

Meanwhile, the tutorial entry has a lot of issues. @fwilliams could you simply remove the tutorial entry from this PR and keep just the new function?

Thanks @alecjacobson for checking back. Could you comment on the problems with the tutorial? One thing I can see is the visualize_tet_wireframe function inserts duplicate edges. I have a fixed version on my local branch if this is an issue. I'll gladly remove the tutorial but I'd love to know what the problems are.

Sure. But let's split this into two PRs: 1) performance tuning for isosurface extracting from tets and 2) a tutorial entry for said extraction.

Regarding 1), I now realize that the master branch already had a function that does isosurface extraction from a tet-mesh igl/slice_tets.h. This PR is still cool because it's about 4x faster (on my machine). However, the current slice_tets has much richer output (not just the raw mesh). So, this PR should be accepted as a performance optimization on slice_tets but we should maintain the features we already have (doesn't look too hard).

Regarding 2), the issues I immediately saw were:

1. reimplementation of igl/edges in the tutorial code (tutorial code should be simple and short; new code shouldn't reimplement code we already have),
2. calls tetgen rather than reading an existing .mesh (again, complicates tutorial code; makes TetGen a dependency on this entry; makes this example slower)

(btw, thanks for your patience guys. It's way easier to straighten out duplicate code/functions beforehand than months later)>

Ahah, so much for not reinventing the wheel :p (the slice_tets implementations seems rather long compared to @fwilliams's version).

Btw, is there any function in libigl for subdividing tets? That could be a way to produce a surface with some resolution, without having to load a gigantic .mesh file for the tutorial.

No, and wouldn't be a good idea, even if we did. Just load shared/bunny.mesh.

@alecjacobson did you manage to sort out the code duplication? Can we help you out with this? I would like to integrate this function soon since we need it for a few projects.

I had most of it done locally. Here's the version with feature parity to slice_tets which is much faster and includes a bunch of useful overloads.

@alecjacobson Does everything look good to merge?

This looks great. I need to find time to check that it runs smoothly on some examples.

Alternatively, you could add a unit test for slice_tets to https://github.com/libigl/libigl-unit-tests and test that this PR doesn't break it.

Bump.

Any news on this?

Thanks for the bump. I took a peek a while back and convinced myself this is ready to merge. I'm trying to find some time to wrangle function names a bit. I will merge soon.