RFE: optimize surface zone calculations

[Tristan Croll]

(Not urgent, but it would be nice in many ways)

I've been spending quite a bit of time on parallelising Clipper's structure factor calculations for the live-updating crystallographic maps, and have got it down to the point where for the demo model I get about 2 map updates a second on my Xeon workstation. The map calculations themselves and contouring of the portion "seen" by ChimeraX all happen in separate C++ threads, so the time overhead at those steps are fairly minimal (a few ms, mostly in the copying of ChimeraX atoms to Clipper before launching the map recalculation thread).

Pretty cool to watch (and practically useful to boot) - but now the rate-limiting step is VolumeSurface._set_surface() (or, really, Drawing.set_geometry()): if I have the whole demo model covered and simulating with the standard four surfaces shown (standard map, sharpened map, and +/- difference map), then the grand total cost of _set_surface() for all four surfaces is about 100ms - 44ms each for the standard and sharpened maps, and ~6ms each for the difference maps (simpler; fewer triangles). It's obviously substantially better when smaller regions are displayed, but still enough to introduce a noticeable "hitch" to the display.

I can imagine that there are many cases where it would be nice to improve on this (essentially anything to do with animating volumes, for a start).

[Tom Goddard]

The VolumeSurface._set_surface() routine does almost nothing -- sets a a handful of variables. So I'm not sure what you are actually timing. I guess the times you cite are for the marching cubes contour surface calculation (VolumeSurface.update_surface) that creates the surface triangulation from the density map and a threshold level. That code already has multi-thread support, see the "use_thread" argument of VolumeSurface.set_level(), although it is not used by default because it can delay surface calculations for many frames which can cause problems if subsequent commands (e.g. calculating surfaces properties) incorrectly assume the surface has already been calculated.

[tic20@…]

I put timing calls at each step of update_surface(), and it was _set_surface() that by far took the bulk of the time. Looking at it, this is where Drawing.set_geometry() gets called, which in turn triggers its redraw callback - I assume that this is where the cost is. I know it’s very fast for lots of instances of simple geometry, but when it’s a single instance of a complex geometry it starts to slow down.

 
 

 

[tic20@…]

Ah! My mistake - after some more profiling the culprit is actually the 
rezoning callback.

On 2019-05-22 06:08, ChimeraX wrote:

[Tristan Croll]

Almost all the time here is going in to find_close_points() - which isn't hugely surprising, I suppose... for this test case there are 3353 atoms and 64436 vertices to get through. From some preliminary timings it looks like a more brute-force approach may be quite a bit faster here. If this volume is v and its surface is s, then:

%timeit v.interpolated_values(s.vertices, out_of_bounds_list=True)
1.67 ms ± 144 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

... so precalculating a binary mask wrapped in a Volume object should be fairly straightforward (and trivially parallelisable if necessary), and would save a whole lot of time (particularly in my case where the same mask will be applied to all volumes). Working on a test implementation at the moment.

}}}

[Tom Goddard]

I see. Yes, Drawing.set_geometry() does do the mask and per-vertex color updating. Using the ChimeraX zone calculation every time the map changes, when the zone always stays the same is going to slow things down, and precomputing the zone mask seems like a reasonable way to optimize that.

[Tristan Croll]

I have some actual timings on this now. If I create a Volume of datatype numpy.uint8 with the same resolution as the map, and just big enough to cover the coordinates + mask, then for the same test case as above (1658 atoms, 4A cutoff radius) creating the mask takes 35 ms, and interpolating the values at all vertex positions in the target surface takes 1.9ms. The combination is still faster than a single call to find_close_points(), and in most cases the one mask will be re-used across multiple surfaces.

There's lots of scope to speed up the mask calculation itself. I wrote this version to use arbitrary cell angles, but that's pretty silly - limit it to strictly cubic voxels and a huge number of transform operations will go away.

... Actually, I just tried that. Brings the mask recalculation down to 10.6ms. Could still be optimised further by setting its resolution according to the desired cutoff, rather than the resolution of the target map.

That should make it just about fast enough so the mask *can* be updated with every coordinate update (I'd actually like to make that happen, so that the map responds naturally when atoms move a long way).

[Tristan Croll]

4A mask, 1.33A grid

[Tristan Croll]

... and a little more tinkering gets it down to 2.7ms for a 4A map (assuming the volume extent stays the same so it just has to be refilled - obviously slightly longer if a new box needs to be created). I'll tidy it up a bit more tomorrow - happy to pass it to you after that if you want it for more general use.

[Tristan Croll]

OK, I have a working implementation of this approach that looks about the same as before, and gets the real-world time for updating each surface in my test case down from 45ms to ~5ms, by re-using the one mask across multiple surfaces. Good enough for my purposes for now. Feel free to close this if it's not a priority for you.

[Tristan Croll]

Incidentally, perhaps it's just late in the day but getting a mask to apply on first drawing of a Volume seems fiendishly tricky. A minor issue that I'm sure I'll work out, but thought I'd mention it anyway.

[goddard@…]

I'm not clear which approach you are using to mask.  You could be masking the surface, basically just hiding some vertices, that is what the volume zone code was doing.  Or you could actually be masking the volume data itself, setting map values to zero outside the region of interest.  I can't comment on why it would be hard get it to be used when the volume is first shown since I don't know what you are doing.

I reassigned the ticket to you, so you should close it when you are satisfied it is resolved.

[tic20@…]

I’m masking the surface as per the existing Zone code. It just seems that on the first time the Volume is drawn, the triangle mask is either lost or ignored. I’ve tried delaying application of the mask until the following “frame drawn” trigger, but that still has no effect. Yet everything works fine from that point on. Can’t see what I’m doing wrong, but I’m sure I’ll figure it out.
 

[goddard@…]

Strange.  The Drawing.set_geometry() code runs the auto mask updating, so it should mask the first time assuming you setup the masking before the surface is computed.

[tic20@…]

Will dig into it properly tomorrow. 
 

 

[Tristan Croll]

Ah - I see the problem now. Since the surface calculations are threaded, the initial creation of the ZoneMask object happens before the surface actually has any vertices or triangles, so set_surface_mask() returns early and surface.auto_remask_triangles is never set. Note that this will also affect your code in surface/zone.py if you use threaded surface calculations. The below amendment seems to work fine.

By the way, the Clipper plugin is now out in a public GitHub repository at ​https://github.com/tristanic/chimerax-clipper - since I designated it as LGPL it's about time I did that.

class ZoneMask(State):
    def __init__(self, surface, points, distance, max_components):
        self.surface = surface
        self.points = points
        self.distance = distance
        self.max_components = max_components

    def __call__(self):
        self.set_surface_mask()

    def set_surface_mask(self):
        surface = self.surface
        v, t = surface.vertices, surface.triangles
        if t is None:
>           surface.auto_remask_triangles = self
            return

[Tristan Croll]

Sorry - meant to add: I'm happy that everything's working on my end. Reassigned to you for the minor issue with ZoneMask. Feel free to close.

[Tom Goddard]

Ok. I made setting a surface zone before a surface has been computed correctly update.