- hpc.geometry.average_numba module.

Author: andrewtliew

src/compas/geometry/average.py

@@ -1,6 +1,6 @@
-from __future__ import print_function
 from __future__ import absolute_import
 from __future__ import division
+from __future__ import print_function
 
 from compas.utilities import window
 
@@ -14,7 +14,7 @@
 
 
 __author__    = ['Tom Van Mele', ]
-__copyright__ = 'Copyright 2016 - Block Research Group, ETH Zurich'
+__copyright__ = 'Copyright 2017 - Block Research Group, ETH Zurich'
 __license__   = 'MIT License'
 __email__     = '[email protected]'
 

src/compas/hpc/__init__.py

@@ -33,6 +33,9 @@
 geometry
 ========
 
+basic
+-----
+
 .. autosummary::
     :toctree: generated/
     :nosignatures:
@@ -78,6 +81,20 @@
     circle_from_points_numba
     circle_from_points_xy_numba
 
+average
+-------
+
+.. autosummary::
+    :toctree: generated/
+    :nosignatures:
+
+    centroid_points_numba
+    centroid_points_xy_numba
+    midpoint_point_point_numba
+    midpoint_point_point_xy_numba
+    center_of_mass_polyline_numba
+    center_of_mass_polyline_xy_numba
+
 
 core
 ====

src/compas/hpc/geometry/average_numba.py

@@ -0,0 +1,264 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+
+from compas.hpc.geometry import add_vectors_numba
+from compas.hpc.geometry import add_vectors_xy_numba
+from compas.hpc.geometry import cross_vectors_numba
+from compas.hpc.geometry import dot_vectors_numba
+from compas.hpc.geometry import length_vector_numba
+from compas.hpc.geometry import length_vector_xy_numba
+from compas.hpc.geometry import scale_vector_numba
+from compas.hpc.geometry import scale_vector_xy_numba
+from compas.hpc.geometry import subtract_vectors_numba
+from compas.hpc.geometry import subtract_vectors_xy_numba
+from compas.hpc.geometry import sum_vectors_numba
+
+from numba import f8
+from numba import i8
+from numba import jit
+from numba import prange
+
+from numpy import array
+
+
+__author__    = ['Andrew Liew <[email protected]>']
+__copyright__ = 'Copyright 2017, BLOCK Research Group - ETH Zurich'
+__license__   = 'MIT License'
+__email__     = '[email protected]'
+
+
+__all__ = [
+    'centroid_points_numba',
+    'centroid_points_xy_numba',
+    'midpoint_point_point_numba',
+    'midpoint_point_point_xy_numba',
+    'center_of_mass_polyline_numba',
+    'center_of_mass_polyline_xy_numba',
+    # 'center_of_mass_polyhedron_numba',
+]
+
+
+@jit(f8[:](f8[:, :]), nogil=True, nopython=True, parallel=True)
+def centroid_points_numba(u):
+    """Compute the centroid of a set of points.
+
+    Parameters
+    ----------
+    u : array
+        An array of XYZ coordinates.
+
+    Returns
+    -------
+    array
+        Centroid of the points.
+    """
+    m = u.shape[0]
+    f = 1 / m
+    return scale_vector_numba(sum_vectors_numba(u, axis=0), factor=f)
+
+
+@jit(f8[:](f8[:, :]), nogil=True, nopython=True, parallel=True)
+def centroid_points_xy_numba(u):
+    """Compute the centroid of a set of points lying in the XY plane.
+
+    Parameters
+    ----------
+    u : array
+        An array of XY(Z) coordinates in the XY plane.
+
+    Returns
+    -------
+    array
+        Centroid of the points (Z = 0.0).
+    """
+    u[:, 2] = 0
+    return centroid_points_numba(u)
+
+
+@jit(f8[:](f8[:], f8[:]), nogil=True, nopython=True, parallel=True)
+def midpoint_point_point_numba(u, v):
+    """Compute the midpoint of two points.
+
+    Parameters
+    ----------
+    u : array
+        XYZ coordinates of the first point.
+    v : array
+        XYZ coordinates of the second point.
+
+    Returns
+    -------
+    array
+        XYZ coordinates of the midpoint.
+    """
+    return scale_vector_numba(add_vectors_numba(u, v), factor=0.5)
+
+
+@jit(f8[:](f8[:], f8[:]), nogil=True, nopython=True, parallel=True)
+def midpoint_point_point_xy_numba(u, v):
+    """Compute the midpoint of two points lying in the XY-plane.
+
+    Parameters
+    ----------
+    u : array
+        XY(Z) coordinates of the first point in the XY plane.
+    v : array
+        XY(Z) coordinates of the second point in the XY plane.
+
+    Returns
+    -------
+    array
+        XYZ (Z = 0.0) coordinates of the midpoint.
+    """
+    return scale_vector_xy_numba(add_vectors_xy_numba(u, v), factor=0.5)
+
+
+@jit(f8[:](f8[:, :]), nogil=True, nopython=True, parallel=False)
+def center_of_mass_polyline_numba(polyline):
+    """Compute the center of mass of polyline edges defined as an array of points.
+
+    Parameters
+    ----------
+    polyline : array
+        XYZ coordinates representing the corners of a polyline (m x 3).
+
+    Returns
+    -------
+    array
+        The XYZ coordinates of the center of mass.
+    """
+    L = 0
+    cx = 0
+    cy = 0
+    cz = 0
+    m = polyline.shape[0]
+    ii = list(range(1, m)) + [0]
+    for i in prange(m):
+        p1 = polyline[i, :]
+        p2 = polyline[ii[i], :]
+        d = length_vector_numba(subtract_vectors_numba(p2, p1))
+        cx += 0.5 * d * (p1[0] + p2[0])
+        cy += 0.5 * d * (p1[1] + p2[1])
+        cz += 0.5 * d * (p1[2] + p2[2])
+        L += d
+    c = array([cx, cy, cz])
+    f = 1 / L
+    return scale_vector_numba(c, factor=f)
+
+
+@jit(f8[:](f8[:, :]), nogil=True, nopython=True, parallel=False)
+def center_of_mass_polyline_xy_numba(polyline):
+    """Compute the center of mass of polyline edges in the XY plane, defined as an array of points.
+
+    Parameters
+    ----------
+    polyline : array
+        XY(Z) coordinates representing the corners of a polyline (m x 3).
+
+    Returns
+    -------
+    array
+        The XY(Z) coordinates of the center of mass.
+    """
+    L = 0
+    cx = 0
+    cy = 0
+    m = polyline.shape[0]
+    ii = list(range(1, m)) + [0]
+    for i in prange(m):
+        p1 = polyline[i, :]
+        p2 = polyline[ii[i], :]
+        d = length_vector_xy_numba(subtract_vectors_xy_numba(p2, p1))
+        cx += 0.5 * d * (p1[0] + p2[0])
+        cy += 0.5 * d * (p1[1] + p2[1])
+        L += d
+    c = array([cx, cy, 0.])
+    f = 1 / L
+    return scale_vector_numba(c, factor=f)
+
+
+@jit(f8[:](f8[:, :], i8[:, :]), nogil=True, nopython=True, parallel=False)  # Needs checking
+def center_of_mass_polyhedron_numba(vertices, faces):
+    """Compute the center of mass of the edges of a polyhedron.
+
+    Parameters
+    ----------
+    vertices : array
+        Array containing the XYZ coordinates of the polyhedron vertices (n x 3).
+    faces : array
+        Array contianing the indices of triangle faces of the polyhedron (m x 3).
+
+    Return
+    ------
+    array
+        The XYZ coordinates of the polyhedron edges' centre of mass.
+    """
+    m = faces.shape[0]
+    V = 0.
+    x  = 0.
+    y  = 0.
+    z  = 0.
+    ex = array([1., 0., 0.])
+    ey = array([0., 1., 0.])
+    ez = array([0., 0., 1.])
+    ii = [1, 2, 0]
+
+    for i in prange(m):
+        a = vertices[faces[i, 0]]
+        b = vertices[faces[i, 1]]
+        c = vertices[faces[i, 2]]
+        ab = subtract_vectors_numba(b, a)
+        ac = subtract_vectors_numba(c, a)
+        n  = cross_vectors_numba(ab, ac)
+        V += dot_vectors_numba(a, n)
+        nx = dot_vectors_numba(n, ex)
+        ny = dot_vectors_numba(n, ey)
+        nz = dot_vectors_numba(n, ez)
+
+        for k in prange(3):
+            ab = add_vectors_numba(vertices[faces[i, k]], vertices[faces[i, ii[k]]])
+            x += nx * dot_vectors_numba(ab, ex)**2
+            y += ny * dot_vectors_numba(ab, ey)**2
+            z += nz * dot_vectors_numba(ab, ez)**2
+
+    if V < 10**(-9):
+        V = 0.
+        d = 1. / 48.
+    else:
+        V /= 6.
+        d = 1. / 48. / V
+    x *= d
+    y *= d
+    z *= d
+
+    return array([x, y, z])
+
+
+# ==============================================================================
+# Main
+# ==============================================================================
+
+if __name__ == "__main__":
+
+    from time import time
+
+    u = array([1., 2., 3.])
+    v = array([4., 5., 6.])
+    c = array([[0., 0., 1.], [3., 4., 1.], [6., 0., 1.]])
+    d = array([[0, 1, 2], [1, 2, 0]])
+    e = array([[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]])
+
+    tic = time()
+
+    for i in range(10**5):
+
+        # centroid_points_numba(c)
+        # centroid_points_xy_numba(c)
+        # midpoint_point_point_numba(u, v)
+        # midpoint_point_point_xy_numba(u, v)
+        # center_of_mass_polyline_numba(c)
+        # center_of_mass_polyline_xy_numba(c)
+        center_of_mass_polyhedron_numba(c, array([[0, 1, 2]]))
+
+    print(time() - tic)