From 056b8badc3a9abc52df560c0be89f6d6e7fdb87d Mon Sep 17 00:00:00 2001 From: Eivind Fonn Date: Wed, 19 Apr 2017 10:46:15 +0200 Subject: [PATCH] Add option to compute topology based on corners only --- splipy/SplineModel.py | 53 +++++++++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/splipy/SplineModel.py b/splipy/SplineModel.py index 6a723885..6b473576 100644 --- a/splipy/SplineModel.py +++ b/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,10 +136,23 @@ 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): @@ -147,14 +160,14 @@ def compute(cls, cpa, cpb=None): # 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):