@@ -63,6 +63,29 @@ class Triangulation(object):
6363
6464 hull -- list of point_id's giving the nodes which form the convex hull
6565 of the point set. This list is sorted in counter-clockwise order.
66+
67+ Duplicate points.
68+ If there are no duplicate points, Triangulation stores the specified
69+ x and y arrays and there is no difference between the client's and
70+ Triangulation's understanding of point indices used in edge_db,
71+ triangle_nodes and hull.
72+
73+ If there are duplicate points, they are removed from the stored
74+ self.x and self.y as the underlying delaunay code cannot deal with
75+ duplicates. len(self.x) is therefore equal to len(x) minus the
76+ number of duplicate points. Triangulation's edge_db, triangle_nodes
77+ and hull refer to point indices in self.x and self.y, for internal
78+ consistency within Triangulation and the corresponding Interpolator
79+ classes. Client code must take care to deal with this in one of
80+ two ways:
81+
82+ 1. Ignore the x,y it specified in Triangulation's constructor and
83+ use triangulation.x and triangulation.y instead, as these are
84+ consistent with edge_db, triangle_nodes and hull.
85+
86+ 2. If using the x,y the client specified then edge_db,
87+ triangle_nodes and hull should be passed through the function
88+ to_client_point_indices() first.
6689 """
6790 def __init__ (self , x , y ):
6891 self .x = np .asarray (x , dtype = np .float64 )
@@ -72,38 +95,46 @@ def __init__(self, x, y):
7295 raise ValueError ("x,y must be equal-length 1-D arrays" )
7396
7497 self .old_shape = self .x .shape
75- j_unique = self ._collapse_duplicate_points ()
98+ duplicates = self ._get_duplicate_point_indices ()
7699
77- if j_unique . shape != self . x . shape :
100+ if len ( duplicates ) > 0 :
78101 warnings .warn (
79102 "Input data contains duplicate x,y points; some values are ignored." ,
80103 DuplicatePointWarning ,
81104 )
82- self .j_unique = j_unique
105+
106+ # self.j_unique is the array of non-duplicate indices, in
107+ # increasing order.
108+ self .j_unique = np .delete (np .arange (len (self .x )), duplicates )
83109 self .x = self .x [self .j_unique ]
84110 self .y = self .y [self .j_unique ]
85111 else :
86112 self .j_unique = None
87113
114+ # If there are duplicate points, need a map of point indices used
115+ # by delaunay to those used by client. If there are no duplicate
116+ # points then the map is not needed. Either way, the map is
117+ # conveniently the same as j_unique, so share it.
118+ self ._client_point_index_map = self .j_unique
88119
89120 self .circumcenters , self .edge_db , self .triangle_nodes , \
90121 self .triangle_neighbors = delaunay (self .x , self .y )
91122
92123 self .hull = self ._compute_convex_hull ()
93124
94- def _collapse_duplicate_points (self ):
95- """Generate index array that picks out unique x,y points.
96-
97- This appears to be required by the underlying delaunay triangulation
98- code.
125+ def _get_duplicate_point_indices (self ):
126+ """Return array of indices of x,y points that are duplicates of
127+ previous points. Indices are in no particular order.
99128 """
100- # Find the indices of the unique entries
129+ # Indices of sorted x,y points.
101130 j_sorted = np .lexsort (keys = (self .x , self .y ))
102- mask_unique = np .hstack ([
103- True ,
104- (np .diff (self .x [j_sorted ]) != 0 ) | (np .diff (self .y [j_sorted ]) ! = 0 ),
131+ mask_duplicates = np .hstack ([
132+ False ,
133+ (np .diff (self .x [j_sorted ]) == 0 ) & (np .diff (self .y [j_sorted ]) = = 0 ),
105134 ])
106- return j_sorted [mask_unique ]
135+
136+ # Array of duplicate point indices, in no particular order.
137+ return j_sorted [mask_duplicates ]
107138
108139 def _compute_convex_hull (self ):
109140 """Extract the convex hull from the triangulation information.
@@ -131,6 +162,16 @@ def _compute_convex_hull(self):
131162
132163 return hull
133164
165+ def to_client_point_indices (self , array ):
166+ """Converts any array of point indices used within this class to
167+ refer to point indices within the (x,y) arrays specified in the
168+ constructor before duplicates were removed.
169+ """
170+ if self ._client_point_index_map is not None :
171+ return self ._client_point_index_map [array ]
172+ else :
173+ return array
174+
134175 def linear_interpolator (self , z , default_value = np .nan ):
135176 """Get an object which can interpolate within the convex hull by
136177 assigning a plane to each triangle.
0 commit comments