From 6893bc71a64cff988087d45631fbf6d40b8e7b2e Mon Sep 17 00:00:00 2001 From: Benjamin Date: Sat, 21 Sep 2019 10:50:23 +0200 Subject: [PATCH 01/10] Adds bounding hull geometric and image manipulation routines and test --- africanus/install/requirements.py | 2 +- africanus/linalg/geometry.py | 574 +++++++++++++++++++++++ africanus/testing/tests/test_geometry.py | 219 +++++++++ 3 files changed, 794 insertions(+), 1 deletion(-) create mode 100644 africanus/linalg/geometry.py create mode 100644 africanus/testing/tests/test_geometry.py diff --git a/africanus/install/requirements.py b/africanus/install/requirements.py index f29ff5ae6..2dfd20084 100644 --- a/africanus/install/requirements.py +++ b/africanus/install/requirements.py @@ -30,7 +30,7 @@ 'astropy': ['astropy >= 2.0.0, < 3.0; python_version <= "2.7"', 'astropy >= 3.0; python_version >= "3.0"'], 'python-casacore': ['python-casacore == 3.0.0'], - 'testing': ['pytest', 'flaky', 'pytest-flake8'] + 'testing': ['pytest', 'flaky', 'pytest-flake8', 'pathlib'] } _non_cuda_extras = [er for n, er in extras_require.items() if n != "cuda"] diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py new file mode 100644 index 000000000..b4da4d92c --- /dev/null +++ b/africanus/linalg/geometry.py @@ -0,0 +1,574 @@ +import numpy as np +import scipy.stats as sstats +import scipy.signal as ssig +import scipy.spatial as spat +import copy + +import time +def timeit(method): + def timed(*args, **kw): + ts = time.time() + result = method(*args, **kw) + te = time.time() + if 'log_time' in kw: + name = kw.get('log_name', method.__name__.upper()) + kw['log_time'][name] = int((te - ts) * 1000) + else: + print(('%r %2.2f ms' % \ + (method.__name__, (te - ts) * 1000))) + + return result + return timed + +class BoundingConvexHull(object): + def __init__(self, list_hulls, name="unnamed", mask = None, check_mask_outofbounds=True): + """ Initializes a bounding convex hull around a list of bounding convex hulls or series of + points. A unity-weighted mask is computed for the region that falls within this convex hull + if a mask of (y, x) coordinates is not provided. Otherwise if a mask is provided and the + check_mask_outofbounds value is set the masked coordinates are not verified to fall within + the hull. The latter should thus be used with some caution by the user, but can potentially + significantly speed up the mask creation process for axis aligned regions. + """ + self._name = name + self._check_mask_outofbounds = check_mask_outofbounds + self._cached_filled_mask = None + self._vertices = points = np.vstack([b.corners + if hasattr(b, "corners") else [b[0], b[1]] for b in list_hulls]) + self._hull = spat.ConvexHull(points) + if mask is None: + self._mask, self._mask_weights = self.init_mask() + else: + self.sparse_mask = mask + + def invalidate_cached_masks(self): + """ Invalidates the cached masks (sparse or regular) """ + self._cached_filled_mask = None + self._mask, self._mask_weights = self.init_mask() + + def __str__(self): + return ",".join(["({0:d},{1:d})".format(x,y) for (x,y) in self.corners]) + + def init_mask(self): + """ creates a sparse mask of the convex hull of the form (y, x) tuples """ + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + x = np.arange(minx, maxx + 1, 1) #upper limit inclusive + y = np.arange(miny, maxy + 1, 1) + meshgrid = np.meshgrid(y, x) + bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) + + sparse_mask = bounding_mesh if not self._check_mask_outofbounds else \ + [c for c in bounding_mesh if c[::-1] in self] + + mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding + return sparse_mask, mask_weights + + @property + def sprase_mask_weights(self): + """ returns sparse mask weights """ + return self._mask_weights + + @property + def sparse_mask(self): + """ returns a sparse mask (y, x) values of all points in the masked region """ + return self._mask + + @sparse_mask.setter + def sparse_mask(self, mask): + """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ + if not isinstance(mask, list): + raise TypeError("Mask must be list") + if not (hasattr(mask, "__len__") and (len(mask) == 0 or (hasattr(mask[0], "__len__") and len(mask[0]) == 2))): + raise TypeError("Mask must be a sparse mask of 2 element values") + if self._check_mask_outofbounds: + self._mask = copy.deepcopy([c for c in mask if (c[1], c[0]) in self]) + else: + self._mask = copy.deepcopy(mask) + self._mask_weights = np.ones(len(self._mask)) + + @property + def mask(self, dtype=np.float64): + """ Creates a filled rectangular mask grid of size y, x """ + if self._cached_filled_mask is not None: + return self._cached_filled_mask + + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + nx = maxx - minx + 1 # inclusive + ny = maxy - miny + 1 + mesh = np.zeros(nx*ny, dtype=dtype) + if nx==0 or ny==0 or len(self.sparse_mask) == 0: + self._cached_filled_mask = mesh.reshape((ny, nx)) + else: + sparse_mask = np.array(self.sparse_mask) + sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= minx, + sparse_mask[:, 1] <= maxx), + np.logical_and(sparse_mask[:, 0] >= miny, + sparse_mask[:, 0] <= maxy)) + flat_index = (sparse_mask[sel][:, 0] - miny)*nx + (sparse_mask[sel][:, 1] - minx) + mesh[flat_index] = self._mask_weights[sel] + self._cached_filled_mask = mesh.reshape((ny, nx)) + return self._cached_filled_mask + + @classmethod + def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): + """ 2D array containing all values within convex hull + sliced out along axes provided as argument. Portions of sel_region + that are outside of the data_cube is set to oob_value + + assumes the last value of axes is the fastest varying axis + """ + if not isinstance(sel_region, BoundingConvexHull): + raise TypeError("Object passed in is not of type BoundingConvexHull") + if not (hasattr(axes, "__len__") and len(axes) == 2): + raise ValueError("Expected a tupple of axes along which to slice out a region") + axes = sorted(axes) + + lines = np.hstack([sel_region.corners, np.roll(sel_region.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + x = np.arange(minx, maxx + 1, 1) + y = np.arange(miny, maxy + 1, 1) + + pad_left = max(0, 0 - minx) + pad_bottom = max(0, 0 - miny) + pad_right = max(0, maxx - data_cube.shape[axes[1]] + 1) #inclusive of upper limit + pad_top = max(0, maxy - data_cube.shape[axes[0]] + 1) + + if minx > data_cube.shape[axes[0]] or miny > data_cube.shape[axes[1]] or \ + maxy < 0 or maxx < 0: + raise ValueError("Expected a bounding hull that is at least partially within the image") + + # extract data, pad if necessary + slc_data = [slice(None)] * len(data_cube.shape) + for (start, end), axis in zip([(miny + pad_bottom, maxy - pad_top + 1), + (minx + pad_left, maxx - pad_right + 1)], axes): + slc_data[axis] = slice(start, end) + slc_padded = [slice(None)] * len(data_cube.shape) + for (start, end), axis in zip([(pad_bottom, -miny + maxy + 1 - pad_top), + (pad_left, -minx + maxx + 1 - pad_right)], axes): + slc_padded[axis] = slice(start, end) + + selected_data = data_cube[tuple(slc_data)] + new_shape = list(data_cube.shape) + new_shape[axes[0]] = (maxy - miny + 1) + new_shape[axes[1]] = (maxx - minx + 1) + + if any(np.array([pad_left, pad_bottom, pad_right, pad_top]) > 0): + padded_data = np.zeros(tuple(new_shape), dtype=selected_data.dtype) * oob_value + padded_data[tuple(slc_padded)] = selected_data.copy() + else: + padded_data = selected_data.copy() + + # finally apply mask + slc_padded_data = [slice(None)] * len(padded_data.shape) + for (start, end), axis in zip([(0, maxy - miny + 1), #mask starts at origin in the padded image + (0, maxx - minx + 1)], axes): + slc_padded_data[axis] = slice(start, end) + slc_mask = [None] * len(padded_data.shape) + for (start, end), axis in zip([(0, sel_region.mask.shape[0]), #mask starts at origin in the padded image + (0, sel_region.mask.shape[1])], axes): + slc_mask[axis] = slice(start, end) + mask = sel_region.mask.copy() + mask[mask == 0] = oob_value + padded_data[tuple(slc_padded_data)] *= mask[tuple(slc_mask)] + window_extents = [minx, maxx, + miny, maxy] + return padded_data, window_extents + + @classmethod + def normalize_masks(cls, regions, only_overlapped_regions=True): + """ Normalizes region masks for overlapping pixels. This is necessary to properly coadd + overlapping facets. If masks are guarenteed to be initialized to unity (e.g. after + bounding region creation) the user can skip normalizing non-overlapping regions with + flag only_overlapped_regions. + """ + if not all([isinstance(reg, BoundingConvexHull) for reg in regions]): + raise TypeError("Expected a list of bounding convex hulls") + # Implements painters-like algorithm to + # count the number of times a pixel coordinate falls within masks + # The overlapping sections of regions can then be normalized + # For now all regions have equal contribution + allmasks = [] + for reg in regions: + allmasks += list(reg.sparse_mask) if isinstance(reg.sparse_mask, np.ndarray) else reg.sparse_mask + + # flatten for faster comparisons + allmasks = np.array(allmasks) + maxx = np.max(allmasks[:, 1]) + nx = maxx + 1 + allmasks_flatten = allmasks[:, 0] * nx + allmasks[:, 1] + + # now count the number of times a pixel is painted onto + unique_pxls_flatten, paint_count = np.unique(allmasks_flatten, return_counts=True) + paint_count = paint_count.astype(np.float) + + if only_overlapped_regions: + sel = paint_count > 1 + unique_pxls_flatten = unique_pxls_flatten[sel] + paint_count = paint_count[sel] + + # with the reduced number of overlap pixels unflatten + unique_pxls = np.vstack([unique_pxls_flatten // nx, + unique_pxls_flatten % nx]).T + unique_pxls = list(map(tuple, unique_pxls)) + paint_count[...] = 1.0 / paint_count + + # and finally update mask weights + for reg in regions: + reg._cached_filled_mask = None # invalidate + overlap = [x for x in zip(paint_count, unique_pxls) if x[1] in reg.sparse_mask] + for px_pc, px in overlap: + sel = reg.sparse_mask.index(px) if isinstance(reg.sparse_mask, list) else \ + np.all(reg.sparse_mask - px == 0, axis=1) + reg._mask_weights[sel] = px_pc + + @property + def circumference(self): + """ area contained in hull """ + lines = self.edges + return np.sum(np.linalg.norm(lines[:, 1, :] - lines[:, 0, :], axis=1) + 1) + @property + def area(self): + """ area contained in hull """ + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + return 0.5 * np.abs(np.sum([x1*(y2)-(x2)*y1 for x1,y1,x2,y2 in lines])) + 0.5 * self.circumference - 1 + + @property + def name(self): + return self._name + + @name.setter + def name(self, v): + self._name = v + + @property + def corners(self): + """ Returns vertices and guarentees clockwise winding """ + return self._vertices[self._hull.vertices][::-1] + + def normals(self, left = True): + """ return a list of left normals to the hull """ + normals = [] + for i in range(self.corners.shape[0]): + # assuming clockwise winding + j = (i + 1) % self.corners.shape[0] + edge = self.corners[j, :] - self.corners[i, :] + if left: + normals.append((-edge[1], edge[0])) + else: + normals.append((edge[1], -edge[0])) + return np.asarray(normals, dtype=np.double) + + @property + def edges(self): + """ return edge segments of the hull (clockwise wound) """ + edges = [] + for i in range(self.corners.shape[0]): + # assuming clockwise winding + j = (i + 1) % self.corners.shape[0] + edge = tuple([self.corners[i, :], self.corners[j, :]]) + edges.append(edge) + return np.asarray(edges, dtype=np.double) + + @property + def edge_midpoints(self): + """ return edge midpoints of the hull (clockwise wound) """ + edges = self.edges + return np.mean(edges, axis=1) + + @property + def lnormals(self): + """ left normals to the edges of the hull """ + return self.normals(left = True) + + @property + def rnormals(self): + """ right normals to the edges of the hull """ + return self.normals(left=False) + + def overlaps_with(self, other, min_sep_dist=0.5): #less than half a pixel away + """ + Implements the separating lines collision detection theorem + to test whether the hull intersects with 'other' hull + """ + if not isinstance(other, BoundingConvexHull): + raise TypeError("rhs must be a BoundingConvexHull") + + # get the projection axes + normals = np.vstack([self.lnormals, other.lnormals]) + norms = np.linalg.norm(normals, axis=1) + normals = normals / norms[None, 2] + + # compute vectors to corners from origin + vecs_reg1 = self.corners + vecs_reg2 = other.corners + + # compute projections onto normals + for ni, n in enumerate(normals): + projs = np.dot(vecs_reg1, n.T) + minproj_reg1 = np.min(projs) + maxproj_reg1 = np.max(projs) + projs = np.dot(vecs_reg2, n.T) + minproj_reg2 = np.min(projs) + maxproj_reg2 = np.max(projs) + if minproj_reg2 - maxproj_reg1 > min_sep_dist or minproj_reg1 - maxproj_reg2 > min_sep_dist: + return False + return True + + @property + def centre(self, integral=True): + """ Barycentre of hull """ + if integral: + def rnd(x): + return int(np.floor(x) if x >= 0 else np.ceil(x)) + return [rnd(x) for x in np.mean(self._vertices, axis=0)] + else: + return np.mean(self._vertices, axis=0) + + def __contains__(self, s, tolerance=0.5): #less than half a pixel away + """ tests whether a point s(x,y) is in the convex hull """ + # there are three cases to consider + # CASE 1: + # scalar projection between all inner pointing right normals (clockwise winding) + # and the point must be positive if the point were to lie inside + # the region (true) + # CASE 2: + # point is on an edge - the scalar projection onto the axis is 0 for that edge + # and greater than 0 for the other edges (true) + # CASE 3: + # it is outside (false) + x, y = s + isin = True + normals = self.rnormals + xyvec = np.array([x, y])[None, :] - np.array(self.corners) + + dot = np.einsum("ij,ij->i", normals, xyvec) + return np.all(dot > -tolerance) + +class BoundingBox(BoundingConvexHull): + def __init__(self, xl, xu, yl, yu, name="unnamed", mask=None, **kwargs): + if not all([isinstance(x, (int, np.int64, np.int32, np.int16)) for x in [xl, xu, yl, yu]]): + raise ValueError("Box limits must be integers") + self.__xnpx = abs(xu - xl + 1) #inclusive of the upper pixel + self.__ynpx = abs(yu - yl + 1) + BoundingConvexHull.__init__(self, + [[xl,yl],[xl,yu],[xu,yu],[xu,yl]], + name, + mask=mask, + **kwargs) + + def init_mask(self): + """ creates a sparse mask of the convex hull of the form (y, x) tuples """ + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + x = np.arange(minx, maxx + 1, 1) #upper limit inclusive + y = np.arange(miny, maxy + 1, 1) + meshgrid = np.meshgrid(y, x) + bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) + + sparse_mask = np.asarray(bounding_mesh) # by default for a BB region the mask is always going to be the entire region + + mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding + return sparse_mask, mask_weights + + def __contains__(self, s): + """ tests whether a point s(x,y) is in the box""" + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + return s[0] >= minx and s[0] <= maxx and s[1] >= miny and s[1] <= maxy + + @property + def box_npx(self): + return (self.__xnpx, self.__ynpx) + + @property + def sparse_mask(self): + """ returns a sparse mask (y, x) values of all points in the masked region """ + return self._mask + + @sparse_mask.setter + def sparse_mask(self, mask): + """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ + if not isinstance(mask, list) and not isinstance(mask, np.ndarray): + raise TypeError("Mask must be list") + if not (hasattr(mask, "__len__") and (len(mask) == 0 or (hasattr(mask[0], "__len__") and len(mask[0]) == 2))): + raise TypeError("Mask must be a sparse mask of 2 element values") + if mask == []: + self._mask = [] + else: + lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + nx = maxx - minx + 1 # inclusive + ny = maxy - miny + 1 + sparse_mask = np.asarray(mask) + sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= minx, + sparse_mask[:, 1] <= maxx), + np.logical_and(sparse_mask[:, 0] >= miny, + sparse_mask[:, 0] <= maxy)) + self._mask = sparse_mask[sel] + self._mask_weights = np.ones(len(self._mask)) + + @classmethod + def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np.float64, **kwargs): + """ Projects individial regions back onto a single contiguous cube """ + if not (hasattr(regional_data_list, "__len__") and hasattr(regions_list, "__len__") and \ + len(regions_list) == len(regional_data_list)): + raise TypeError("Region data list and regions lists must be lists of equal length") + if not all([isinstance(x, np.ndarray) for x in regional_data_list]): + raise TypeError("Region data list must be a list of ndarrays") + if not all([isinstance(x, BoundingBox) for x in regions_list]): + raise TypeError("Region list must be a list of Axis Aligned Bounding Boxes") + if regions_list == []: + return np.empty((0)) + if not all([reg.ndim == regional_data_list[0].ndim for reg in regional_data_list]): + raise ValueError("All data cubes must be of equal dimension") + axes = tuple(sorted(axes)) + + minx = np.min([np.min(f.corners[:, 0]) for f in regions_list]) + maxx = np.max([np.max(f.corners[:, 0]) for f in regions_list]) + miny = np.min([np.min(f.corners[:, 1]) for f in regions_list]) + maxy = np.max([np.max(f.corners[:, 1]) for f in regions_list]) + npxx = maxx - minx + 1 + npxy = maxy - miny + 1 + global_offsetx = -minx #-min(0, minx) + global_offsety = -miny #-min(0, miny) + + projected_image_size = list(regional_data_list[0].shape) + projected_image_size[axes[0]] = npxy + projected_image_size[axes[1]] = npxx + stitched_img = np.zeros(tuple(projected_image_size), dtype=dtype) + + combined_mask = [] + for f, freg in zip(regional_data_list, regions_list): + f[np.isnan(f)] = 0 + xl = max(0, global_offsetx+np.min(freg.corners[:, 0])) + xu = min(global_offsetx+np.max(freg.corners[:, 0]) + 1, npxx) + yl = max(0, global_offsety+np.min(freg.corners[:, 1])) + yu = min(global_offsety+np.max(freg.corners[:, 1]) + 1, npxy) + fnx = xu - xl + 1 # inclusive + fny = yu - yl + 1 # inclusive + if f.shape[axes[0]] != fny - 1 or f.shape[axes[1]] != fnx - 1: + raise ValueError("One or more bounding box descriptors does not match shape of corresponding data cubes") + slc_data = [slice(None)] * len(stitched_img.shape) + for (start, end), axis in zip([(yl, yu), (xl, xu)], axes): + slc_data[axis] = slice(start, end) + + stitched_img[tuple(slc_data)] += f + combined_mask += list(freg.sparse_mask) + + return stitched_img, BoundingBox(minx, maxx, miny, maxy, mask=combined_mask, **kwargs) + +######################################################################## +## Factories +######################################################################## + +class BoundingBoxFactory(object): + @classmethod + def AxisAlignedBoundingBox(cls, convex_hull_object, square=False, enforce_odd=True, **kwargs): + """ Constructs an axis aligned bounding box around convex hull """ + if not isinstance(convex_hull_object, BoundingConvexHull): + raise TypeError("Convex hull object passed in constructor is not of type BoundingConvexHull") + if square: + nx = np.max(convex_hull_object.corners[:, 0]) - np.min(convex_hull_object.corners[:, 0]) + 1 #inclusive + ny = np.max(convex_hull_object.corners[:, 1]) - np.min(convex_hull_object.corners[:, 1]) + 1 #inclusive + boxdiam = max(nx, ny) + boxrad = boxdiam // 2 + cx, cy = convex_hull_object.centre + xl = cx - boxrad + xu = cx + boxdiam - boxrad - 1 + yl = cy - boxrad + yu = cy + boxdiam - boxrad - 1 + else: + xl = np.min(convex_hull_object.corners[:, 0]) + xu = np.max(convex_hull_object.corners[:, 0]) + yl = np.min(convex_hull_object.corners[:, 1]) + yu = np.max(convex_hull_object.corners[:, 1]) + + xu += (xu - xl) % 2 if enforce_odd else 0 + yu += (yu - yl) % 2 if enforce_odd else 0 + + return BoundingBox(xl, xu, yl, yu, + convex_hull_object.name, + mask=convex_hull_object.sparse_mask, + **kwargs) + + @classmethod + def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): + """ Split a axis-aligned bounding box into smaller boxes """ + if not isinstance(bounding_box_object, BoundingBox): + raise TypeError("Expected bounding box object") + if not (isinstance(nsubboxes, int) and nsubboxes >= 1): + raise ValueError("nsubboxes must be integral type and be 1 or more") + xl = np.min(bounding_box_object.corners[:, 0]) + xu = np.max(bounding_box_object.corners[:, 0]) + yl = np.min(bounding_box_object.corners[:, 1]) + yu = np.max(bounding_box_object.corners[:, 1]) + + # construct a nonregular meshgrid bound to xu and yu + x = xl + np.arange(0, nsubboxes + 1) * int(np.ceil((xu - xl + 1) / float(nsubboxes))) + y = yl + np.arange(0, nsubboxes + 1) * int(np.ceil((yu - yl + 1) / float(nsubboxes))) + xx, yy = np.meshgrid(x, y) + + # split into boxes + xls = xx[0:-1, 0:-1].copy() + xus = xx[1:, 1:].copy() + yls = yy[0:-1, 0:-1].copy() + yus = yy[1:, 1:].copy() + + # make sure no boxes overlap + xus = xus - 1 + yus = yus - 1 + + # clamp the final coordinate to the upper end (may result in rectanglular box at the end) + xus[:, -1] = max(xu, min(xus[0, -1], xu)) + yus[-1, :] = max(yu, min(yus[-1, 0], yu)) + + #coordinates for all the contained boxes, anti-clockwise wound + xls = xls.ravel() + yls = yls.ravel() + xus = xus.ravel() + yus = yus.ravel() + bl = list(zip(xls, yls)) + br = list(zip(xus, yls)) + ur = list(zip(xus, yus)) + ul = list(zip(xls, yus)) + contained_boxes = list(zip(bl, br, ur, ul)) + + #finally create bbs for each of the contained boxes with the mask + #chopped up between the boxes by the convex hull initializer + new_regions = [BoundingBox(bl[0], br[0], bl[1], ul[1], + bounding_box_object.name, + mask=bounding_box_object.sparse_mask, + **kwargs) + for bl, br, ur, ul in contained_boxes] + + return new_regions + + @classmethod + def PadBox(cls, bounding_box_object, desired_nx, desired_ny, **kwargs): + """ Creates a box with a padded border around a axis-aligned bounding box """ + if not isinstance(bounding_box_object, BoundingBox): + raise TypeError("Expected bounding box object") + nx, ny = bounding_box_object.box_npx + if desired_nx - nx < 0 or desired_ny - ny < 0: + raise ValueError("Padded box must be bigger than original box") + pad_left = desired_nx // 2 + pad_right = desired_nx - pad_left - 1 + pad_bottom = desired_ny // 2 + pad_top = desired_ny - pad_bottom - 1 + cx, cy = bounding_box_object.centre + xl = cx - pad_left + xu = cx + pad_right + yl = cy - pad_bottom + yu = cy + pad_top + return BoundingBox(xl, xu, yl, yu, + bounding_box_object.name, + mask=bounding_box_object.sparse_mask, + **kwargs) #mask unchanged in the new shape, border frame discarded + diff --git a/africanus/testing/tests/test_geometry.py b/africanus/testing/tests/test_geometry.py new file mode 100644 index 000000000..315a436ac --- /dev/null +++ b/africanus/testing/tests/test_geometry.py @@ -0,0 +1,219 @@ +import numpy as np +import scipy.stats as sstats +import scipy.signal as ssig +import scipy.spatial as spat +import copy +from africanus.linalg.geometry import (BoundingConvexHull, + BoundingBox, + BoundingBoxFactory) +DEBUG=True +def test_hull_construction(): + # test case 1 + vals = np.array([[50, 60], [20, 40], [-74, 50], [-95, +10], [20, 60]]) + bh = BoundingConvexHull(vals) + mask = bh.mask + assert mask.shape == (np.max(vals[:, 1]) - np.min(vals[:, 1]) + 1, np.max(vals[:, 0]) - np.min(vals[:, 0]) + 1) + assert np.abs(mask.sum() - bh.area) / bh.area < 0.05 # integral mask area needs to be close to true area + normalized_normals = bh.rnormals / np.linalg.norm(bh.rnormals, axis=1)[:, None] + # test case 2 + for e, n in zip(bh.edges, normalized_normals): + edge_vec = e[1] - e[0] + assert np.all(np.abs(np.dot(edge_vec, n)) < 1.0e-8) + + # test case 3 + valsextract = np.array([[-10, 120], [90, 268], [293, 110],[40, -30]]) + bh_extract = BoundingConvexHull(valsextract) + sinc_npx = 255 + sinc = np.sinc(np.linspace(-7, 7, sinc_npx)) + sinc2d = np.outer(sinc, sinc).reshape((1, 1, sinc_npx, sinc_npx)) + extracted_data, extracted_window_extents = BoundingConvexHull.regional_data(bh_extract, sinc2d, oob_value=np.nan) + assert extracted_window_extents == [-10, 293, -30, 268] + sparse_mask = np.array(bh_extract.sparse_mask) + lines = np.hstack([bh_extract.corners, np.roll(bh_extract.corners, -1, axis=0)]) + minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= 0, + sparse_mask[:, 1] < 255), + np.logical_and(sparse_mask[:, 0] >= 0, + sparse_mask[:, 0] < 255)) + + flat_index = (sparse_mask[sel][:, 0])*sinc_npx + (sparse_mask[sel][:, 1]) + sinc_integral = np.sum(sinc2d.ravel()[flat_index]) + assert np.abs(sinc_integral - np.nansum(extracted_data.ravel())) < 1.0e-8 + v = np.nanargmax(extracted_data) + vx = v % extracted_data.shape[3]; vy = v // extracted_data.shape[3] + cextracted = (extracted_window_extents[0] + vx, + extracted_window_extents[2] + vy) + v = np.nanargmax(sinc2d) + sincvx = v % sinc_npx; sincvy = v // sinc_npx + csinc = tuple([sincvx, sincvy]) + assert csinc == cextracted + + # test case 4 + vals2 = np.array([[-20, -120], [0, 60], [40, -60]]) + vals3 = np.array([[-20, 58], [-40, 80], [20, 100]]) + bh2 = BoundingConvexHull(vals2) + bh3 = BoundingConvexHull(vals3) + assert bh.overlaps_with(bh2) + assert not bh.overlaps_with(bh3) + assert not bh2.overlaps_with(bh3) + + # test case 5 + assert (-1000, -1000) not in bh + assert (30, 0) not in bh + assert (0, 0) not in bh + assert (-40, 30) in bh + + # test case 6 + bb = BoundingBox(-14, 20, 30, 49) + assert bb.centre == [3, 39] + assert bb.box_npx == (35, 20) + assert bb.mask.shape == bb.box_npx[::-1] + assert bb.area == 35 * 20 + + assert np.sum(bb.mask) == bb.area + assert (-15, 35) not in bb + assert (0, 35) in bb + + bb2 = BoundingBoxFactory.AxisAlignedBoundingBox(bb) #enforce odd + assert bb2.box_npx == (35, 21) + assert bb2.area == 35 * 21 + assert (bb.sparse_mask == bb2.sparse_mask).all() + assert (-15, 35) not in bb2 + assert (0, 35) in bb2 + + bb3 = BoundingBoxFactory.AxisAlignedBoundingBox(bb, square=True) #enforce odd + assert bb3.box_npx[0] == bb3.box_npx[1] + assert bb3.box_npx[0] % 2 == 1 #enforce odd + assert bb3.area == bb3.box_npx[0]**2 + assert (bb.sparse_mask == bb3.sparse_mask).all() + assert (-15, 35) not in bb2 + assert (0, 35) in bb2 + + # test case 7 + bb4s = BoundingBoxFactory.SplitBox(bb, nsubboxes=3) + assert len(bb4s) == 9 + xlims = [(np.min(c.corners[:, 0]), np.max(c.corners[:, 0])) for c in bb4s][0:3] + ylims = [(np.min(c.corners[:, 1]), np.max(c.corners[:, 1])) for c in bb4s][0::3] + assert np.all(xlims == np.array([(-14, -3), (-2, 9), (10, 20)])) + assert np.all(ylims == np.array([(30, 36), (37, 43), (44, 49)])) + assert np.sum([b.area for b in bb4s]) == bb.area + + for bb4 in bb4s: + assert bb4.area == np.sum(bb4.mask) + + # test case 8 + bb5 = BoundingBox(-14, 20, 30, 50) + assert bb5.box_npx == (35, 21) + bb6 = BoundingBoxFactory.PadBox(bb5, 41, 27) + assert bb6.box_npx == (41, 27) + assert bb5.centre == bb6.centre + assert np.sum(bb5.mask) == np.sum(bb6.mask) + bb7s = list(map(lambda x: BoundingBoxFactory.PadBox(x, 17, 11), bb4s)) + assert all([b.box_npx == (17, 11) for b in bb7s]) + assert np.sum([np.sum(b.mask) for b in bb7s]) == np.sum([np.sum(b.mask) for b in bb4s]) + + # test case 9 + facet_regions = list(map(lambda f: BoundingBoxFactory.PadBox(f, 63, 63), + BoundingBoxFactory.SplitBox(BoundingBoxFactory.AxisAlignedBoundingBox(bh_extract), nsubboxes=5))) + facets = list(map(lambda pf: BoundingConvexHull.regional_data(pf, sinc2d, oob_value=np.nan), + facet_regions)) + stitched_image, stitched_region = BoundingBox.project_regions([f[0] for f in facets], facet_regions) + assert np.abs(sinc_integral - np.nansum([np.nansum(f[0]) for f in facets])) < 1.0e-8 + assert np.abs(sinc_integral - np.sum(stitched_image)) < 1.0e-8 + v = np.argmax(stitched_image) + vx = v % stitched_image.shape[3]; vy = v // stitched_image.shape[3] + cstitched = (np.min(stitched_region.corners[:, 0]) + vx, np.min(stitched_region.corners[:, 1]) + vy) + assert cstitched == csinc + + # test case 10 + olap_box1 = BoundingBox(110, 138, 110, 135) + olap_box2 = BoundingBox(115, 150, 109, 150) + olap_box3 = BoundingBox(125, 130, 125, 130) + BoundingConvexHull.normalize_masks([olap_box1, olap_box2, olap_box3]) + ext1 = BoundingConvexHull.regional_data(olap_box1, sinc2d)[0] + ext2 = BoundingConvexHull.regional_data(olap_box2, sinc2d)[0] + ext3 = BoundingConvexHull.regional_data(olap_box3, sinc2d)[0] + olaps_stitched_image, olaps_stitched_region = BoundingBox.project_regions([ext1, ext2, ext3], + [olap_box1, olap_box2, olap_box3]) + v = np.nanargmax(olaps_stitched_image) + vx = v % olaps_stitched_image.shape[3]; vy = v // olaps_stitched_image.shape[3] + cstitched_olap = (np.min(olaps_stitched_region.corners[:, 0]) + vx, + np.min(olaps_stitched_region.corners[:, 1]) + vy) + assert cstitched_olap == csinc + assert np.abs(1.0 - np.nanmax(olaps_stitched_image)) < 1.0e-8 + + # visual inspection + if DEBUG: + from matplotlib import pyplot as plt + plt.figure(figsize=(7, 2.5)) + plt.title("Winding, normals and masking check") + for h in [bh, bh2, bh3]: + for ei, e in enumerate(h.edges): + plt.plot(e[:, 0], e[:, 1], "r--") + plt.text(e[0, 0], e[0, 1], str(ei)) + + plt.plot(bh.edge_midpoints[:, 0], bh.edge_midpoints[:, 1], "ko") + for e, n in zip(bh.edge_midpoints, normalized_normals): + p0 = e + p = e + n*6 + plt.plot([p0[0], p[0]], [p0[1], p[1]], "b--", lw=2) + + plt.scatter(vals[:, 0], vals[:, 1]) + plt.imshow(mask, extent=[np.min(vals[:, 0]), np.max(vals[:, 0]), np.max(vals[:, 1]), np.min(vals[:, 1])]) + + plt.grid(True) + plt.savefig("/tmp/winding.png") + + plt.figure(figsize=(7, 2.5)) + plt.title("Data extraction check (global)") + for h in [bh_extract]: + for ei, e in enumerate(h.edges): + plt.plot(e[:, 0], e[:, 1], "r--") + plt.imshow(sinc2d[0, 0, :, :], extent=[0, sinc_npx, sinc_npx, 0]) + plt.grid(True) + plt.savefig("/tmp/extract_global.png") + + plt.figure(figsize=(7, 2.5)) + plt.title("Data extraction check (local)") + for h in [bh_extract]: + for ei, e in enumerate(h.edges): + plt.plot(e[:, 0], e[:, 1], "r--") + plt.imshow(extracted_data[0, 0, :, :], + extent=[extracted_window_extents[0], extracted_window_extents[1], + extracted_window_extents[3], extracted_window_extents[2]]) + plt.savefig("/tmp/extract_local.png") + + plt.figure(figsize=(7, 2.5)) + plt.title("Faceting check") + for h in [bh_extract]: + for ei, e in enumerate(h.edges): + plt.plot(e[:, 0], e[:, 1], "r--") + for f in facet_regions: + for ei, e in enumerate(f.edges): + plt.plot(e[:, 0], e[:, 1], "co--") + + + plt.imshow(stitched_image[0, 0, :, :], + extent=[np.min(stitched_region.corners[:, 0]), np.max(stitched_region.corners[:, 0]), + np.max(stitched_region.corners[:, 1]), np.min(stitched_region.corners[:, 1])]) + plt.savefig("/tmp/facet.png") + + plt.figure(figsize=(7, 2.5)) + plt.title("Overlapping faceting check") + for f in [olap_box1, olap_box2, olap_box3]: + for ei, e in enumerate(f.edges): + plt.plot(e[:, 0], e[:, 1], "co--") + + + plt.imshow(olaps_stitched_image[0, 0, :, :], + extent=[np.min(olaps_stitched_region.corners[:, 0]), np.max(olaps_stitched_region.corners[:, 0]), + np.max(olaps_stitched_region.corners[:, 1]), np.min(olaps_stitched_region.corners[:, 1])]) + plt.xlim((np.min(olaps_stitched_region.corners[:, 0]) - 15, + np.max(olaps_stitched_region.corners[:, 0]) + 15)) + plt.ylim((np.min(olaps_stitched_region.corners[:, 1]) - 15, + np.max(olaps_stitched_region.corners[:, 1]) + 15)) + plt.savefig("/tmp/overlap_facet.png") + +if __name__ == "__main__": + test_hull_construction() From b9183edf10d539256a8fc36bc039d8c6e0b59752 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:10:25 +0200 Subject: [PATCH 02/10] Handle optional scipy requirements --- africanus/linalg/geometry.py | 86 ++++++++++++++++++++---------------- 1 file changed, 48 insertions(+), 38 deletions(-) diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py index b4da4d92c..7769c7b3d 100644 --- a/africanus/linalg/geometry.py +++ b/africanus/linalg/geometry.py @@ -1,10 +1,19 @@ import numpy as np -import scipy.stats as sstats -import scipy.signal as ssig -import scipy.spatial as spat import copy - import time + +try: + # import scipy.stats as sstats + # import scipy.signal as ssig + import scipy.spatial as spat +except ImportError as e: + opt_import_err = e +else: + opt_import_err = None + +from africanus.util.requirements import requires_optional + + def timeit(method): def timed(*args, **kw): ts = time.time() @@ -16,15 +25,16 @@ def timed(*args, **kw): else: print(('%r %2.2f ms' % \ (method.__name__, (te - ts) * 1000))) - + return result return timed +@requires_optional("scipy.stats", opt_import_err) class BoundingConvexHull(object): def __init__(self, list_hulls, name="unnamed", mask = None, check_mask_outofbounds=True): """ Initializes a bounding convex hull around a list of bounding convex hulls or series of points. A unity-weighted mask is computed for the region that falls within this convex hull - if a mask of (y, x) coordinates is not provided. Otherwise if a mask is provided and the + if a mask of (y, x) coordinates is not provided. Otherwise if a mask is provided and the check_mask_outofbounds value is set the masked coordinates are not verified to fall within the hull. The latter should thus be used with some caution by the user, but can potentially significantly speed up the mask creation process for axis aligned regions. @@ -37,7 +47,7 @@ def __init__(self, list_hulls, name="unnamed", mask = None, check_mask_outofboun self._hull = spat.ConvexHull(points) if mask is None: self._mask, self._mask_weights = self.init_mask() - else: + else: self.sparse_mask = mask def invalidate_cached_masks(self): @@ -47,7 +57,7 @@ def invalidate_cached_masks(self): def __str__(self): return ",".join(["({0:d},{1:d})".format(x,y) for (x,y) in self.corners]) - + def init_mask(self): """ creates a sparse mask of the convex hull of the form (y, x) tuples """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) @@ -57,10 +67,10 @@ def init_mask(self): y = np.arange(miny, maxy + 1, 1) meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) - + sparse_mask = bounding_mesh if not self._check_mask_outofbounds else \ [c for c in bounding_mesh if c[::-1] in self] - + mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding return sparse_mask, mask_weights @@ -86,7 +96,7 @@ def sparse_mask(self, mask): else: self._mask = copy.deepcopy(mask) self._mask_weights = np.ones(len(self._mask)) - + @property def mask(self, dtype=np.float64): """ Creates a filled rectangular mask grid of size y, x """ @@ -99,7 +109,7 @@ def mask(self, dtype=np.float64): nx = maxx - minx + 1 # inclusive ny = maxy - miny + 1 mesh = np.zeros(nx*ny, dtype=dtype) - if nx==0 or ny==0 or len(self.sparse_mask) == 0: + if nx==0 or ny==0 or len(self.sparse_mask) == 0: self._cached_filled_mask = mesh.reshape((ny, nx)) else: sparse_mask = np.array(self.sparse_mask) @@ -114,7 +124,7 @@ def mask(self, dtype=np.float64): @classmethod def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): - """ 2D array containing all values within convex hull + """ 2D array containing all values within convex hull sliced out along axes provided as argument. Portions of sel_region that are outside of the data_cube is set to oob_value @@ -143,11 +153,11 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): # extract data, pad if necessary slc_data = [slice(None)] * len(data_cube.shape) - for (start, end), axis in zip([(miny + pad_bottom, maxy - pad_top + 1), + for (start, end), axis in zip([(miny + pad_bottom, maxy - pad_top + 1), (minx + pad_left, maxx - pad_right + 1)], axes): slc_data[axis] = slice(start, end) slc_padded = [slice(None)] * len(data_cube.shape) - for (start, end), axis in zip([(pad_bottom, -miny + maxy + 1 - pad_top), + for (start, end), axis in zip([(pad_bottom, -miny + maxy + 1 - pad_top), (pad_left, -minx + maxx + 1 - pad_right)], axes): slc_padded[axis] = slice(start, end) @@ -174,10 +184,10 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): mask = sel_region.mask.copy() mask[mask == 0] = oob_value padded_data[tuple(slc_padded_data)] *= mask[tuple(slc_mask)] - window_extents = [minx, maxx, + window_extents = [minx, maxx, miny, maxy] return padded_data, window_extents - + @classmethod def normalize_masks(cls, regions, only_overlapped_regions=True): """ Normalizes region masks for overlapping pixels. This is necessary to properly coadd @@ -194,7 +204,7 @@ def normalize_masks(cls, regions, only_overlapped_regions=True): allmasks = [] for reg in regions: allmasks += list(reg.sparse_mask) if isinstance(reg.sparse_mask, np.ndarray) else reg.sparse_mask - + # flatten for faster comparisons allmasks = np.array(allmasks) maxx = np.max(allmasks[:, 1]) @@ -204,12 +214,12 @@ def normalize_masks(cls, regions, only_overlapped_regions=True): # now count the number of times a pixel is painted onto unique_pxls_flatten, paint_count = np.unique(allmasks_flatten, return_counts=True) paint_count = paint_count.astype(np.float) - + if only_overlapped_regions: sel = paint_count > 1 unique_pxls_flatten = unique_pxls_flatten[sel] paint_count = paint_count[sel] - + # with the reduced number of overlap pixels unflatten unique_pxls = np.vstack([unique_pxls_flatten // nx, unique_pxls_flatten % nx]).T @@ -239,7 +249,7 @@ def area(self): @property def name(self): return self._name - + @name.setter def name(self, v): self._name = v @@ -288,10 +298,10 @@ def lnormals(self): def rnormals(self): """ right normals to the edges of the hull """ return self.normals(left=False) - + def overlaps_with(self, other, min_sep_dist=0.5): #less than half a pixel away - """ - Implements the separating lines collision detection theorem + """ + Implements the separating lines collision detection theorem to test whether the hull intersects with 'other' hull """ if not isinstance(other, BoundingConvexHull): @@ -332,7 +342,7 @@ def __contains__(self, s, tolerance=0.5): #less than half a pixel away """ tests whether a point s(x,y) is in the convex hull """ # there are three cases to consider # CASE 1: - # scalar projection between all inner pointing right normals (clockwise winding) + # scalar projection between all inner pointing right normals (clockwise winding) # and the point must be positive if the point were to lie inside # the region (true) # CASE 2: @@ -342,9 +352,9 @@ def __contains__(self, s, tolerance=0.5): #less than half a pixel away # it is outside (false) x, y = s isin = True - normals = self.rnormals + normals = self.rnormals xyvec = np.array([x, y])[None, :] - np.array(self.corners) - + dot = np.einsum("ij,ij->i", normals, xyvec) return np.all(dot > -tolerance) @@ -369,9 +379,9 @@ def init_mask(self): y = np.arange(miny, maxy + 1, 1) meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) - + sparse_mask = np.asarray(bounding_mesh) # by default for a BB region the mask is always going to be the entire region - + mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding return sparse_mask, mask_weights @@ -458,7 +468,7 @@ def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np slc_data = [slice(None)] * len(stitched_img.shape) for (start, end), axis in zip([(yl, yu), (xl, xu)], axes): slc_data[axis] = slice(start, end) - + stitched_img[tuple(slc_data)] += f combined_mask += list(freg.sparse_mask) @@ -469,7 +479,7 @@ def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np ######################################################################## class BoundingBoxFactory(object): - @classmethod + @classmethod def AxisAlignedBoundingBox(cls, convex_hull_object, square=False, enforce_odd=True, **kwargs): """ Constructs an axis aligned bounding box around convex hull """ if not isinstance(convex_hull_object, BoundingConvexHull): @@ -495,7 +505,7 @@ def AxisAlignedBoundingBox(cls, convex_hull_object, square=False, enforce_odd=Tr return BoundingBox(xl, xu, yl, yu, convex_hull_object.name, - mask=convex_hull_object.sparse_mask, + mask=convex_hull_object.sparse_mask, **kwargs) @classmethod @@ -509,7 +519,7 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): xu = np.max(bounding_box_object.corners[:, 0]) yl = np.min(bounding_box_object.corners[:, 1]) yu = np.max(bounding_box_object.corners[:, 1]) - + # construct a nonregular meshgrid bound to xu and yu x = xl + np.arange(0, nsubboxes + 1) * int(np.ceil((xu - xl + 1) / float(nsubboxes))) y = yl + np.arange(0, nsubboxes + 1) * int(np.ceil((yu - yl + 1) / float(nsubboxes))) @@ -526,9 +536,9 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): yus = yus - 1 # clamp the final coordinate to the upper end (may result in rectanglular box at the end) - xus[:, -1] = max(xu, min(xus[0, -1], xu)) + xus[:, -1] = max(xu, min(xus[0, -1], xu)) yus[-1, :] = max(yu, min(yus[-1, 0], yu)) - + #coordinates for all the contained boxes, anti-clockwise wound xls = xls.ravel() yls = yls.ravel() @@ -544,7 +554,7 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): #chopped up between the boxes by the convex hull initializer new_regions = [BoundingBox(bl[0], br[0], bl[1], ul[1], bounding_box_object.name, - mask=bounding_box_object.sparse_mask, + mask=bounding_box_object.sparse_mask, **kwargs) for bl, br, ur, ul in contained_boxes] @@ -568,7 +578,7 @@ def PadBox(cls, bounding_box_object, desired_nx, desired_ny, **kwargs): yl = cy - pad_bottom yu = cy + pad_top return BoundingBox(xl, xu, yl, yu, - bounding_box_object.name, - mask=bounding_box_object.sparse_mask, + bounding_box_object.name, + mask=bounding_box_object.sparse_mask, **kwargs) #mask unchanged in the new shape, border frame discarded From 86b9294ab2bfcac32051e26c3921ce430b248b4c Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:14:04 +0200 Subject: [PATCH 03/10] Move geometry test --- africanus/{testing/tests => linalg/test}/test_geometry.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename africanus/{testing/tests => linalg/test}/test_geometry.py (100%) diff --git a/africanus/testing/tests/test_geometry.py b/africanus/linalg/test/test_geometry.py similarity index 100% rename from africanus/testing/tests/test_geometry.py rename to africanus/linalg/test/test_geometry.py From 8ef62bfdbe261ce751cb44d7e02feff91dd3683d Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:14:32 +0200 Subject: [PATCH 04/10] Make debug a parametrization --- africanus/linalg/test/test_geometry.py | 56 +++++++++++++------------- 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/africanus/linalg/test/test_geometry.py b/africanus/linalg/test/test_geometry.py index 315a436ac..f4f90e7d3 100644 --- a/africanus/linalg/test/test_geometry.py +++ b/africanus/linalg/test/test_geometry.py @@ -1,17 +1,15 @@ import numpy as np -import scipy.stats as sstats -import scipy.signal as ssig -import scipy.spatial as spat -import copy +import pytest + from africanus.linalg.geometry import (BoundingConvexHull, BoundingBox, BoundingBoxFactory) -DEBUG=True -def test_hull_construction(): +@pytest.mark.parametrize("debug", [False]) +def test_hull_construction(debug): # test case 1 vals = np.array([[50, 60], [20, 40], [-74, 50], [-95, +10], [20, 60]]) bh = BoundingConvexHull(vals) - mask = bh.mask + mask = bh.mask assert mask.shape == (np.max(vals[:, 1]) - np.min(vals[:, 1]) + 1, np.max(vals[:, 0]) - np.min(vals[:, 0]) + 1) assert np.abs(mask.sum() - bh.area) / bh.area < 0.05 # integral mask area needs to be close to true area normalized_normals = bh.rnormals / np.linalg.norm(bh.rnormals, axis=1)[:, None] @@ -36,9 +34,9 @@ def test_hull_construction(): sparse_mask[:, 1] < 255), np.logical_and(sparse_mask[:, 0] >= 0, sparse_mask[:, 0] < 255)) - + flat_index = (sparse_mask[sel][:, 0])*sinc_npx + (sparse_mask[sel][:, 1]) - sinc_integral = np.sum(sinc2d.ravel()[flat_index]) + sinc_integral = np.sum(sinc2d.ravel()[flat_index]) assert np.abs(sinc_integral - np.nansum(extracted_data.ravel())) < 1.0e-8 v = np.nanargmax(extracted_data) vx = v % extracted_data.shape[3]; vy = v // extracted_data.shape[3] @@ -46,9 +44,9 @@ def test_hull_construction(): extracted_window_extents[2] + vy) v = np.nanargmax(sinc2d) sincvx = v % sinc_npx; sincvy = v // sinc_npx - csinc = tuple([sincvx, sincvy]) + csinc = tuple([sincvx, sincvy]) assert csinc == cextracted - + # test case 4 vals2 = np.array([[-20, -120], [0, 60], [40, -60]]) vals3 = np.array([[-20, 58], [-40, 80], [20, 100]]) @@ -70,7 +68,7 @@ def test_hull_construction(): assert bb.box_npx == (35, 20) assert bb.mask.shape == bb.box_npx[::-1] assert bb.area == 35 * 20 - + assert np.sum(bb.mask) == bb.area assert (-15, 35) not in bb assert (0, 35) in bb @@ -114,7 +112,7 @@ def test_hull_construction(): assert np.sum([np.sum(b.mask) for b in bb7s]) == np.sum([np.sum(b.mask) for b in bb4s]) # test case 9 - facet_regions = list(map(lambda f: BoundingBoxFactory.PadBox(f, 63, 63), + facet_regions = list(map(lambda f: BoundingBoxFactory.PadBox(f, 63, 63), BoundingBoxFactory.SplitBox(BoundingBoxFactory.AxisAlignedBoundingBox(bh_extract), nsubboxes=5))) facets = list(map(lambda pf: BoundingConvexHull.regional_data(pf, sinc2d, oob_value=np.nan), facet_regions)) @@ -134,17 +132,17 @@ def test_hull_construction(): ext1 = BoundingConvexHull.regional_data(olap_box1, sinc2d)[0] ext2 = BoundingConvexHull.regional_data(olap_box2, sinc2d)[0] ext3 = BoundingConvexHull.regional_data(olap_box3, sinc2d)[0] - olaps_stitched_image, olaps_stitched_region = BoundingBox.project_regions([ext1, ext2, ext3], + olaps_stitched_image, olaps_stitched_region = BoundingBox.project_regions([ext1, ext2, ext3], [olap_box1, olap_box2, olap_box3]) v = np.nanargmax(olaps_stitched_image) vx = v % olaps_stitched_image.shape[3]; vy = v // olaps_stitched_image.shape[3] - cstitched_olap = (np.min(olaps_stitched_region.corners[:, 0]) + vx, + cstitched_olap = (np.min(olaps_stitched_region.corners[:, 0]) + vx, np.min(olaps_stitched_region.corners[:, 1]) + vy) assert cstitched_olap == csinc assert np.abs(1.0 - np.nanmax(olaps_stitched_image)) < 1.0e-8 - + # visual inspection - if DEBUG: + if debug: from matplotlib import pyplot as plt plt.figure(figsize=(7, 2.5)) plt.title("Winding, normals and masking check") @@ -152,14 +150,14 @@ def test_hull_construction(): for ei, e in enumerate(h.edges): plt.plot(e[:, 0], e[:, 1], "r--") plt.text(e[0, 0], e[0, 1], str(ei)) - + plt.plot(bh.edge_midpoints[:, 0], bh.edge_midpoints[:, 1], "ko") for e, n in zip(bh.edge_midpoints, normalized_normals): - p0 = e + p0 = e p = e + n*6 plt.plot([p0[0], p[0]], [p0[1], p[1]], "b--", lw=2) - - plt.scatter(vals[:, 0], vals[:, 1]) + + plt.scatter(vals[:, 0], vals[:, 1]) plt.imshow(mask, extent=[np.min(vals[:, 0]), np.max(vals[:, 0]), np.max(vals[:, 1]), np.min(vals[:, 1])]) plt.grid(True) @@ -181,7 +179,7 @@ def test_hull_construction(): plt.plot(e[:, 0], e[:, 1], "r--") plt.imshow(extracted_data[0, 0, :, :], extent=[extracted_window_extents[0], extracted_window_extents[1], - extracted_window_extents[3], extracted_window_extents[2]]) + extracted_window_extents[3], extracted_window_extents[2]]) plt.savefig("/tmp/extract_local.png") plt.figure(figsize=(7, 2.5)) @@ -192,26 +190,26 @@ def test_hull_construction(): for f in facet_regions: for ei, e in enumerate(f.edges): plt.plot(e[:, 0], e[:, 1], "co--") - - plt.imshow(stitched_image[0, 0, :, :], + + plt.imshow(stitched_image[0, 0, :, :], extent=[np.min(stitched_region.corners[:, 0]), np.max(stitched_region.corners[:, 0]), np.max(stitched_region.corners[:, 1]), np.min(stitched_region.corners[:, 1])]) plt.savefig("/tmp/facet.png") - + plt.figure(figsize=(7, 2.5)) plt.title("Overlapping faceting check") for f in [olap_box1, olap_box2, olap_box3]: for ei, e in enumerate(f.edges): plt.plot(e[:, 0], e[:, 1], "co--") - - plt.imshow(olaps_stitched_image[0, 0, :, :], + + plt.imshow(olaps_stitched_image[0, 0, :, :], extent=[np.min(olaps_stitched_region.corners[:, 0]), np.max(olaps_stitched_region.corners[:, 0]), np.max(olaps_stitched_region.corners[:, 1]), np.min(olaps_stitched_region.corners[:, 1])]) - plt.xlim((np.min(olaps_stitched_region.corners[:, 0]) - 15, + plt.xlim((np.min(olaps_stitched_region.corners[:, 0]) - 15, np.max(olaps_stitched_region.corners[:, 0]) + 15)) - plt.ylim((np.min(olaps_stitched_region.corners[:, 1]) - 15, + plt.ylim((np.min(olaps_stitched_region.corners[:, 1]) - 15, np.max(olaps_stitched_region.corners[:, 1]) + 15)) plt.savefig("/tmp/overlap_facet.png") From 162edd4a1d6130a9fa8f25dd984f77b74badf40a Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:30:41 +0200 Subject: [PATCH 05/10] Apply black formatting --- africanus/linalg/geometry.py | 441 ++++++++++++++++++------- africanus/linalg/test/test_geometry.py | 204 ++++++++---- 2 files changed, 461 insertions(+), 184 deletions(-) diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py index 7769c7b3d..1546becb3 100644 --- a/africanus/linalg/geometry.py +++ b/africanus/linalg/geometry.py @@ -19,19 +19,22 @@ def timed(*args, **kw): ts = time.time() result = method(*args, **kw) te = time.time() - if 'log_time' in kw: - name = kw.get('log_name', method.__name__.upper()) - kw['log_time'][name] = int((te - ts) * 1000) + if "log_time" in kw: + name = kw.get("log_name", method.__name__.upper()) + kw["log_time"][name] = int((te - ts) * 1000) else: - print(('%r %2.2f ms' % \ - (method.__name__, (te - ts) * 1000))) + print(("%r %2.2f ms" % (method.__name__, (te - ts) * 1000))) return result + return timed + @requires_optional("scipy.stats", opt_import_err) class BoundingConvexHull(object): - def __init__(self, list_hulls, name="unnamed", mask = None, check_mask_outofbounds=True): + def __init__( + self, list_hulls, name="unnamed", mask=None, check_mask_outofbounds=True + ): """ Initializes a bounding convex hull around a list of bounding convex hulls or series of points. A unity-weighted mask is computed for the region that falls within this convex hull if a mask of (y, x) coordinates is not provided. Otherwise if a mask is provided and the @@ -42,8 +45,12 @@ def __init__(self, list_hulls, name="unnamed", mask = None, check_mask_outofboun self._name = name self._check_mask_outofbounds = check_mask_outofbounds self._cached_filled_mask = None - self._vertices = points = np.vstack([b.corners - if hasattr(b, "corners") else [b[0], b[1]] for b in list_hulls]) + self._vertices = points = np.vstack( + [ + b.corners if hasattr(b, "corners") else [b[0], b[1]] + for b in list_hulls + ] + ) self._hull = spat.ConvexHull(points) if mask is None: self._mask, self._mask_weights = self.init_mask() @@ -56,22 +63,30 @@ def invalidate_cached_masks(self): self._mask, self._mask_weights = self.init_mask() def __str__(self): - return ",".join(["({0:d},{1:d})".format(x,y) for (x,y) in self.corners]) + return ",".join( + ["({0:d},{1:d})".format(x, y) for (x, y) in self.corners] + ) def init_mask(self): """ creates a sparse mask of the convex hull of the form (y, x) tuples """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) - x = np.arange(minx, maxx + 1, 1) #upper limit inclusive + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) + x = np.arange(minx, maxx + 1, 1) # upper limit inclusive y = np.arange(miny, maxy + 1, 1) meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) - sparse_mask = bounding_mesh if not self._check_mask_outofbounds else \ - [c for c in bounding_mesh if c[::-1] in self] + sparse_mask = ( + bounding_mesh + if not self._check_mask_outofbounds + else [c for c in bounding_mesh if c[::-1] in self] + ) - mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding + # initialize to unity, this should be modified when coadding + mask_weights = np.ones(len(sparse_mask)) return sparse_mask, mask_weights @property @@ -89,10 +104,18 @@ def sparse_mask(self, mask): """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ if not isinstance(mask, list): raise TypeError("Mask must be list") - if not (hasattr(mask, "__len__") and (len(mask) == 0 or (hasattr(mask[0], "__len__") and len(mask[0]) == 2))): + if not ( + hasattr(mask, "__len__") + and ( + len(mask) == 0 + or (hasattr(mask[0], "__len__") and len(mask[0]) == 2) + ) + ): raise TypeError("Mask must be a sparse mask of 2 element values") if self._check_mask_outofbounds: - self._mask = copy.deepcopy([c for c in mask if (c[1], c[0]) in self]) + self._mask = copy.deepcopy( + [c for c in mask if (c[1], c[0]) in self] + ) else: self._mask = copy.deepcopy(mask) self._mask_weights = np.ones(len(self._mask)) @@ -104,20 +127,28 @@ def mask(self, dtype=np.float64): return self._cached_filled_mask lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) - nx = maxx - minx + 1 # inclusive + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) + nx = maxx - minx + 1 # inclusive ny = maxy - miny + 1 - mesh = np.zeros(nx*ny, dtype=dtype) - if nx==0 or ny==0 or len(self.sparse_mask) == 0: + mesh = np.zeros(nx * ny, dtype=dtype) + if nx == 0 or ny == 0 or len(self.sparse_mask) == 0: self._cached_filled_mask = mesh.reshape((ny, nx)) else: sparse_mask = np.array(self.sparse_mask) - sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= minx, - sparse_mask[:, 1] <= maxx), - np.logical_and(sparse_mask[:, 0] >= miny, - sparse_mask[:, 0] <= maxy)) - flat_index = (sparse_mask[sel][:, 0] - miny)*nx + (sparse_mask[sel][:, 1] - minx) + sel = np.logical_and( + np.logical_and( + sparse_mask[:, 1] >= minx, sparse_mask[:, 1] <= maxx + ), + np.logical_and( + sparse_mask[:, 0] >= miny, sparse_mask[:, 0] <= maxy + ), + ) + flat_index = (sparse_mask[sel][:, 0] - miny) * nx + ( + sparse_mask[sel][:, 1] - minx + ) mesh[flat_index] = self._mask_weights[sel] self._cached_filled_mask = mesh.reshape((ny, nx)) return self._cached_filled_mask @@ -131,61 +162,104 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): assumes the last value of axes is the fastest varying axis """ if not isinstance(sel_region, BoundingConvexHull): - raise TypeError("Object passed in is not of type BoundingConvexHull") + raise TypeError( + "Object passed in is not of type BoundingConvexHull" + ) if not (hasattr(axes, "__len__") and len(axes) == 2): - raise ValueError("Expected a tupple of axes along which to slice out a region") + raise ValueError( + "Expected a tupple of axes along which to slice out a region" + ) axes = sorted(axes) - lines = np.hstack([sel_region.corners, np.roll(sel_region.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + lines = np.hstack( + [sel_region.corners, np.roll(sel_region.corners, -1, axis=0)] + ) + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) x = np.arange(minx, maxx + 1, 1) y = np.arange(miny, maxy + 1, 1) pad_left = max(0, 0 - minx) pad_bottom = max(0, 0 - miny) - pad_right = max(0, maxx - data_cube.shape[axes[1]] + 1) #inclusive of upper limit + # inclusive of upper limit + pad_right = max(0, maxx - data_cube.shape[axes[1]] + 1) pad_top = max(0, maxy - data_cube.shape[axes[0]] + 1) - if minx > data_cube.shape[axes[0]] or miny > data_cube.shape[axes[1]] or \ - maxy < 0 or maxx < 0: - raise ValueError("Expected a bounding hull that is at least partially within the image") + if ( + minx > data_cube.shape[axes[0]] + or miny > data_cube.shape[axes[1]] + or maxy < 0 + or maxx < 0 + ): + raise ValueError( + "Expected a bounding hull that is at least partially within the image" + ) # extract data, pad if necessary slc_data = [slice(None)] * len(data_cube.shape) - for (start, end), axis in zip([(miny + pad_bottom, maxy - pad_top + 1), - (minx + pad_left, maxx - pad_right + 1)], axes): + for (start, end), axis in zip( + [ + (miny + pad_bottom, maxy - pad_top + 1), + (minx + pad_left, maxx - pad_right + 1), + ], + axes, + ): slc_data[axis] = slice(start, end) slc_padded = [slice(None)] * len(data_cube.shape) - for (start, end), axis in zip([(pad_bottom, -miny + maxy + 1 - pad_top), - (pad_left, -minx + maxx + 1 - pad_right)], axes): + for (start, end), axis in zip( + [ + (pad_bottom, -miny + maxy + 1 - pad_top), + (pad_left, -minx + maxx + 1 - pad_right), + ], + axes, + ): slc_padded[axis] = slice(start, end) selected_data = data_cube[tuple(slc_data)] new_shape = list(data_cube.shape) - new_shape[axes[0]] = (maxy - miny + 1) - new_shape[axes[1]] = (maxx - minx + 1) + new_shape[axes[0]] = maxy - miny + 1 + new_shape[axes[1]] = maxx - minx + 1 if any(np.array([pad_left, pad_bottom, pad_right, pad_top]) > 0): - padded_data = np.zeros(tuple(new_shape), dtype=selected_data.dtype) * oob_value + padded_data = ( + np.zeros(tuple(new_shape), dtype=selected_data.dtype) + * oob_value + ) padded_data[tuple(slc_padded)] = selected_data.copy() else: padded_data = selected_data.copy() # finally apply mask slc_padded_data = [slice(None)] * len(padded_data.shape) - for (start, end), axis in zip([(0, maxy - miny + 1), #mask starts at origin in the padded image - (0, maxx - minx + 1)], axes): + for (start, end), axis in zip( + [ + ( + 0, + maxy - miny + 1, + ), # mask starts at origin in the padded image + (0, maxx - minx + 1), + ], + axes, + ): slc_padded_data[axis] = slice(start, end) slc_mask = [None] * len(padded_data.shape) - for (start, end), axis in zip([(0, sel_region.mask.shape[0]), #mask starts at origin in the padded image - (0, sel_region.mask.shape[1])], axes): + for (start, end), axis in zip( + [ + ( + 0, + sel_region.mask.shape[0], + ), # mask starts at origin in the padded image + (0, sel_region.mask.shape[1]), + ], + axes, + ): slc_mask[axis] = slice(start, end) mask = sel_region.mask.copy() mask[mask == 0] = oob_value padded_data[tuple(slc_padded_data)] *= mask[tuple(slc_mask)] - window_extents = [minx, maxx, - miny, maxy] + window_extents = [minx, maxx, miny, maxy] return padded_data, window_extents @classmethod @@ -203,7 +277,11 @@ def normalize_masks(cls, regions, only_overlapped_regions=True): # For now all regions have equal contribution allmasks = [] for reg in regions: - allmasks += list(reg.sparse_mask) if isinstance(reg.sparse_mask, np.ndarray) else reg.sparse_mask + allmasks += ( + list(reg.sparse_mask) + if isinstance(reg.sparse_mask, np.ndarray) + else reg.sparse_mask + ) # flatten for faster comparisons allmasks = np.array(allmasks) @@ -212,7 +290,9 @@ def normalize_masks(cls, regions, only_overlapped_regions=True): allmasks_flatten = allmasks[:, 0] * nx + allmasks[:, 1] # now count the number of times a pixel is painted onto - unique_pxls_flatten, paint_count = np.unique(allmasks_flatten, return_counts=True) + unique_pxls_flatten, paint_count = np.unique( + allmasks_flatten, return_counts=True + ) paint_count = paint_count.astype(np.float) if only_overlapped_regions: @@ -221,30 +301,48 @@ def normalize_masks(cls, regions, only_overlapped_regions=True): paint_count = paint_count[sel] # with the reduced number of overlap pixels unflatten - unique_pxls = np.vstack([unique_pxls_flatten // nx, - unique_pxls_flatten % nx]).T + unique_pxls = np.vstack( + [unique_pxls_flatten // nx, unique_pxls_flatten % nx] + ).T unique_pxls = list(map(tuple, unique_pxls)) paint_count[...] = 1.0 / paint_count # and finally update mask weights for reg in regions: - reg._cached_filled_mask = None # invalidate - overlap = [x for x in zip(paint_count, unique_pxls) if x[1] in reg.sparse_mask] + reg._cached_filled_mask = None # invalidate + overlap = [ + x + for x in zip(paint_count, unique_pxls) + if x[1] in reg.sparse_mask + ] for px_pc, px in overlap: - sel = reg.sparse_mask.index(px) if isinstance(reg.sparse_mask, list) else \ - np.all(reg.sparse_mask - px == 0, axis=1) + sel = ( + reg.sparse_mask.index(px) + if isinstance(reg.sparse_mask, list) + else np.all(reg.sparse_mask - px == 0, axis=1) + ) reg._mask_weights[sel] = px_pc @property def circumference(self): """ area contained in hull """ lines = self.edges - return np.sum(np.linalg.norm(lines[:, 1, :] - lines[:, 0, :], axis=1) + 1) + return np.sum( + np.linalg.norm(lines[:, 1, :] - lines[:, 0, :], axis=1) + 1 + ) + @property def area(self): """ area contained in hull """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - return 0.5 * np.abs(np.sum([x1*(y2)-(x2)*y1 for x1,y1,x2,y2 in lines])) + 0.5 * self.circumference - 1 + return ( + 0.5 + * np.abs( + np.sum([x1 * (y2) - (x2) * y1 for x1, y1, x2, y2 in lines]) + ) + + 0.5 * self.circumference + - 1 + ) @property def name(self): @@ -259,7 +357,7 @@ def corners(self): """ Returns vertices and guarentees clockwise winding """ return self._vertices[self._hull.vertices][::-1] - def normals(self, left = True): + def normals(self, left=True): """ return a list of left normals to the hull """ normals = [] for i in range(self.corners.shape[0]): @@ -292,14 +390,16 @@ def edge_midpoints(self): @property def lnormals(self): """ left normals to the edges of the hull """ - return self.normals(left = True) + return self.normals(left=True) @property def rnormals(self): """ right normals to the edges of the hull """ return self.normals(left=False) - def overlaps_with(self, other, min_sep_dist=0.5): #less than half a pixel away + def overlaps_with( + self, other, min_sep_dist=0.5 + ): # less than half a pixel away """ Implements the separating lines collision detection theorem to test whether the hull intersects with 'other' hull @@ -324,7 +424,10 @@ def overlaps_with(self, other, min_sep_dist=0.5): #less than half a pixel away projs = np.dot(vecs_reg2, n.T) minproj_reg2 = np.min(projs) maxproj_reg2 = np.max(projs) - if minproj_reg2 - maxproj_reg1 > min_sep_dist or minproj_reg1 - maxproj_reg2 > min_sep_dist: + if ( + minproj_reg2 - maxproj_reg1 > min_sep_dist + or minproj_reg1 - maxproj_reg2 > min_sep_dist + ): return False return True @@ -332,13 +435,15 @@ def overlaps_with(self, other, min_sep_dist=0.5): #less than half a pixel away def centre(self, integral=True): """ Barycentre of hull """ if integral: + def rnd(x): return int(np.floor(x) if x >= 0 else np.ceil(x)) + return [rnd(x) for x in np.mean(self._vertices, axis=0)] else: return np.mean(self._vertices, axis=0) - def __contains__(self, s, tolerance=0.5): #less than half a pixel away + def __contains__(self, s, tolerance=0.5): # less than half a pixel away """ tests whether a point s(x,y) is in the convex hull """ # there are three cases to consider # CASE 1: @@ -358,38 +463,52 @@ def __contains__(self, s, tolerance=0.5): #less than half a pixel away dot = np.einsum("ij,ij->i", normals, xyvec) return np.all(dot > -tolerance) + class BoundingBox(BoundingConvexHull): def __init__(self, xl, xu, yl, yu, name="unnamed", mask=None, **kwargs): - if not all([isinstance(x, (int, np.int64, np.int32, np.int16)) for x in [xl, xu, yl, yu]]): + if not all( + [ + isinstance(x, (int, np.int64, np.int32, np.int16)) + for x in [xl, xu, yl, yu] + ] + ): raise ValueError("Box limits must be integers") - self.__xnpx = abs(xu - xl + 1) #inclusive of the upper pixel + self.__xnpx = abs(xu - xl + 1) # inclusive of the upper pixel self.__ynpx = abs(yu - yl + 1) - BoundingConvexHull.__init__(self, - [[xl,yl],[xl,yu],[xu,yu],[xu,yl]], - name, - mask=mask, - **kwargs) + BoundingConvexHull.__init__( + self, + [[xl, yl], [xl, yu], [xu, yu], [xu, yl]], + name, + mask=mask, + **kwargs + ) def init_mask(self): """ creates a sparse mask of the convex hull of the form (y, x) tuples """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) - x = np.arange(minx, maxx + 1, 1) #upper limit inclusive + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) + x = np.arange(minx, maxx + 1, 1) # upper limit inclusive y = np.arange(miny, maxy + 1, 1) meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) - sparse_mask = np.asarray(bounding_mesh) # by default for a BB region the mask is always going to be the entire region + # by default for a BB region the mask is always going to be the entire region + sparse_mask = np.asarray(bounding_mesh) - mask_weights = np.ones(len(sparse_mask)) #initialize to unity, this should be modified when coadding + # initialize to unity, this should be modified when coadding + mask_weights = np.ones(len(sparse_mask)) return sparse_mask, mask_weights def __contains__(self, s): """ tests whether a point s(x,y) is in the box""" lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) return s[0] >= minx and s[0] <= maxx and s[1] >= miny and s[1] <= maxy @property @@ -406,37 +525,68 @@ def sparse_mask(self, mask): """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ if not isinstance(mask, list) and not isinstance(mask, np.ndarray): raise TypeError("Mask must be list") - if not (hasattr(mask, "__len__") and (len(mask) == 0 or (hasattr(mask[0], "__len__") and len(mask[0]) == 2))): + if not ( + hasattr(mask, "__len__") + and ( + len(mask) == 0 + or (hasattr(mask[0], "__len__") and len(mask[0]) == 2) + ) + ): raise TypeError("Mask must be a sparse mask of 2 element values") if mask == []: self._mask = [] else: lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) - nx = maxx - minx + 1 # inclusive + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) + nx = maxx - minx + 1 # inclusive ny = maxy - miny + 1 sparse_mask = np.asarray(mask) - sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= minx, - sparse_mask[:, 1] <= maxx), - np.logical_and(sparse_mask[:, 0] >= miny, - sparse_mask[:, 0] <= maxy)) + sel = np.logical_and( + np.logical_and( + sparse_mask[:, 1] >= minx, sparse_mask[:, 1] <= maxx + ), + np.logical_and( + sparse_mask[:, 0] >= miny, sparse_mask[:, 0] <= maxy + ), + ) self._mask = sparse_mask[sel] self._mask_weights = np.ones(len(self._mask)) @classmethod - def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np.float64, **kwargs): + def project_regions( + cls, + regional_data_list, + regions_list, + axes=(2, 3), + dtype=np.float64, + **kwargs + ): """ Projects individial regions back onto a single contiguous cube """ - if not (hasattr(regional_data_list, "__len__") and hasattr(regions_list, "__len__") and \ - len(regions_list) == len(regional_data_list)): - raise TypeError("Region data list and regions lists must be lists of equal length") + if not ( + hasattr(regional_data_list, "__len__") + and hasattr(regions_list, "__len__") + and len(regions_list) == len(regional_data_list) + ): + raise TypeError( + "Region data list and regions lists must be lists of equal length" + ) if not all([isinstance(x, np.ndarray) for x in regional_data_list]): raise TypeError("Region data list must be a list of ndarrays") if not all([isinstance(x, BoundingBox) for x in regions_list]): - raise TypeError("Region list must be a list of Axis Aligned Bounding Boxes") + raise TypeError( + "Region list must be a list of Axis Aligned Bounding Boxes" + ) if regions_list == []: return np.empty((0)) - if not all([reg.ndim == regional_data_list[0].ndim for reg in regional_data_list]): + if not all( + [ + reg.ndim == regional_data_list[0].ndim + for reg in regional_data_list + ] + ): raise ValueError("All data cubes must be of equal dimension") axes = tuple(sorted(axes)) @@ -446,8 +596,8 @@ def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np maxy = np.max([np.max(f.corners[:, 1]) for f in regions_list]) npxx = maxx - minx + 1 npxy = maxy - miny + 1 - global_offsetx = -minx #-min(0, minx) - global_offsety = -miny #-min(0, miny) + global_offsetx = -minx # -min(0, minx) + global_offsety = -miny # -min(0, miny) projected_image_size = list(regional_data_list[0].shape) projected_image_size[axes[0]] = npxy @@ -457,14 +607,16 @@ def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np combined_mask = [] for f, freg in zip(regional_data_list, regions_list): f[np.isnan(f)] = 0 - xl = max(0, global_offsetx+np.min(freg.corners[:, 0])) - xu = min(global_offsetx+np.max(freg.corners[:, 0]) + 1, npxx) - yl = max(0, global_offsety+np.min(freg.corners[:, 1])) - yu = min(global_offsety+np.max(freg.corners[:, 1]) + 1, npxy) - fnx = xu - xl + 1 # inclusive - fny = yu - yl + 1 # inclusive + xl = max(0, global_offsetx + np.min(freg.corners[:, 0])) + xu = min(global_offsetx + np.max(freg.corners[:, 0]) + 1, npxx) + yl = max(0, global_offsety + np.min(freg.corners[:, 1])) + yu = min(global_offsety + np.max(freg.corners[:, 1]) + 1, npxy) + fnx = xu - xl + 1 # inclusive + fny = yu - yl + 1 # inclusive if f.shape[axes[0]] != fny - 1 or f.shape[axes[1]] != fnx - 1: - raise ValueError("One or more bounding box descriptors does not match shape of corresponding data cubes") + raise ValueError( + "One or more bounding box descriptors does not match shape of corresponding data cubes" + ) slc_data = [slice(None)] * len(stitched_img.shape) for (start, end), axis in zip([(yl, yu), (xl, xu)], axes): slc_data[axis] = slice(start, end) @@ -472,21 +624,38 @@ def project_regions(cls, regional_data_list, regions_list, axes=(2, 3), dtype=np stitched_img[tuple(slc_data)] += f combined_mask += list(freg.sparse_mask) - return stitched_img, BoundingBox(minx, maxx, miny, maxy, mask=combined_mask, **kwargs) + return ( + stitched_img, + BoundingBox(minx, maxx, miny, maxy, mask=combined_mask, **kwargs), + ) + ######################################################################## -## Factories +# Factories ######################################################################## + class BoundingBoxFactory(object): @classmethod - def AxisAlignedBoundingBox(cls, convex_hull_object, square=False, enforce_odd=True, **kwargs): + def AxisAlignedBoundingBox( + cls, convex_hull_object, square=False, enforce_odd=True, **kwargs + ): """ Constructs an axis aligned bounding box around convex hull """ if not isinstance(convex_hull_object, BoundingConvexHull): - raise TypeError("Convex hull object passed in constructor is not of type BoundingConvexHull") + raise TypeError( + "Convex hull object passed in constructor is not of type BoundingConvexHull" + ) if square: - nx = np.max(convex_hull_object.corners[:, 0]) - np.min(convex_hull_object.corners[:, 0]) + 1 #inclusive - ny = np.max(convex_hull_object.corners[:, 1]) - np.min(convex_hull_object.corners[:, 1]) + 1 #inclusive + nx = ( + np.max(convex_hull_object.corners[:, 0]) + - np.min(convex_hull_object.corners[:, 0]) + + 1 + ) # inclusive + ny = ( + np.max(convex_hull_object.corners[:, 1]) + - np.min(convex_hull_object.corners[:, 1]) + + 1 + ) # inclusive boxdiam = max(nx, ny) boxrad = boxdiam // 2 cx, cy = convex_hull_object.centre @@ -503,10 +672,15 @@ def AxisAlignedBoundingBox(cls, convex_hull_object, square=False, enforce_odd=Tr xu += (xu - xl) % 2 if enforce_odd else 0 yu += (yu - yl) % 2 if enforce_odd else 0 - return BoundingBox(xl, xu, yl, yu, - convex_hull_object.name, - mask=convex_hull_object.sparse_mask, - **kwargs) + return BoundingBox( + xl, + xu, + yl, + yu, + convex_hull_object.name, + mask=convex_hull_object.sparse_mask, + **kwargs + ) @classmethod def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): @@ -521,8 +695,12 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): yu = np.max(bounding_box_object.corners[:, 1]) # construct a nonregular meshgrid bound to xu and yu - x = xl + np.arange(0, nsubboxes + 1) * int(np.ceil((xu - xl + 1) / float(nsubboxes))) - y = yl + np.arange(0, nsubboxes + 1) * int(np.ceil((yu - yl + 1) / float(nsubboxes))) + x = xl + np.arange(0, nsubboxes + 1) * int( + np.ceil((xu - xl + 1) / float(nsubboxes)) + ) + y = yl + np.arange(0, nsubboxes + 1) * int( + np.ceil((yu - yl + 1) / float(nsubboxes)) + ) xx, yy = np.meshgrid(x, y) # split into boxes @@ -539,7 +717,7 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): xus[:, -1] = max(xu, min(xus[0, -1], xu)) yus[-1, :] = max(yu, min(yus[-1, 0], yu)) - #coordinates for all the contained boxes, anti-clockwise wound + # coordinates for all the contained boxes, anti-clockwise wound xls = xls.ravel() yls = yls.ravel() xus = xus.ravel() @@ -550,13 +728,20 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): ul = list(zip(xls, yus)) contained_boxes = list(zip(bl, br, ur, ul)) - #finally create bbs for each of the contained boxes with the mask - #chopped up between the boxes by the convex hull initializer - new_regions = [BoundingBox(bl[0], br[0], bl[1], ul[1], - bounding_box_object.name, - mask=bounding_box_object.sparse_mask, - **kwargs) - for bl, br, ur, ul in contained_boxes] + # finally create bbs for each of the contained boxes with the mask + # chopped up between the boxes by the convex hull initializer + new_regions = [ + BoundingBox( + bl[0], + br[0], + bl[1], + ul[1], + bounding_box_object.name, + mask=bounding_box_object.sparse_mask, + **kwargs + ) + for bl, br, ur, ul in contained_boxes + ] return new_regions @@ -577,8 +762,12 @@ def PadBox(cls, bounding_box_object, desired_nx, desired_ny, **kwargs): xu = cx + pad_right yl = cy - pad_bottom yu = cy + pad_top - return BoundingBox(xl, xu, yl, yu, - bounding_box_object.name, - mask=bounding_box_object.sparse_mask, - **kwargs) #mask unchanged in the new shape, border frame discarded - + return BoundingBox( + xl, + xu, + yl, + yu, + bounding_box_object.name, + mask=bounding_box_object.sparse_mask, + **kwargs + ) # mask unchanged in the new shape, border frame discarded diff --git a/africanus/linalg/test/test_geometry.py b/africanus/linalg/test/test_geometry.py index f4f90e7d3..48809c2f8 100644 --- a/africanus/linalg/test/test_geometry.py +++ b/africanus/linalg/test/test_geometry.py @@ -1,49 +1,69 @@ import numpy as np import pytest -from africanus.linalg.geometry import (BoundingConvexHull, - BoundingBox, - BoundingBoxFactory) +from africanus.linalg.geometry import ( + BoundingConvexHull, + BoundingBox, + BoundingBoxFactory, +) + + @pytest.mark.parametrize("debug", [False]) def test_hull_construction(debug): # test case 1 vals = np.array([[50, 60], [20, 40], [-74, 50], [-95, +10], [20, 60]]) bh = BoundingConvexHull(vals) mask = bh.mask - assert mask.shape == (np.max(vals[:, 1]) - np.min(vals[:, 1]) + 1, np.max(vals[:, 0]) - np.min(vals[:, 0]) + 1) - assert np.abs(mask.sum() - bh.area) / bh.area < 0.05 # integral mask area needs to be close to true area - normalized_normals = bh.rnormals / np.linalg.norm(bh.rnormals, axis=1)[:, None] + assert mask.shape == ( + np.max(vals[:, 1]) - np.min(vals[:, 1]) + 1, + np.max(vals[:, 0]) - np.min(vals[:, 0]) + 1, + ) + # integral mask area needs to be close to true area + assert np.abs(mask.sum() - bh.area) / bh.area < 0.05 + normalized_normals = ( + bh.rnormals / np.linalg.norm(bh.rnormals, axis=1)[:, None] + ) # test case 2 for e, n in zip(bh.edges, normalized_normals): edge_vec = e[1] - e[0] assert np.all(np.abs(np.dot(edge_vec, n)) < 1.0e-8) # test case 3 - valsextract = np.array([[-10, 120], [90, 268], [293, 110],[40, -30]]) + valsextract = np.array([[-10, 120], [90, 268], [293, 110], [40, -30]]) bh_extract = BoundingConvexHull(valsextract) sinc_npx = 255 sinc = np.sinc(np.linspace(-7, 7, sinc_npx)) sinc2d = np.outer(sinc, sinc).reshape((1, 1, sinc_npx, sinc_npx)) - extracted_data, extracted_window_extents = BoundingConvexHull.regional_data(bh_extract, sinc2d, oob_value=np.nan) + extracted_data, extracted_window_extents = BoundingConvexHull.regional_data( + bh_extract, sinc2d, oob_value=np.nan + ) assert extracted_window_extents == [-10, 293, -30, 268] sparse_mask = np.array(bh_extract.sparse_mask) - lines = np.hstack([bh_extract.corners, np.roll(bh_extract.corners, -1, axis=0)]) - minx = np.min(lines[:, 0:4:2]); maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]); maxy = np.max(lines[:, 1:4:2]) - sel = np.logical_and(np.logical_and(sparse_mask[:, 1] >= 0, - sparse_mask[:, 1] < 255), - np.logical_and(sparse_mask[:, 0] >= 0, - sparse_mask[:, 0] < 255)) - - flat_index = (sparse_mask[sel][:, 0])*sinc_npx + (sparse_mask[sel][:, 1]) + lines = np.hstack( + [bh_extract.corners, np.roll(bh_extract.corners, -1, axis=0)] + ) + minx = np.min(lines[:, 0:4:2]) + maxx = np.max(lines[:, 0:4:2]) + miny = np.min(lines[:, 1:4:2]) + maxy = np.max(lines[:, 1:4:2]) + sel = np.logical_and( + np.logical_and(sparse_mask[:, 1] >= 0, sparse_mask[:, 1] < 255), + np.logical_and(sparse_mask[:, 0] >= 0, sparse_mask[:, 0] < 255), + ) + + flat_index = (sparse_mask[sel][:, 0]) * sinc_npx + (sparse_mask[sel][:, 1]) sinc_integral = np.sum(sinc2d.ravel()[flat_index]) assert np.abs(sinc_integral - np.nansum(extracted_data.ravel())) < 1.0e-8 v = np.nanargmax(extracted_data) - vx = v % extracted_data.shape[3]; vy = v // extracted_data.shape[3] - cextracted = (extracted_window_extents[0] + vx, - extracted_window_extents[2] + vy) + vx = v % extracted_data.shape[3] + vy = v // extracted_data.shape[3] + cextracted = ( + extracted_window_extents[0] + vx, + extracted_window_extents[2] + vy, + ) v = np.nanargmax(sinc2d) - sincvx = v % sinc_npx; sincvy = v // sinc_npx + sincvx = v % sinc_npx + sincvy = v // sinc_npx csinc = tuple([sincvx, sincvy]) assert csinc == cextracted @@ -73,17 +93,19 @@ def test_hull_construction(debug): assert (-15, 35) not in bb assert (0, 35) in bb - bb2 = BoundingBoxFactory.AxisAlignedBoundingBox(bb) #enforce odd + bb2 = BoundingBoxFactory.AxisAlignedBoundingBox(bb) # enforce odd assert bb2.box_npx == (35, 21) assert bb2.area == 35 * 21 assert (bb.sparse_mask == bb2.sparse_mask).all() assert (-15, 35) not in bb2 assert (0, 35) in bb2 - bb3 = BoundingBoxFactory.AxisAlignedBoundingBox(bb, square=True) #enforce odd + bb3 = BoundingBoxFactory.AxisAlignedBoundingBox( + bb, square=True + ) # enforce odd assert bb3.box_npx[0] == bb3.box_npx[1] - assert bb3.box_npx[0] % 2 == 1 #enforce odd - assert bb3.area == bb3.box_npx[0]**2 + assert bb3.box_npx[0] % 2 == 1 # enforce odd + assert bb3.area == bb3.box_npx[0] ** 2 assert (bb.sparse_mask == bb3.sparse_mask).all() assert (-15, 35) not in bb2 assert (0, 35) in bb2 @@ -91,8 +113,12 @@ def test_hull_construction(debug): # test case 7 bb4s = BoundingBoxFactory.SplitBox(bb, nsubboxes=3) assert len(bb4s) == 9 - xlims = [(np.min(c.corners[:, 0]), np.max(c.corners[:, 0])) for c in bb4s][0:3] - ylims = [(np.min(c.corners[:, 1]), np.max(c.corners[:, 1])) for c in bb4s][0::3] + xlims = [(np.min(c.corners[:, 0]), np.max(c.corners[:, 0])) for c in bb4s][ + 0:3 + ] + ylims = [(np.min(c.corners[:, 1]), np.max(c.corners[:, 1])) for c in bb4s][ + 0::3 + ] assert np.all(xlims == np.array([(-14, -3), (-2, 9), (10, 20)])) assert np.all(ylims == np.array([(30, 36), (37, 43), (44, 49)])) assert np.sum([b.area for b in bb4s]) == bb.area @@ -109,19 +135,43 @@ def test_hull_construction(debug): assert np.sum(bb5.mask) == np.sum(bb6.mask) bb7s = list(map(lambda x: BoundingBoxFactory.PadBox(x, 17, 11), bb4s)) assert all([b.box_npx == (17, 11) for b in bb7s]) - assert np.sum([np.sum(b.mask) for b in bb7s]) == np.sum([np.sum(b.mask) for b in bb4s]) + assert np.sum([np.sum(b.mask) for b in bb7s]) == np.sum( + [np.sum(b.mask) for b in bb4s] + ) # test case 9 - facet_regions = list(map(lambda f: BoundingBoxFactory.PadBox(f, 63, 63), - BoundingBoxFactory.SplitBox(BoundingBoxFactory.AxisAlignedBoundingBox(bh_extract), nsubboxes=5))) - facets = list(map(lambda pf: BoundingConvexHull.regional_data(pf, sinc2d, oob_value=np.nan), - facet_regions)) - stitched_image, stitched_region = BoundingBox.project_regions([f[0] for f in facets], facet_regions) - assert np.abs(sinc_integral - np.nansum([np.nansum(f[0]) for f in facets])) < 1.0e-8 + facet_regions = list( + map( + lambda f: BoundingBoxFactory.PadBox(f, 63, 63), + BoundingBoxFactory.SplitBox( + BoundingBoxFactory.AxisAlignedBoundingBox(bh_extract), + nsubboxes=5, + ), + ) + ) + facets = list( + map( + lambda pf: BoundingConvexHull.regional_data( + pf, sinc2d, oob_value=np.nan + ), + facet_regions, + ) + ) + stitched_image, stitched_region = BoundingBox.project_regions( + [f[0] for f in facets], facet_regions + ) + assert ( + np.abs(sinc_integral - np.nansum([np.nansum(f[0]) for f in facets])) + < 1.0e-8 + ) assert np.abs(sinc_integral - np.sum(stitched_image)) < 1.0e-8 v = np.argmax(stitched_image) - vx = v % stitched_image.shape[3]; vy = v // stitched_image.shape[3] - cstitched = (np.min(stitched_region.corners[:, 0]) + vx, np.min(stitched_region.corners[:, 1]) + vy) + vx = v % stitched_image.shape[3] + vy = v // stitched_image.shape[3] + cstitched = ( + np.min(stitched_region.corners[:, 0]) + vx, + np.min(stitched_region.corners[:, 1]) + vy, + ) assert cstitched == csinc # test case 10 @@ -132,18 +182,23 @@ def test_hull_construction(debug): ext1 = BoundingConvexHull.regional_data(olap_box1, sinc2d)[0] ext2 = BoundingConvexHull.regional_data(olap_box2, sinc2d)[0] ext3 = BoundingConvexHull.regional_data(olap_box3, sinc2d)[0] - olaps_stitched_image, olaps_stitched_region = BoundingBox.project_regions([ext1, ext2, ext3], - [olap_box1, olap_box2, olap_box3]) + olaps_stitched_image, olaps_stitched_region = BoundingBox.project_regions( + [ext1, ext2, ext3], [olap_box1, olap_box2, olap_box3] + ) v = np.nanargmax(olaps_stitched_image) - vx = v % olaps_stitched_image.shape[3]; vy = v // olaps_stitched_image.shape[3] - cstitched_olap = (np.min(olaps_stitched_region.corners[:, 0]) + vx, - np.min(olaps_stitched_region.corners[:, 1]) + vy) + vx = v % olaps_stitched_image.shape[3] + vy = v // olaps_stitched_image.shape[3] + cstitched_olap = ( + np.min(olaps_stitched_region.corners[:, 0]) + vx, + np.min(olaps_stitched_region.corners[:, 1]) + vy, + ) assert cstitched_olap == csinc assert np.abs(1.0 - np.nanmax(olaps_stitched_image)) < 1.0e-8 # visual inspection if debug: from matplotlib import pyplot as plt + plt.figure(figsize=(7, 2.5)) plt.title("Winding, normals and masking check") for h in [bh, bh2, bh3]: @@ -154,11 +209,19 @@ def test_hull_construction(debug): plt.plot(bh.edge_midpoints[:, 0], bh.edge_midpoints[:, 1], "ko") for e, n in zip(bh.edge_midpoints, normalized_normals): p0 = e - p = e + n*6 + p = e + n * 6 plt.plot([p0[0], p[0]], [p0[1], p[1]], "b--", lw=2) plt.scatter(vals[:, 0], vals[:, 1]) - plt.imshow(mask, extent=[np.min(vals[:, 0]), np.max(vals[:, 0]), np.max(vals[:, 1]), np.min(vals[:, 1])]) + plt.imshow( + mask, + extent=[ + np.min(vals[:, 0]), + np.max(vals[:, 0]), + np.max(vals[:, 1]), + np.min(vals[:, 1]), + ], + ) plt.grid(True) plt.savefig("/tmp/winding.png") @@ -177,9 +240,15 @@ def test_hull_construction(debug): for h in [bh_extract]: for ei, e in enumerate(h.edges): plt.plot(e[:, 0], e[:, 1], "r--") - plt.imshow(extracted_data[0, 0, :, :], - extent=[extracted_window_extents[0], extracted_window_extents[1], - extracted_window_extents[3], extracted_window_extents[2]]) + plt.imshow( + extracted_data[0, 0, :, :], + extent=[ + extracted_window_extents[0], + extracted_window_extents[1], + extracted_window_extents[3], + extracted_window_extents[2], + ], + ) plt.savefig("/tmp/extract_local.png") plt.figure(figsize=(7, 2.5)) @@ -191,10 +260,15 @@ def test_hull_construction(debug): for ei, e in enumerate(f.edges): plt.plot(e[:, 0], e[:, 1], "co--") - - plt.imshow(stitched_image[0, 0, :, :], - extent=[np.min(stitched_region.corners[:, 0]), np.max(stitched_region.corners[:, 0]), - np.max(stitched_region.corners[:, 1]), np.min(stitched_region.corners[:, 1])]) + plt.imshow( + stitched_image[0, 0, :, :], + extent=[ + np.min(stitched_region.corners[:, 0]), + np.max(stitched_region.corners[:, 0]), + np.max(stitched_region.corners[:, 1]), + np.min(stitched_region.corners[:, 1]), + ], + ) plt.savefig("/tmp/facet.png") plt.figure(figsize=(7, 2.5)) @@ -203,15 +277,29 @@ def test_hull_construction(debug): for ei, e in enumerate(f.edges): plt.plot(e[:, 0], e[:, 1], "co--") - - plt.imshow(olaps_stitched_image[0, 0, :, :], - extent=[np.min(olaps_stitched_region.corners[:, 0]), np.max(olaps_stitched_region.corners[:, 0]), - np.max(olaps_stitched_region.corners[:, 1]), np.min(olaps_stitched_region.corners[:, 1])]) - plt.xlim((np.min(olaps_stitched_region.corners[:, 0]) - 15, - np.max(olaps_stitched_region.corners[:, 0]) + 15)) - plt.ylim((np.min(olaps_stitched_region.corners[:, 1]) - 15, - np.max(olaps_stitched_region.corners[:, 1]) + 15)) + plt.imshow( + olaps_stitched_image[0, 0, :, :], + extent=[ + np.min(olaps_stitched_region.corners[:, 0]), + np.max(olaps_stitched_region.corners[:, 0]), + np.max(olaps_stitched_region.corners[:, 1]), + np.min(olaps_stitched_region.corners[:, 1]), + ], + ) + plt.xlim( + ( + np.min(olaps_stitched_region.corners[:, 0]) - 15, + np.max(olaps_stitched_region.corners[:, 0]) + 15, + ) + ) + plt.ylim( + ( + np.min(olaps_stitched_region.corners[:, 1]) - 15, + np.max(olaps_stitched_region.corners[:, 1]) + 15, + ) + ) plt.savefig("/tmp/overlap_facet.png") + if __name__ == "__main__": test_hull_construction() From d5646f6a126baeb0443e5a66035bfd0246e85f70 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:51:18 +0200 Subject: [PATCH 06/10] line length fixes, remove unused vars --- africanus/linalg/geometry.py | 120 ++++++++++++++----------- africanus/linalg/test/test_geometry.py | 10 +-- 2 files changed, 72 insertions(+), 58 deletions(-) diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py index 1546becb3..0d95a6b69 100644 --- a/africanus/linalg/geometry.py +++ b/africanus/linalg/geometry.py @@ -32,15 +32,21 @@ def timed(*args, **kw): @requires_optional("scipy.stats", opt_import_err) class BoundingConvexHull(object): - def __init__( - self, list_hulls, name="unnamed", mask=None, check_mask_outofbounds=True - ): - """ Initializes a bounding convex hull around a list of bounding convex hulls or series of - points. A unity-weighted mask is computed for the region that falls within this convex hull - if a mask of (y, x) coordinates is not provided. Otherwise if a mask is provided and the - check_mask_outofbounds value is set the masked coordinates are not verified to fall within - the hull. The latter should thus be used with some caution by the user, but can potentially - significantly speed up the mask creation process for axis aligned regions. + def __init__(self, list_hulls, name="unnamed", + mask=None, check_mask_outofbounds=True): + """ + Initializes a bounding convex hull around a list of bounding + convex hulls or series of points. + A unity-weighted mask is computed for the region that + falls within this convex hull + if a mask of (y, x) coordinates is not provided. + Otherwise if a mask is provided and the + check_mask_outofbounds value is set the + masked coordinates are not verified to fall within + the hull. The latter should thus be used with some + caution by the user, but can potentially + significantly speed up the mask creation process for + axis aligned regions. """ self._name = name self._check_mask_outofbounds = check_mask_outofbounds @@ -68,7 +74,9 @@ def __str__(self): ) def init_mask(self): - """ creates a sparse mask of the convex hull of the form (y, x) tuples """ + """ + Creates a sparse mask of the convex hull of the form (y, x) tuples + """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) minx = np.min(lines[:, 0:4:2]) maxx = np.max(lines[:, 0:4:2]) @@ -76,7 +84,6 @@ def init_mask(self): maxy = np.max(lines[:, 1:4:2]) x = np.arange(minx, maxx + 1, 1) # upper limit inclusive y = np.arange(miny, maxy + 1, 1) - meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) sparse_mask = ( @@ -96,12 +103,17 @@ def sprase_mask_weights(self): @property def sparse_mask(self): - """ returns a sparse mask (y, x) values of all points in the masked region """ + """ + Returns a sparse mask (y, x) values of all points in the masked region + """ return self._mask @sparse_mask.setter def sparse_mask(self, mask): - """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ + """ + Sets the mask of the hull from + a sparse mask - list of (y, x) coordinates + """ if not isinstance(mask, list): raise TypeError("Mask must be list") if not ( @@ -178,8 +190,6 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): maxx = np.max(lines[:, 0:4:2]) miny = np.min(lines[:, 1:4:2]) maxy = np.max(lines[:, 1:4:2]) - x = np.arange(minx, maxx + 1, 1) - y = np.arange(miny, maxy + 1, 1) pad_left = max(0, 0 - minx) pad_bottom = max(0, 0 - miny) @@ -193,9 +203,8 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): or maxy < 0 or maxx < 0 ): - raise ValueError( - "Expected a bounding hull that is at least partially within the image" - ) + raise ValueError("Expected a bounding hull that is " + "at least partially within the image") # extract data, pad if necessary slc_data = [slice(None)] * len(data_cube.shape) @@ -264,10 +273,12 @@ def regional_data(cls, sel_region, data_cube, axes=(2, 3), oob_value=0): @classmethod def normalize_masks(cls, regions, only_overlapped_regions=True): - """ Normalizes region masks for overlapping pixels. This is necessary to properly coadd - overlapping facets. If masks are guarenteed to be initialized to unity (e.g. after - bounding region creation) the user can skip normalizing non-overlapping regions with - flag only_overlapped_regions. + """ + Normalizes region masks for overlapping pixels. This is necessary to + properly coadd overlapping facets. + If masks are guarenteed to be initialized to unity (e.g. after + bounding region creation) the user can skip normalizing non-overlapping + regions with flag only_overlapped_regions. """ if not all([isinstance(reg, BoundingConvexHull) for reg in regions]): raise TypeError("Expected a list of bounding convex hulls") @@ -397,12 +408,11 @@ def rnormals(self): """ right normals to the edges of the hull """ return self.normals(left=False) - def overlaps_with( - self, other, min_sep_dist=0.5 - ): # less than half a pixel away + def overlaps_with(self, other, min_sep_dist=0.5): + # less than half a pixel away """ - Implements the separating lines collision detection theorem - to test whether the hull intersects with 'other' hull + Implements the separating lines collision detection theorem + to test whether the hull intersects with 'other' hull """ if not isinstance(other, BoundingConvexHull): raise TypeError("rhs must be a BoundingConvexHull") @@ -447,16 +457,17 @@ def __contains__(self, s, tolerance=0.5): # less than half a pixel away """ tests whether a point s(x,y) is in the convex hull """ # there are three cases to consider # CASE 1: - # scalar projection between all inner pointing right normals (clockwise winding) + # scalar projection between all + # inner pointing right normals (clockwise winding) # and the point must be positive if the point were to lie inside # the region (true) # CASE 2: - # point is on an edge - the scalar projection onto the axis is 0 for that edge + # point is on an edge - the scalar projection onto + # the axis is 0 for that edge # and greater than 0 for the other edges (true) # CASE 3: # it is outside (false) x, y = s - isin = True normals = self.rnormals xyvec = np.array([x, y])[None, :] - np.array(self.corners) @@ -484,7 +495,9 @@ def __init__(self, xl, xu, yl, yu, name="unnamed", mask=None, **kwargs): ) def init_mask(self): - """ creates a sparse mask of the convex hull of the form (y, x) tuples """ + """ + creates a sparse mask of the convex hull of the form (y, x) tuples + """ lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) minx = np.min(lines[:, 0:4:2]) maxx = np.max(lines[:, 0:4:2]) @@ -492,10 +505,10 @@ def init_mask(self): maxy = np.max(lines[:, 1:4:2]) x = np.arange(minx, maxx + 1, 1) # upper limit inclusive y = np.arange(miny, maxy + 1, 1) - meshgrid = np.meshgrid(y, x) bounding_mesh = list(zip(*[np.ravel(x) for x in np.meshgrid(y, x)])) - # by default for a BB region the mask is always going to be the entire region + # by default for a BB region the mask is + # always going to be the entire region sparse_mask = np.asarray(bounding_mesh) # initialize to unity, this should be modified when coadding @@ -517,12 +530,17 @@ def box_npx(self): @property def sparse_mask(self): - """ returns a sparse mask (y, x) values of all points in the masked region """ + """ + returns a sparse mask (y, x) values of all points in the masked region + """ return self._mask @sparse_mask.setter def sparse_mask(self, mask): - """ Sets the mask of the hull from a sparse mask - list of (y, x) coordinates """ + """ + Sets the mask of the hull from a + sparse mask - list of (y, x) coordinates + """ if not isinstance(mask, list) and not isinstance(mask, np.ndarray): raise TypeError("Mask must be list") if not ( @@ -536,13 +554,12 @@ def sparse_mask(self, mask): if mask == []: self._mask = [] else: - lines = np.hstack([self.corners, np.roll(self.corners, -1, axis=0)]) + lines = np.hstack( + [self.corners, np.roll(self.corners, -1, axis=0)]) minx = np.min(lines[:, 0:4:2]) maxx = np.max(lines[:, 0:4:2]) miny = np.min(lines[:, 1:4:2]) maxy = np.max(lines[:, 1:4:2]) - nx = maxx - minx + 1 # inclusive - ny = maxy - miny + 1 sparse_mask = np.asarray(mask) sel = np.logical_and( np.logical_and( @@ -570,9 +587,8 @@ def project_regions( and hasattr(regions_list, "__len__") and len(regions_list) == len(regional_data_list) ): - raise TypeError( - "Region data list and regions lists must be lists of equal length" - ) + raise TypeError("Region data list and regions lists " + "must be lists of equal length") if not all([isinstance(x, np.ndarray) for x in regional_data_list]): raise TypeError("Region data list must be a list of ndarrays") if not all([isinstance(x, BoundingBox) for x in regions_list]): @@ -614,9 +630,9 @@ def project_regions( fnx = xu - xl + 1 # inclusive fny = yu - yl + 1 # inclusive if f.shape[axes[0]] != fny - 1 or f.shape[axes[1]] != fnx - 1: - raise ValueError( - "One or more bounding box descriptors does not match shape of corresponding data cubes" - ) + raise ValueError("One or more bounding box descriptors " + "does not match shape of corresponding " + "data cubes") slc_data = [slice(None)] * len(stitched_img.shape) for (start, end), axis in zip([(yl, yu), (xl, xu)], axes): slc_data[axis] = slice(start, end) @@ -642,9 +658,8 @@ def AxisAlignedBoundingBox( ): """ Constructs an axis aligned bounding box around convex hull """ if not isinstance(convex_hull_object, BoundingConvexHull): - raise TypeError( - "Convex hull object passed in constructor is not of type BoundingConvexHull" - ) + raise TypeError("Convex hull object passed in constructor " + "is not of type BoundingConvexHull") if square: nx = ( np.max(convex_hull_object.corners[:, 0]) @@ -688,7 +703,8 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): if not isinstance(bounding_box_object, BoundingBox): raise TypeError("Expected bounding box object") if not (isinstance(nsubboxes, int) and nsubboxes >= 1): - raise ValueError("nsubboxes must be integral type and be 1 or more") + raise ValueError( + "nsubboxes must be integral type and be 1 or more") xl = np.min(bounding_box_object.corners[:, 0]) xu = np.max(bounding_box_object.corners[:, 0]) yl = np.min(bounding_box_object.corners[:, 1]) @@ -713,7 +729,8 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): xus = xus - 1 yus = yus - 1 - # clamp the final coordinate to the upper end (may result in rectanglular box at the end) + # clamp the final coordinate to the upper end + # (may result in rectanglular box at the end) xus[:, -1] = max(xu, min(xus[0, -1], xu)) yus[-1, :] = max(yu, min(yus[-1, 0], yu)) @@ -747,7 +764,10 @@ def SplitBox(cls, bounding_box_object, nsubboxes=1, **kwargs): @classmethod def PadBox(cls, bounding_box_object, desired_nx, desired_ny, **kwargs): - """ Creates a box with a padded border around a axis-aligned bounding box """ + """ + Creates a box with a padded border + around a axis-aligned bounding box + """ if not isinstance(bounding_box_object, BoundingBox): raise TypeError("Expected bounding box object") nx, ny = bounding_box_object.box_npx diff --git a/africanus/linalg/test/test_geometry.py b/africanus/linalg/test/test_geometry.py index 48809c2f8..a50cd2831 100644 --- a/africanus/linalg/test/test_geometry.py +++ b/africanus/linalg/test/test_geometry.py @@ -34,18 +34,12 @@ def test_hull_construction(debug): sinc_npx = 255 sinc = np.sinc(np.linspace(-7, 7, sinc_npx)) sinc2d = np.outer(sinc, sinc).reshape((1, 1, sinc_npx, sinc_npx)) - extracted_data, extracted_window_extents = BoundingConvexHull.regional_data( + (extracted_data, + extracted_window_extents) = BoundingConvexHull.regional_data( bh_extract, sinc2d, oob_value=np.nan ) assert extracted_window_extents == [-10, 293, -30, 268] sparse_mask = np.array(bh_extract.sparse_mask) - lines = np.hstack( - [bh_extract.corners, np.roll(bh_extract.corners, -1, axis=0)] - ) - minx = np.min(lines[:, 0:4:2]) - maxx = np.max(lines[:, 0:4:2]) - miny = np.min(lines[:, 1:4:2]) - maxy = np.max(lines[:, 1:4:2]) sel = np.logical_and( np.logical_and(sparse_mask[:, 1] >= 0, sparse_mask[:, 1] < 255), np.logical_and(sparse_mask[:, 0] >= 0, sparse_mask[:, 0] < 255), From e61e5eeee97dbf32279f5c050706977019dd9c66 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:55:50 +0200 Subject: [PATCH 07/10] [skip ci] Claim the PR --- HISTORY.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/HISTORY.rst b/HISTORY.rst index 582fc0f13..9e743596b 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -4,6 +4,7 @@ History X.Y.Z (YYYY-MM-DD) ------------------ +* Bounding hull geometric and image manipulation routines (:pr:`192`, :pr:`154`) * Fix spectral model documentation (:pr:`190`), to match changes in (:pr:`189`) 0.2.3 (2020-05-14) From cf9e69940b5df73ef7b18100a61c81c015988474 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:59:47 +0200 Subject: [PATCH 08/10] Fix optional decorator location --- africanus/linalg/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py index 0d95a6b69..f935c36e8 100644 --- a/africanus/linalg/geometry.py +++ b/africanus/linalg/geometry.py @@ -30,8 +30,8 @@ def timed(*args, **kw): return timed -@requires_optional("scipy.stats", opt_import_err) class BoundingConvexHull(object): + @requires_optional("scipy.stats", opt_import_err) def __init__(self, list_hulls, name="unnamed", mask=None, check_mask_outofbounds=True): """ From de6dcedfeefbb5854abe9eba95d3746ca10722a3 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 13:59:59 +0200 Subject: [PATCH 09/10] Fix deprecation warning --- africanus/linalg/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/africanus/linalg/geometry.py b/africanus/linalg/geometry.py index f935c36e8..934e73598 100644 --- a/africanus/linalg/geometry.py +++ b/africanus/linalg/geometry.py @@ -551,7 +551,7 @@ def sparse_mask(self, mask): ) ): raise TypeError("Mask must be a sparse mask of 2 element values") - if mask == []: + if len(mask) == 0: self._mask = [] else: lines = np.hstack( From 79aa25b9cd29ed10c8b7ef67548fa85a42137980 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Fri, 22 May 2020 14:04:51 +0200 Subject: [PATCH 10/10] Remove requirements.py re-introduced in merge --- africanus/install/requirements.py | 47 ------------------------------- 1 file changed, 47 deletions(-) delete mode 100644 africanus/install/requirements.py diff --git a/africanus/install/requirements.py b/africanus/install/requirements.py deleted file mode 100644 index 2dfd20084..000000000 --- a/africanus/install/requirements.py +++ /dev/null @@ -1,47 +0,0 @@ -# -*- coding: utf-8 -*- - - -# NOTE(sjperkins) -# Non standard library imports should be avoided, -# or should fail gracefully as functionality -# in these modules is called by setup.py -import os - -# requirements -on_rtd = os.environ.get('READTHEDOCS') == 'True' - -# Basic requirements containing no C extensions. -# This is necessary for building on RTD -requirements = ['appdirs >= 1.4.3', - 'decorator'] - -if not on_rtd: - requirements += [ - # astropy breaks with numpy 1.15.3 - # https://github.com/astropy/astropy/issues/7943 - 'numpy >= 1.14.0, != 1.15.3', - 'numba >= 0.43.0'] - -extras_require = { - 'cuda': ['cupy >= 5.0.0', 'jinja2 >= 2.10'], - 'dask': ['dask[array] >= 1.1.0'], - 'jax': ['jax == 0.1.27', 'jaxlib == 0.1.14'], - 'scipy': ['scipy >= 1.0.0'], - 'astropy': ['astropy >= 2.0.0, < 3.0; python_version <= "2.7"', - 'astropy >= 3.0; python_version >= "3.0"'], - 'python-casacore': ['python-casacore == 3.0.0'], - 'testing': ['pytest', 'flaky', 'pytest-flake8', 'pathlib'] -} - -_non_cuda_extras = [er for n, er in extras_require.items() if n != "cuda"] -_all_extras = extras_require.values() - -extras_require['complete'] = sorted(set(sum(_non_cuda_extras, []))) -extras_require['complete-cuda'] = sorted(set(sum(_all_extras, []))) - -setup_requirements = [] -test_requirements = (extras_require['testing'] + - extras_require['astropy'] + - extras_require['python-casacore'] + - extras_require['dask'] + - extras_require['scipy'])