Support for more spaces/geodesics means a new version.

meander/__init__.py

@@ -1,2 +1,4 @@
 from .contour_compute import compute_contours
+from .contour_compute import planar_contours
+from .contour_compute import spherical_contours
 del contour_compute

meander/contour_compute.py

@@ -120,7 +120,7 @@ def _simplex_intersections(samples, levels, simplices):
 
     return simplex_intersections_by_level
 
-def _cartesian_geodesic(p0, p1, x):
+def _planar_geodesic(p0, p1, x):
     # This function describes the geodesic between p0 and p1 in a 2d cartesian space
     # p0 and p1 are assumed to be of the form (x0, y0) and (x1, y1)
     assert(x >= 0.0 and x <= 1.0)
@@ -132,7 +132,7 @@ def _spherical_geodesic(p0, p1, x):
     # p0 and p1 are assumed to be of the form (theta0, phi0) and (theta1, phi1)
     # where theta and phi are the polar and azimuthal angles respectively measured in radians
     assert(x >= 0.0 and x <= 1.0)
-    v0, v0 = [np.array((np.cos(phi)*np.sin(theta), np.sin(phi)*np.sin(theta), np.cos(theta))) for theta, phi in [p0, p1]]
+    v0, v1 = [np.array((np.cos(phi)*np.sin(theta), np.sin(phi)*np.sin(theta), np.cos(theta))) for theta, phi in [p0, p1]]
     c = np.sqrt(np.sum((v1-v0)**2))
     theta = 2.0*np.arcsin(c/2.0)
     theta0 = x*theta
@@ -149,7 +149,7 @@ def _spherical_geodesic(p0, p1, x):
     p = (np.arccos(v3[2]), np.arctan2(v3[1], v3[0]))
     return p
 
-def _interpolate_simplex_intersections(sample_points, samples, levels, simplex_intersections_by_level, geodesic=_cartesian_geodesic):
+def _interpolate_simplex_intersections(sample_points, samples, levels, simplex_intersections_by_level, geodesic=_planar_geodesic):
     contours_by_level = []
     for level_i, level in enumerate(levels):
         connected_line_segments = simplex_intersections_by_level[level_i]
@@ -158,76 +158,133 @@ def _interpolate_simplex_intersections(sample_points, samples, levels, simplex_i
         line_ps = dict()
 
         for line in lines_in_contours:
+            line = tuple(line)
             p0, p1 = line
             (p0, y0), (p1, y1) = sorted([(p0, samples[p0]), (p1, samples[p1])], key=lambda x: x[1])
             x = (level-y0) / (y1-y0)
             p = geodesic(sample_points[p0], sample_points[p1], x)
-            line_ps[li] = p
-        line_ps = np.array(line_ps)
+            line_ps[line] = p
 
-        contours_by_level.append([np.array([line_ps[li] for li in index_contour]) for index_contour in connect_line_segments(simplex_lines)])
+        contours_by_level.append([np.array([line_ps[tuple(li)] for li in index_contour]) for index_contour in connected_line_segments])
     return contours_by_level
 
-def compute_contours(sample_points, samples, level):
-    """ Compute contours (iso-lines of a scalar field) in a 2d cartesian space """
-    # compute the delaunay triangulation so we have a set of
-    # triangles to run the meandering triangles algorithm on
+def _compute_planar_simplices(sample_points):
     tri = scipy.spatial.Delaunay(sample_points)
     simplices = tri.simplices
+    return simplices
+
+def _compute_spherical_simplices(sample_points):
+    sample_points = [(np.cos(phi)*np.sin(theta), np.sin(phi)*np.sin(theta), np.cos(theta)) for theta, phi in sample_points]
+    hull = scipy.spatial.ConvexHull(sample_points)
+    simplices = hull.simplices
+    return simplices
+
+def _compute_simplices(geodesic, simplex_function, sample_points):
+    try:
+        simplex_array = np.asarray(simplex_function)
+        simplices_shape = simplex_array.shape
+        assert(len(simplices_shape) == 2)
+        n_samples = len(sample_points)
+        if simplices_shape[0] == n_samples:
+            assert(simplices_shape[1] == 3)
+            return simplex_array
+        elif simplices_shape[1] == n_samples:
+            assert(simplices_shape[0] == 3)
+            return simplex_array.T
+    except:
+        pass
+
+    if callable(simplex_function):
+        try:
+            simplices = simplex_function(sample_points)
+            return simplices
+        except:
+            pass
+
+    if simplex_function is None:
+        if type(geodesic) is str:
+            simplex_function = geodesic
+
+    if simplex_function == 'planar':
+        return _compute_planar_simplices(sample_points)
+    elif simplex_function == 'spherical':
+        return _compute_spherical_simplices(sample_points)
+    else:
+        raise ValueError("Not sure what to do with simplex_function = %s" % str(simplex_function))
+
+def _get_geodesic_function(geodesic):
+    if callable(geodesic):
+        return geodesic
+    elif geodesic == 'planar':
+        return _planar_geodesic
+    elif geodesic == 'spherical':
+        return _spherical_geodesic
+    else:
+        raise ValueError("Not sure what to do with geodesic = %s" % str(geodesic))
+
+def compute_contours(sample_points, samples, levels, geodesic='planar', simplices=None):
+    """ Compute contours (iso-lines of a scalar field) in a 2d space
+
+        Parameters
+        ----------
+        sample_points : array_like
+            The points in the 2d space for which there are samples
+            of the scalar field. This should be structured as
+            [(x0,y0), (x1,y1), ...]
+        samples : array_like
+            The values of the scalar field at the locations of sample_points
+        levels : array_like
+            The values of the scalar field at which the contours are computed
+        geodesic : str or callable, optional
+            The geodesic of the space in which the sample_points exist.
+            If geodesic is a str, it can have values ['planar', 'spherical']
+            to specify if the points exist in a planar space or on the surface
+            of a sphere.
+            If geodesic is callable it should have the form geodesic(p0, p1, x)
+            where p0 and p1 are points in the space and x is a number between 0
+            and 1 that describes a point between p0 and p1 along the geodesic. In
+            this case geodesic should return the point corresponding to x.
+            The default is 'planar'.
+        simplices : None or array_like or callable, optional
+            If None, the simplices are computed via the Delaunay triangulation
+            on either a planar or spherical surface depending on the geodesic
+            If array_like, the simplices are assumed to be of the form
+            [(pi0, pj0, pk0), (pi1, pj1, pk1), ...] where the pi,pj,pk refer to indices
+            in sample_points that specify the points that compose each simplex
+            If callable, simplices is assumed to take the sample_points and return the
+            an array_like version of simplices
+            The default is None
+        Returns
+        -------
+        contours_by_level : list(list(list(point)))
+            The contours for each level
+            Outermost list indexes by level
+            Second list indexes by contours at a particular level
+            Third list indexes by points in each contour
+            Points are of the same form as sample_points
+    """
+
+    simplices = _compute_simplices(geodesic, simplices, sample_points)
+
+    simplex_intersections_by_level = _simplex_intersections(samples, levels, simplices)
+
+    geodesic_function = _get_geodesic_function(geodesic)
+
+    contours_by_level = _interpolate_simplex_intersections(sample_points, samples, levels, simplex_intersections_by_level, geodesic=geodesic_function)
 
-    # to simplify the problem each point has an index
-    # and each line has an index
+    return contours_by_level
 
-    simplex_lines = [] # lines between the simplex edges that cross the level
-    lines = [] # the lines that compose the simplex edges
-    li_to_si = dict() # maps line index -> simplex indices
-    l_to_li = dict() # maps line -> line index
-    for si, points in enumerate(simplices):
-        points = sorted(points)
-        # points and lines are ordered so lx does not contain px (px is opposite lx)
-        p0, p1, p2 = points
-        these_lines = [(p1,p2), (p0,p2), (p0,p1)]
+def planar_contours(sample_points, samples, levels, simplices=None):
+    """ Compute contours (iso-lines of a scalar field) in a 2d planar space
+        See compute_contours for more information.
+        compute_contours(sample_points, samples, levels, geodesic='planar', simplices=simplices)
+    """
+    return compute_contours(sample_points, samples, levels, geodesic='planar', simplices=simplices)
 
-        # fill the maps
-        for l in these_lines:
-            if l not in l_to_li:
-                li_to_si[len(lines)] = []
-                l_to_li[l] = len(lines)
-                lines.append(l)
-            li = l_to_li[l]
-            li_to_si[li].append(si)
+def spherical_contours(sample_points, samples, levels, simplices=None):
+    """ Compute contours (iso-lines of a scalar field) in a 2d space on the surface of a unit sphere
+        See compute_contours for more information.
+        compute_contours(sample_points, samples, levels, geodesic='spherical', simplices=simplices)
+    """
+    return compute_contours(sample_points, samples, levels, geodesic='spherical', simplices=simplices)
 
-        # which points are above / below the level
-        b = np.array([samples[p] >= level for p in points]).astype(bool)
-        b_sum = np.sum(b)
-
-        # skip triangles thar are all above or all below the level
-        # for triangles spanning the level, get the lines that cross the level
-        if b_sum == 1:
-            line_indices = np.array([l_to_li[l] for l in these_lines])[np.array([i for i in xrange(3) if i != np.argmax(b)])]
-        elif b_sum == 2:
-            line_indices = np.array([l_to_li[l] for l in these_lines])[np.array([i for i in xrange(3) if i != np.argmin(b)])]
-        else:
-            continue
-        simplex_lines.append(tuple(sorted(line_indices)))
-
-    # simplex_lines specifies which simplex lines to draw the contour lines between
-    # but we still need to know where on those lines to put the points
-    # we can estimate where the level is by linearly interpolating between the sample points
-
-    line_ps = np.zeros(len(lines)).tolist()
-    #line_xs = np.zeros(len(lines)).tolist()
-
-    for li in xrange(len(lines)):
-        p0, p1 = lines[li]
-        (p0, y0), (p1, y1) = sorted([(p0, samples[p0]), (p1, samples[p1])], key=lambda x: x[1])
-        x = (level-y0) / (y1-y0)
-        p = np.array(sample_points[p0]) * (1.0-x) + np.array(sample_points[p1]) * x
-        #line_xs[li] = x
-        line_ps[li] = p
-    line_ps = np.array(line_ps)
-    #line_xs = np.array(line_xs)
-    #mask = np.logical_and(line_xs < 1.0, line_xs > 0.0)
-    #simplex_lines = [(l0, l1) for l0, l1 in simplex_lines if mask[l0] and mask[l1]]
-
-    return [np.array([line_ps[li] for li in index_contour]) for index_contour in connect_line_segments(simplex_lines)]

setup.py

@@ -1,9 +1,14 @@
 from setuptools import setup
 
+def readme():
+    with open('README.md') as f:
+        return f.read()
+
 setup(
         name='meander',
-        version='0.0.1',
+        version='0.0.2',
         description='Tools for computing contours on 2d surfaces.',
+        long_description=readme(),
         url='https://github.com/austinschneider/meander',
         author='Austin Schneider',
         author_email='[email protected]',