bugGNU Octave - Bugs: bug #55642, isosurface is slow

 
 

bug #55642: isosurface is slow

Submitter:  Guillaume <gyom>
Submitted:  Mon 04 Feb 2019 01:59:26 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Guillaume Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 23 Feb 2019 04:57:14 PM UTC, comment #8: 

That performance increase really is impressive!

I pushed the patch under your credentials here:
http://hg.savannah.gnu.org/hgweb/octave/rev/845aba35fca7

Good news that we got rid of the previous code :-)

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Fri 08 Feb 2019 02:39:56 PM UTC, comment #7: 

Thank you, Guillaume, for applying the principle to the helper function. Actually, last night I was also active. I attach a version of _unite_shared_vertices_.m that is now perhaps a bit more efficient because the faces are updated only once, and I have now a tolerance criterion that should be the tightest possible (it is now again vector-valued). And further, it is stable in the sense that the order of vertices as they are returned is the same as they were provided. This was not the case with my previous version, but it was the case with Markus' initial version, and this was a condition for his deletion of facevertexcdata to work. If you add the lines at the end


dum=true(size(lut));
dum(J)=false;
J=dum;


then you should be able to use this as a drop-in replacement for the previous version. More efficient would be to return J as it is and select the corresponding elements from facevertexcdata (then it would be more flexible if in the future this function is used for other purposes), or as you did it to pass it in and out and select from it in the loop.

Rik, sorry, I am doing this here on Octave 4.0.3 (greetings from Debian Stable) by pulling in the relevant .m-files from the web repository. I do not have a full mercurial set-up here so that I could prepare a changeset.

(file #46213)

Michael Leitner <mleitner>
Fri 08 Feb 2019 12:04:48 PM UTC, comment #6: 

I attach a version of Michael's patch where the changes now take place in _unite_shared_vertices_.m. As discussed here, the tolerance is now different than before and I don't know what is best; it passes all tests.

(file #46212)

Guillaume <gyom>
Thu 07 Feb 2019 04:49:53 PM UTC, comment #5: 

@Guillaume or Michael: Can you work the patch in to a complete changeset?  The speed-up is very impressive.

Rik <rik5>
Group administrator
Tue 05 Feb 2019 05:17:03 PM UTC, comment #4: 

Answers to your comments:
- yes, then it would be best to apply my patch just to _unite_shared_vertices_.m (and compute J according to its previous definition), because it indeed is the vectorization Markus asked for in the FIXME. Thus the speed improvement would apply also to reducepatch.
- I was initially thinking that marching_cubes works sequentially, thus that every cube can affect every other one. Then the eps of the largest vertex overall would indeed apply. But yes, seemingly it works in parallel, so only the largest vertex in the cube defines the achievable precision. So yes, it would become a vector, but I think the previous choice was too optimistic (when one vertex results very close to zero, its position is still affected by the finite precision of the other vertices of the cube it resulted from). But there is no one true solution --- you will always either keep returning vertices as distinct that with perfect arithmetic would be just one, or you will merge vertices that actually should be separated. So I would say, it is good enough.
- You want to drop all faces that have at least two equal vertices. And if neither f(1)==f(2) nor f(2)==f(3) nor f(3)==f(1), then necessarily all three are distinct. If you knew that the faces were sorted, then you would need only two comparisons (as in the last lines of _unite_shared_vertices_.m).

Michael Leitner <mleitner>
Tue 05 Feb 2019 04:52:59 PM UTC, comment #3: 

Thanks Michael, the speed improvement is impressive! To be fair, the speed issue was acknowledged by the authors in the FIXME note in _unite_shared_vertices_.m and it was probably not so important for small examples.

A few comments/questions:

  • _unite_shared_vertices_.m probably still has to be kept as it is also used by reducepatch.m (patch #8912).


  • Concerning tolerance, this was discussed in bug #46946, see Marco's comment #8. This would make the tolerance a vector instead of a scalar in the loop?


  • The one line I don't understand is this one. Why is it enough to only compare with a [2 3 1] reordering?



fvc.faces = fvc.faces(all (fvc.faces-fvc.faces(:,[2 3 1])!=0,2),:);


Guillaume <gyom>
Mon 04 Feb 2019 06:17:34 PM UTC, comment #2: 

This is a great opportunity to use the built-in profiler.


octave:5> profile on
octave:6> fv = isosurface (x, y, z, v, max (v(:))/2);
octave:7> profile off


If I then use profexplore I see


Top
  isosurface: 1 calls, 14.937 total, 0.001 self
    1) __unite_shared_vertices__: 1 calls, 14.856 total, 5.115 self
    2) __marching_cube__: 1 calls, 0.079 total, 0.027 self
    3) isosurface>__get_check_isosurface_args__: 1 calls, 0.001 total, 0.001 self
    4) nargin: 2 calls, 0.000 total, 0.000 self
    5) nargout: 2 calls, 0.000 total, 0.000 self
    6) isempty: 3 calls, 0.000 total, 0.000 self
    7) prefix !: 2 calls, 0.000 total, 0.000 self
    8) binary <: 1 calls, 0.000 total, 0.000 self
    9) binary >: 1 calls, 0.000 total, 0.000 self


So, Michael is correct that _marching_cube_ is way more efficient than _unite_shared_vertices_.

Michael's patch definitely decreases the time required.  This isn't my area of expertise so someone else should evaluate the mathematical logic (5 eps, for example as the tolerance).

Rik <rik5>
Group administrator
Mon 04 Feb 2019 05:41:30 PM UTC, comment #1: 

Yes, it is embarrassing if the main problem is handled by an efficient algorithm like marching cubes, and the post-processing is extremely inefficient -- to give numbers, on my computer the four lines with n = 64 take 15 seconds, while adding "noshare" to the invocation of isosurface cuts this down to 0.08 sec.

Your idea is obviously correct. Below I attach a patch that incorporates this idea. I also reworked the actual deletion of the vertices, which are now only seven lines (excluding comments) and thus probably do not warrant an outsourcing to _unite_shared_vertices_.m any more.

The example runs now in 0.10 sec. This could perhaps be decreased still by some 10 percent, along the lines of the FIXME in the code. But no need for C++.

(file #46178)

Michael Leitner <mleitner>
Mon 04 Feb 2019 01:59:26 PM UTC, original submission:  

isosurface becomes very slow when the size of the input data increases:


n = 64;
[x, y, z] = meshgrid (1:n, 1:n, 1:n);
v = (x-n/2).^2 + (y-n/2).^2 + (z-n/2).^2;
tic; fv = isosurface (x, y, z, v, max (v(:))/2); toc


Octave takes 11.34s while Matlab takes 0.04s.

If you increase n to 128, Octave takes 200s and Matlab 0.25s

And for n = 256, Matlab takes 2s while Octave was still running after an hour.

Most of the time is spent in _unite_shared_vertices_ where the loop over vertices is known to be slow (see bug #46946 and patch #8912).

One way would be to move _unite_shared_vertices_ into C++. Otherwise I wonder if one can implement 'unique with tolerance' from sorting the vertices:


v = sortrows (fvc.vertices);
d = diff (v, 1, 1);
all (abs (d) < tol)


I also notice that the help text of isosurface still mentions the use of 'unique':


## If given the string input argument @qcode{"noshare"}, vertices may be
## returned multiple times for different faces.  The default behavior is to
## eliminate vertices shared by adjacent faces with @code{unique} which may be
## time consuming.


Guillaume <gyom>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #46212:  isosurface.patch added by gyom (5KiB - text/x-patch)
file #46178:  patch added by mleitner (1KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mleitner (Updated the item)
  • -email is unavailable- added by gyom (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-26 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2019-02-23 mmuetzel StatusConfirmed Ready For Test
    2019-02-08 mleitner Attached File- Added __unite_shared_vertices__n.m, #46213
    2019-02-08 gyom Attached File- Added isosurface.patch, #46212
    2019-02-04 rik5 StatusNone Confirmed
    2019-02-04 mleitner Attached File- Added patch, #46178

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code