Add option to compute topology based on corners only

splipy/SplineModel.py

@@ -114,7 +114,7 @@ def __init__(self, perm, flip):
         self.perm_inv = tuple(perm.index(d) for d in range(len(perm)))
 
     @classmethod
-    def compute(cls, cpa, cpb=None):
+    def compute(cls, cpa, cpb=None, interior=True):
         """Compute and return a new orientation object representing the mapping between
         `cpa` (the reference system) and `cpb` (the mapped system).
 
@@ -136,25 +136,38 @@ def compute(cls, cpa, cpb=None):
 
         shape_b = cpb.controlpoints.shape
 
-        # Deal with the easy cases: dimension mismatch, and
-        # comparing the shapes as multisets
+        # Easy error checking
         if len(shape_a) != len(shape_b):
             raise OrientationError("Mismatching parametric dimensions")
+
+        cpa = cpa.controlpoints
+        cpb = cpb.controlpoints
+        if not interior:
+            # Pick out just the corners
+            indexes = np.ix_(*[[0,-1] for _ in range(pardim)])
+            indexes = list(indexes) + [slice(None)]
+            cpa = cpa[indexes]
+            cpb = cpb[indexes]
+            shape_a = cpa.shape
+            shape_b = cpb.shape
+
+        # Deal with the rest of the easy cases: dimension mismatch, and
+        # comparing the shapes as multisets
         if shape_a[-1] != shape_b[-1]:
             raise OrientationError("Mismatching physical dimensions")
         if Counter(shape_a) != Counter(shape_b):
             raise OrientationError("Non-matching objects")
 
         # Enumerate all permutations of directions
         for perm in permutations(range(pardim)):
-            transposed = cpb.controlpoints.transpose(perm + (pardim,))
+            transposed = cpb.transpose(perm + (pardim,))
             if transposed.shape != shape_a:
                 continue
             # Enumerate all possible direction reversals
             for flip in product([False, True], repeat=pardim):
                 slices = tuple(slice(None, None, -1) if f else slice(None) for f in flip)
                 test_b = transposed[slices + (slice(None),)]
-                if np.allclose(cpa.controlpoints, test_b,
+                if np.allclose(cpa, test_b,
                                rtol=state.controlpoint_relative_tolerance,
                                atol=state.controlpoint_absolute_tolerance):
                     if all([cpa.bases[i].matches(cpb.bases[perm[i]], reverse=flip[i]) for i in range(pardim)]):
@@ -254,13 +267,16 @@ class TopologicalNode(object):
         of any kind.
     """
 
-    def __init__(self, obj, lower_nodes):
+    def __init__(self, catalogue, obj, lower_nodes):
         """Initialize a `TopologicalNode` object associated with the given
         `SplineObject` and lower order nodes.
 
+        :param ObjectCatalogue catalogue: The catalogue to which this node
+            belongs
         :param SplineObject obj: The underlying spline object
         :param lower_nodes: A nested list of lower order nodes
         """
+        self.catalogue = catalogue
         self.obj = obj
         self.lower_nodes = lower_nodes
         self.higher_nodes = {}
@@ -288,9 +304,9 @@ def view(self, other_obj=None):
             underlying object
         """
         if other_obj:
-            orientation = Orientation.compute(self.obj, other_obj)
+            orientation = self.catalogue.make_orientation(self.obj, other_obj)
         else:
-            orientation = Orientation.compute(self.obj)
+            orientation = self.catalogue.make_orientation(self.obj)
         return NodeView(self, orientation)
 
 
@@ -334,7 +350,9 @@ def section(self, *args, **kwargs):
 
         # The underlying lower-order node may not have an orientation that
         # matches the higher-order node, so we need to compose two orientations
-        ref_ori = Orientation.compute(node.obj, self.node.obj.section(*section, unwrap_points=False))
+        ref_ori = self.node.catalogue.make_orientation(
+            node.obj, self.node.obj.section(*section, unwrap_points=False)
+        )
         my_ori = self.orientation.view_section(section)
 
         return NodeView(node, ref_ori * my_ori)
@@ -372,7 +390,7 @@ class ObjectCatalogue(object):
     at most `pardim` parametric directions.
     """
 
-    def __init__(self, pardim):
+    def __init__(self, pardim, interior=True):
         """Initialize a catalogue for objects of parametric dimension
         `pardim`.
         """
@@ -381,10 +399,15 @@ def __init__(self, pardim):
         # Internal mapping from tuples of lower-order nodes to lists of nodes
         self.internal = {}
 
+        # Function for computing orientations
+        self.make_orientation = (
+            lambda *args, **kwargs: Orientation.compute(*args, interior=interior, **kwargs)
+        )
+
         # Each catalogue has a catalogue of lower dimension
         # For points, we use a VertexDict
         if pardim > 0:
-            self.lower = ObjectCatalogue(pardim - 1)
+            self.lower = ObjectCatalogue(pardim - 1, interior=interior)
         else:
             self.lower = VertexDict()
 
@@ -408,7 +431,7 @@ def lookup(self, obj, add=False):
         # Special case for points: self.lower is a mapping from array to node
         if self.pardim == 0:
             if add:
-                node = TopologicalNode(obj, [])
+                node = TopologicalNode(self, obj, [])
                 return self.lower.setdefault(obj.controlpoints, node).view()
             return self.lower[obj.controlpoints].view()
 
@@ -436,7 +459,7 @@ def lookup(self, obj, add=False):
         except (KeyError, OrientationError):
             if not add:
                 raise KeyError("No such object found")
-            node = TopologicalNode(obj, lower_nodes)
+            node = TopologicalNode(self, obj, lower_nodes)
             # Assign the new node to each possible permutation of lower-order
             # nodes. This is slight overkill since some of these permutations
             # are invalid, but c'est la vie.
@@ -480,11 +503,11 @@ def nodes(self, pardim):
 
 class SplineModel(object):
 
-    def __init__(self, pardim=3, dimension=3, objs=[]):
+    def __init__(self, pardim=3, dimension=3, objs=[], interior=True):
         self.pardim = pardim
         self.dimension = dimension
 
-        self.catalogue = ObjectCatalogue(pardim)
+        self.catalogue = ObjectCatalogue(pardim, interior=interior)
         self.add_patches(objs)
 
     def add_patch(self, obj):