MAINT: Unify calculation of normal vectors from polygons

[eric-wieser]

This combines get_normals and _generate_normals, and eliminates all other calls to np.cross. This came up in #9990, but I wanted to postpone it to (this) later PR.

get_normals and _generate_normals were profiled, and it was found that vectorizing np.cross like in get_normals was faster:

import numpy as np

def get_normals(polygons):
    v1 = np.empty((len(polygons), 3))
    v2 = np.empty((len(polygons), 3))
    for poly_i, ps in enumerate(polygons):
        # pick three points around the polygon at which to find the
        # normal doesn't vectorize because polygons is jagged
        i1, i2, i3 = 0, len(ps)//3, 2*len(ps)//3
        v1[poly_i, :] = ps[i1, :] - ps[i2, :]
        v2[poly_i, :] = ps[i2, :] - ps[i3, :]
    return np.cross(v1, v2)

def _generate_normals(polygons):
    normals = []
    for verts in polygons:
        v1 = np.array(verts[0]) - np.array(verts[1])
        v2 = np.array(verts[2]) - np.array(verts[0])
        normals.append(np.cross(v1, v2))
    return np.array(normals)

def _generate_normals_new(polygons):
    if isinstance(polygons, np.ndarray):
        # optimization: polygons all have the same number of points, so can
        # vectorize
        n = polygons.shape[-2]
        i1, i2, i3 = 0, n//3, 2*n//3
        v1 = polygons[..., i1, :] - polygons[..., i2, :]
        v2 = polygons[..., i2, :] - polygons[..., i3, :]
    else:
        # The subtraction doesn't vectorize because polygons is jagged.
        v1 = np.empty((len(polygons), 3))
        v2 = np.empty((len(polygons), 3))
        for poly_i, ps in enumerate(polygons):
            n = len(ps)
            i1, i2, i3 = 0, n//3, 2*n//3
            v1[poly_i, :] = ps[i1, :] - ps[i2, :]
            v2[poly_i, :] = ps[i2, :] - ps[i3, :]
    return np.cross(v1, v2)

polygons_jagged = [
    np.random.rand(np.random.randint(10, 1000), 3)
    for i in range(100)
]

polygons_uniform = np.random.rand(100, 100, 3)

In [10]: %timeit _generate_normals(polygons_jagged)
8.43 ms ± 148 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [11]: %timeit get_normals(polygons_jagged)
583 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [12]: %timeit _generate_normals_new(polygons_jagged)
593 µs ± 26 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


In [14]: %timeit _generate_normals(polygons_uniform)
8.45 ms ± 156 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [15]: %timeit get_normals(polygons_uniform)
599 µs ± 6.13 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [13]: %timeit _generate_normals_new(polygons_uniform)
110 µs ± 10.1 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

PR Checklist

• Code is Flake 8 compliant

Do not apply, since this is just a refactor:

• Has Pytest style unit tests
• New features are documented, with examples if plot related
• Documentation is sphinx and numpydoc compliant
• Added an entry to doc/users/next_whats_new/ if major new feature (follow instructions in README.rst there)
• Documented in doc/api/api_changes.rst if API changed in a backward-incompatible way

[eric-wieser]

Branch name hints that my goal here is to add shading to ax3d.voxels, in a later PR

[eric-wieser]

An alternate implementation with less duplication, but with less clarity, could look like:

        def _inplane_vectors(polygons, v1_out=None, v2_out=None):
            n = polygons.shape[-2]
            i1, i2, i3 = 0, n//3, 2*n//3
            v1 = np.subtract(ps[..., i1, :], ps[..., i2, :], out=v1_out)
            v2 = np.subtract(ps[..., i2, :], ps[..., i3, :], out=v2_out)
            return v1, v2

        if isinstance(polygons, np.ndarray):
            # optimization: polygons all have the same number of points, so can
            # vectorize
            v1, v2 = _inplane_vectors(polygons)
        else:
            # The subtraction doesn't vectorize because polygons is jagged.
            v1 = np.empty((len(polygons), 3))
            v2 = np.empty((len(polygons), 3))
            for poly_i, ps in enumerate(polygons):
                _inplane_vectors(ps, v1_out=v1[poly_i, :], v2_out=v2[poly_i, :])
        return np.cross(v1, v2)

Let me know which is preferable

[eric-wieser]

This fails because it changes the sense of the polygons. It seems that some functions use a right-hand-rule for the "front" face of a polygon, while others use a left-hand rule. Related: #12138

[eric-wieser]

Just to give an update - this is now waiting on #12259

[eric-wieser]

Rebased and ready to go

[timhoffm]

Looks good otherwise. Please rebase.

[timhoffm]

Why the i index ?

[eric-wieser]

Because the value of M can change between each list item - polygons[0] has shape (M_0, 3), polygons[1] has shape (M_1, 3), and it is likely that M_1 != M_0

[eric-wieser]

@timhoffm: merged using the web UI, since it was doc-only. Feel free to squash to elide the merge commit

[timhoffm]

Is there the sign change here intended? This is the convention of get_normals. _generate_normals had it the other way round. I know that there have been some orientation issues. But I don‘t know the state. Just want to make sure the sign change is notbslipping in unintendedly.

[eric-wieser]

Sign change is only in v1 and v2, it cancels out in the cross product, so it doesn't affect the return value.

I picked the convention from get_normals because that was easiest. If you want me to flip both subtractions here, I can, but it won't make any difference.