Python skimage.measure.find_contours() Examples

def critical_curve_via_magnification_from_tracer_and_grid(tracer, grid):
    magnification = tracer.magnification_from_grid(grid=grid)

    inverse_magnification = 1 / magnification

    critical_curves_indices = measure.find_contours(inverse_magnification.in_2d, 0)

    no_critical_curves = len(critical_curves_indices)
    contours = []
    critical_curves = []

    for jj in np.arange(no_critical_curves):
        contours.append(critical_curves_indices[jj])
        contour_x, contour_y = contours[jj].T
        pixel_coord = np.stack((contour_x, contour_y), axis=-1)

        critical_curve = grid.geometry.grid_scaled_from_grid_pixels_1d_for_marching_squares(
            grid_pixels_1d=pixel_coord, shape_2d=magnification.sub_shape_2d
        )

        critical_curve = np.array(grid=critical_curve)

        critical_curves.append(critical_curve)

    return critical_curves 

def binary_mask_to_polygon(binary_mask, tolerance=0):
    """Converts a binary mask to COCO polygon representation
    Args:
        binary_mask: a 2D binary numpy array where '1's represent the object
        tolerance: Maximum distance from original points of polygon to approximated
            polygonal chain. If tolerance is 0, the original coordinate array is returned.
    """
    polygons = []
    # pad mask to close contours of shapes which start and end at an edge
    padded_binary_mask = np.pad(binary_mask, pad_width=1, mode='constant', constant_values=0)
    contours = measure.find_contours(padded_binary_mask, 0.5)
    contours = np.subtract(contours, 1)
    for contour in contours:
        contour = close_contour(contour)
        contour = measure.approximate_polygon(contour, tolerance)
        if len(contour) < 3:
            continue
        contour = np.flip(contour, axis=1)
        segmentation = contour
        # after padding and subtracting 1 we may get -0.5 points in our segmentation
        segmentation = [np.clip(i,0.0,i).tolist() for i in segmentation]
        polygons.append(segmentation)

    return polygons 

def binary_mask_to_polygon(binary_mask, tolerance=0):
    """Converts a binary mask to COCO polygon representation
    Args:
        binary_mask: a 2D binary numpy array where '1's represent the object
        tolerance: Maximum distance from original points of polygon to approximated
            polygonal chain. If tolerance is 0, the original coordinate array is returned.
    """
    polygons = []
    # pad mask to close contours of shapes which start and end at an edge
    padded_binary_mask = np.pad(binary_mask, pad_width=1, mode='constant', constant_values=0)
    contours = measure.find_contours(padded_binary_mask, 0.5)
    contours = np.subtract(contours, 1)
    for contour in contours:
        contour = close_contour(contour)
        contour = measure.approximate_polygon(contour, tolerance)
        if len(contour) < 3:
            continue
        contour = np.flip(contour, axis=1)
        segmentation = contour.ravel().tolist()
        # after padding and subtracting 1 we may get -0.5 points in our segmentation
        segmentation = [0 if i < 0 else i for i in segmentation]
        polygons.append(segmentation)

    return polygons 

def vectorize_regions(im: np.ndarray, threshold: float = 0.5):
    """
    Vectorizes lines from a binarized array.

    Args:
        im (np.ndarray): Array of shape (H, W) with the first dimension
                         being a probability distribution over the region.
        threshold (float): Threshold for binarization

    Returns:
        [[x0, y0, ... xn, yn], [xm, ym, ..., xk, yk], ... ]
        A list of lists containing the region polygons.
    """
    bin = im > threshold
    contours = find_contours(bin, 0.5, fully_connected='high', positive_orientation='high')
    if len(contours) == 0:
        return contours
    approx_contours = []
    for contour in contours:
        approx_contours.append(approximate_polygon(contour[:,[1,0]], 1).astype('uint').tolist())
    return approx_contours 

def __call__(self, img_binary):
        cs = measure.find_contours(img_binary, 0.5)

        # Collect contours into RotatedBoxes
        results = []
        for c in cs:
            # Now examine the bounding box. If it is too small, we ignore the contour
            ll, ur = np.min(c, 0), np.max(c, 0)
            wh = ur - ll
            if wh[0] * wh[1] < self.min_area:
                continue

            # Finally, construct the rotatedbox. If its aspect ratio is too small, we ignore it
            rb = RotatedBox.from_points(c, self.box_type)
            if rb.height == 0 or rb.width / rb.height < self.min_box_aspect:
                continue

            # All tests fine, add to the list
            results.append(rb)

        # Next sort and leave only max_boxes largest boxes by area
        results.sort(key=lambda x: -x.area)
        return self._merge_boxes(results[0:self.max_boxes]) 

def make_contours(self, level=0.):
        """
        This function returns the contours as a List of coordinate positions
        at one given level - run this more than once with different levels 
        if desired

        :type level: float
        :param level: Level (z value) of grid to contour

        :return:
            contours: List of contours with coordinate positions

        """
        try:
            from skimage import measure
        except:
            raise(Exception("Package 'scikit-image' not available to make contours."))

        return measure.find_contours(self.data, level) 

def level_curves(fname, npoints = 200, smoothing = 10, level = 0.5) :
	"Loads regularly sampled curves from a .PNG image."
	# Find the contour lines
	img = misc.imread(fname, flatten = True) # Grayscale
	img = (img.T[:, ::-1])  / 255.
	img = gaussian_filter(img, smoothing, mode='nearest')
	lines = find_contours(img, level)
	
	# Compute the sampling ratio for every contour line
	lengths = np.array( [arclength(line) for line in lines] )
	points_per_line = np.ceil( npoints * lengths / np.sum(lengths) )
	
	# Interpolate accordingly
	points = [] ; connec = [] ; index_offset = 0
	for ppl, line in zip(points_per_line, lines) :
		(p, c) = resample(line, ppl)
		points.append(p)
		connec.append(c + index_offset)
		index_offset += len(p)
	
	size   = np.maximum(img.shape[0], img.shape[1])
	points = np.vstack(points) / size
	connec = np.vstack(connec)
	return Curve(points, connec)
# Pyplot Output ================================================================================= 

def level_curves(fname, npoints = 200, smoothing = 10, level = 0.5) :
	"Loads regularly sampled curves from a .PNG image."
	# Find the contour lines
	img = misc.imread(fname, flatten = True) # Grayscale
	img = (img.T[:, ::-1])  / 255.
	img = gaussian_filter(img, smoothing, mode='nearest')
	lines = find_contours(img, level)
	
	# Compute the sampling ratio for every contour line
	lengths = np.array( [arclength(line) for line in lines] )
	points_per_line = np.ceil( npoints * lengths / np.sum(lengths) )
	
	# Interpolate accordingly
	points = [] ; connec = [] ; index_offset = 0
	for ppl, line in zip(points_per_line, lines) :
		(p, c) = resample(line, ppl)
		points.append(p)
		connec.append(c + index_offset)
		index_offset += len(p)
	
	size   = np.maximum(img.shape[0], img.shape[1])
	points = np.vstack(points) / size
	connec = np.vstack(connec)
	return Curve(points, connec)
# Pyplot Output ================================================================================= 

def vectorize(self,  **pixelBlocks):
        image = pixelBlocks['raster_pixels']
        _, height, width = image.shape
        image = np.transpose(image, [1,2,0])

        with self.graph.as_default():
            results = self.model.detect([image], verbose=1)

        masks = results[0]['masks']
        mask_count = masks.shape[-1]
        coord_list = []
        for m in range(mask_count):
            mask = masks[:, :, m]
            padded_mask = np.zeros((mask.shape[0]+2, mask.shape[1]+2), dtype=np.uint8)
            padded_mask[1:-1, 1:-1] = mask
            contours = find_contours(padded_mask, 0.5, fully_connected='high')

            if len(contours) != 0:
                verts = contours[0] - 1
                coord_list.append(verts)

        if self.padding != 0:
            coord_list[:] = [item - self.padding for item in coord_list]

        return coord_list, results[0]['scores'], results[0]['class_ids'] 

def fill_convex (image):
    H, W = image.shape
    padded = np.zeros((H+20, W+20), dtype=np.uint8)
    padded[10:(10+H),10:(10+W)] = image

    contours = measure.find_contours(padded, 0.5)
    if len(contours) == 0:
        return image
    if len(contours) == 1:
        contour = contours[0]
    else:
        contour = np.vstack(contours)
    cc = np.zeros_like(contour, dtype=np.int32)
    cc[:,0] = contour[:, 1]
    cc[:,1] = contour[:, 0]
    hull = cv2.convexHull(cc)
    contour = hull.reshape((1, -1, 2)) 
    cv2.fillPoly(padded, contour, 1)
    return padded[10:(10+H),10:(10+W)] 

def fill_convex (image):
    H, W = image.shape
    padded = np.zeros((H+20, W+20), dtype=np.uint8)
    padded[10:(10+H),10:(10+W)] = image

    contours = measure.find_contours(padded, 0.5)
    if len(contours) == 0:
        return image
    if len(contours) == 1:
        contour = contours[0]
    else:
        contour = np.vstack(contours)
    cc = np.zeros_like(contour, dtype=np.int32)
    cc[:,0] = contour[:, 1]
    cc[:,1] = contour[:, 0]
    hull = cv2.convexHull(cc)
    contour = hull.reshape((1, -1, 2)) 
    cv2.fillPoly(padded, contour, 1)
    return padded[10:(10+H),10:(10+W)] 

def get_points(mask, spacing):
    points = []
    if np.amax(mask) > 0:
        contours = measure.find_contours(mask.astype(np.float32), 0.5)
        for contour in contours:
            con = contour.ravel().reshape((-1, 2)) # con[0] is along x (last dimension of mask)
            for p in con:
                points.append([p[1] * spacing[0], p[0] * spacing[1]])
    return np.asarray(points, dtype=np.float32) 

def getRegressionData(self, mapDim, numContours=15, coreBlack=True):
        masked = self.getFaceImage(mapDim)
        regressionMap = None
        contours = []
        contourLengths = {}

        try:
            cellSize = 1 / mapDim
            regressionMap = skfmm.distance(masked, dx=cellSize) * 2
            regmax = np.amax(regressionMap)
            regressionMap = regressionMap[:, :].copy()
            if coreBlack:
                regressionMap[np.where(masked == 0)] = regmax # Make the core black

            for dist in np.linspace(0, regmax, numContours):
                contours.append([])
                contourLengths[dist] = 0
                layerContours = measure.find_contours(regressionMap, dist, fully_connected='high')
                for contour in layerContours:
                    contours[-1].append(contour)
                    contourLengths[dist] += geometry.length(contour, mapDim)

        except ValueError as exc: # If there aren't any contours, do nothing
            print(exc)

        return (masked, regressionMap, contours, contourLengths) 

def binary_mask_to_polygon(binary_mask, tolerance=0):
        """Converts a binary mask to COCO polygon representation

         :param binary_mask: a 2D binary numpy array where '1's represent the object
         :param tolerance: Maximum distance from original points of polygon to approximated polygonal chain. If tolerance is 0, the original coordinate array is returned.
        """
        polygons = []
        # pad mask to close contours of shapes which start and end at an edge
        padded_binary_mask = np.pad(binary_mask, pad_width=1, mode='constant', constant_values=0)
        contours = np.array(measure.find_contours(padded_binary_mask, 0.5))
        # Reverse padding
        contours = contours - 1
        for contour in contours:
            # Make sure contour is closed
            contour = CocoUtility.close_contour(contour)
            # Approximate contour by polygon
            polygon = measure.approximate_polygon(contour, tolerance)
            # Skip invalid polygons
            if len(polygon) < 3:
                continue
            # Flip xy to yx point representation
            polygon = np.flip(polygon, axis=1)
            # Flatten
            polygon = polygon.ravel()
            # after padding and subtracting 1 we may get -0.5 points in our segmentation
            polygon[polygon < 0] = 0
            polygons.append(polygon.tolist())

        return polygons 

def find_contours(*args):
    
    contours = FC(*args)
    
    simplified_contours = [np.array(simplify_coords(x, 1), dtype=np.int32) \
                                for x in contours]
    
    return simplified_contours 

def level_curves(fname, npoints, smoothing = 10, level = 0.5) :
	# Find the contour lines
	img = misc.imread(fname, flatten = True) # Grayscale
	img = img.T[:, ::-1]
	img = img / 255.
	img = gaussian_filter(img, smoothing, mode='nearest')
	lines = find_contours(img, level)
	
	# Compute the sampling ratio
	lengths = []
	for line in lines :
		lengths.append( arclength(line) )
	lengths = array(lengths)
	points_per_line = ceil( npoints * lengths / sum(lengths) )
	
	# Interpolate accordingly
	points = []
	connec = []
	index_offset = 0
	for ppl, line in zip(points_per_line, lines) :
		(p, c) = resample(line, ppl)
		points.append(p)
		connec.append(c + index_offset)
		index_offset += len(p)
	
	points = vstack(points)
	connec = vstack(connec)
	return Curve(points.ravel(), connec, 2) # Dimension 2 ! 

def getRegressionData(self, mapDim, numContours=15, coreBlack=True):
        self.initGeometry(mapDim)
        self.generateCoreMap()

        masked = np.ma.MaskedArray(self.coreMap, self.mask)
        regressionMap = None
        contours = []
        contourLengths = {}

        try:
            self.generateRegressionMap()

            regmax = np.amax(self.regressionMap)

            regressionMap = self.regressionMap[:, :].copy()
            if coreBlack:
                regressionMap[np.where(self.coreMap == 0)] = regmax # Make the core black
            regressionMap = np.ma.MaskedArray(regressionMap, self.mask)

            for dist in np.linspace(0, regmax, numContours):
                contours.append([])
                contourLengths[dist] = 0
                layerContours = measure.find_contours(self.regressionMap, dist, fully_connected='low')
                for contour in layerContours:
                    contours[-1].append(geometry.clean(contour, self.mapDim, 3))
                    contourLengths[dist] += geometry.length(contour, self.mapDim)

        except ValueError as exc: # If there aren't any contours, do nothing
            print(exc)

        return (masked, regressionMap, contours, contourLengths) 

def getCorePerimeter(self, regDist):
        mapDist = self.normalize(regDist)

        corePerimeter = 0
        contours = measure.find_contours(self.regressionMap, mapDist, fully_connected='low')
        for contour in contours:
            corePerimeter += self.mapToLength(geometry.length(contour, self.mapDim))

        return corePerimeter 

def getRegressionData(self, mapDim, numContours=15, coreBlack=True):
        masked = self.getFaceImage(mapDim)
        regressionMap = None
        contours = []
        contourLengths = {}

        try:
            cellSize = 1 / mapDim
            regressionMap = skfmm.distance(masked, dx=cellSize) * 2
            regmax = np.amax(regressionMap)
            regressionMap = regressionMap[:, :].copy()
            if coreBlack:
                regressionMap[np.where(masked == 0)] = regmax # Make the core black

            for dist in np.linspace(0, regmax, numContours):
                contours.append([])
                contourLengths[dist] = 0
                layerContours = measure.find_contours(regressionMap, dist, fully_connected='high')
                for contour in layerContours:
                    contours[-1].append(contour)
                    contourLengths[dist] += geometry.length(contour, mapDim)

        except ValueError as exc: # If there aren't any contours, do nothing
            print(exc)

        return (masked, regressionMap, contours, contourLengths) 

def convex_hull (binary):
    swap_sequence = [(0, 1),  # 102
                     (0, 2),  # 201
                     (0, 2)]  # 102

    output = np.ndarray(binary.shape, dtype=binary.dtype)
    for swp1, swp2 in swap_sequence:
        N = binary.shape[0]
        print 'shape', binary.shape
        for i in range(N):
            contours = measure.find_contours(binary[i], 0.5)
            if len(contours) == 0:
                continue
            if len(contours) == 1:
                contour = contours[0]
            else:
                contour = np.vstack(contours)
            cc = np.zeros_like(contour, dtype=np.int32)
            cc[:,0] = contour[:, 1]
            cc[:,1] = contour[:, 0]
            hull = cv2.convexHull(cc)
            contour = hull.reshape((1, -1, 2)) 
            cv2.fillPoly(binary[i], contour, 1)
            #binary[i] = skimage.morphology.convex_hull_image(binary[i])
            pass
        print 'swap', swp1, swp2
        nb = np.swapaxes(binary, swp1, swp2)
        binary = np.ndarray(nb.shape, dtype=nb.dtype)
        binary[:,:] = nb[:,:]
        pass
    binary = np.swapaxes(binary, 0, 1)
    output[:,:] = binary[:,:]
    return output;
    #binary = binary_dilation(output, iterations=dilate)
    #return binary 

def _calculate_contours_centroids(self):
        cnts_df_list = []
        cts_df_list = []
        A = self.cnmf['A'].load()
        for uid in range(self._u):
            cur_A = A.sel(unit_id=uid)
            cur_idxs = cur_A.squeeze().dims
            cur_thres = dask.delayed(cur_A.max)()
            cur_thres = dask.delayed(float)(cur_thres * .3)
            cur_cnts = dask.delayed(find_contours)(cur_A, cur_thres)
            cur_cnts = dask.delayed(np.concatenate)(cur_cnts)
            cur_cnts = dask.delayed(pd.DataFrame)(cur_cnts, columns=cur_idxs)
            cur_cnts = cur_cnts.assign(unit_id=uid)
            cur_cts = dask.delayed(center_of_mass)(cur_A.values)
            cur_cts = dask.delayed(pd.Series)(cur_cts, index=cur_idxs)
            cur_cts = cur_cts.append(pd.Series(dict(unit_id=uid)))
            cnts_df_list.append(cur_cnts)
            cts_df_list.append(cur_cts)
        cnts_df_list = dask.compute(*cnts_df_list)
        cts_df_list = dask.compute(*cts_df_list)
        cnts_df = pd.concat(cnts_df_list)
        cts_df = pd.concat(cts_df_list, axis=1).T
        for dim in cur_idxs:
            cnts_df[dim].update(cnts_df[dim] / A.sizes[dim] * self._dims[dim])
            cts_df[dim].update(cts_df[dim] / A.sizes[dim] * self._dims[dim])
        return cnts_df, cts_df 

def binary_mask_to_polygon(binary_mask, tolerance=0):
    """Converts a binary mask to COCO polygon representation

    Args:
        binary_mask: a 2D binary numpy array where '1's represent the object
        tolerance: Maximum distance from original points of polygon to approximated
            polygonal chain. If tolerance is 0, the original coordinate array is returned.

    """
    polygons = []
    # pad mask to close contours of shapes which start and end at an edge
    padded_binary_mask = np.pad(binary_mask, pad_width=1, mode='constant', constant_values=0)
    contours = measure.find_contours(padded_binary_mask, 0.5)
    contours = np.subtract(contours, 1)
    for contour in contours:
        contour = close_contour(contour)
        contour = measure.approximate_polygon(contour, tolerance)
        if len(contour) < 3:
            continue
        contour = np.flip(contour, axis=1)
        segmentation = contour.ravel().tolist()
        # after padding and subtracting 1 we may get -0.5 points in our segmentation 
        segmentation = [0 if i < 0 else i for i in segmentation]
        polygons.append(segmentation)

    return polygons 

def show_fig2():
    contours = measure.find_contours(phi, 0)
    ax2 = fig2.add_subplot(111)
    ax2.imshow(img, interpolation='nearest', cmap=plt.cm.gray)
    for n, contour in enumerate(contours):
        ax2.plot(contour[:, 1], contour[:, 0], linewidth=2) 

def convex_hull (binary):
    swap_sequence = [(0, 1),  # 102
                     (0, 2),  # 201
                     (0, 2)]  # 102

    output = np.ndarray(binary.shape, dtype=binary.dtype)
    for swp1, swp2 in swap_sequence:
        N = binary.shape[0]
        print 'shape', binary.shape
        for i in range(N):
            contours = measure.find_contours(binary[i], 0.5)
            if len(contours) == 0:
                continue
            if len(contours) == 1:
                contour = contours[0]
            else:
                contour = np.vstack(contours)
            cc = np.zeros_like(contour, dtype=np.int32)
            cc[:,0] = contour[:, 1]
            cc[:,1] = contour[:, 0]
            hull = cv2.convexHull(cc)
            contour = hull.reshape((1, -1, 2)) 
            cv2.fillPoly(binary[i], contour, 1)
            #binary[i] = skimage.morphology.convex_hull_image(binary[i])
            pass
        print 'swap', swp1, swp2
        nb = np.swapaxes(binary, swp1, swp2)
        binary = np.ndarray(nb.shape, dtype=nb.dtype)
        binary[:,:] = nb[:,:]
        pass
    binary = np.swapaxes(binary, 0, 1)
    output[:,:] = binary[:,:]
    return output;
    #binary = binary_dilation(output, iterations=dilate)
    #return binary 

def find_wl_contours1(im_ndwi, cloud_mask, im_ref_buffer):
    """
    Traditional method for shoreline detection using a global threshold.
    Finds the water line by thresholding the Normalized Difference Water Index 
    and applying the Marching Squares Algorithm to contour the iso-value 
    corresponding to the threshold.

    KV WRL 2018

    Arguments:
    -----------
    im_ndwi: np.ndarray
        Image (2D) with the NDWI (water index)
    cloud_mask: np.ndarray
        2D cloud mask with True where cloud pixels are
    im_ref_buffer: np.array
        Binary image containing a buffer around the reference shoreline

    Returns:    
    -----------
    contours_wl: list of np.arrays
        contains the coordinates of the contour lines

    """

    # reshape image to vector
    vec_ndwi = im_ndwi.reshape(im_ndwi.shape[0] * im_ndwi.shape[1])
    vec_mask = cloud_mask.reshape(cloud_mask.shape[0] * cloud_mask.shape[1])
    vec = vec_ndwi[~vec_mask]
    # apply otsu's threshold
    vec = vec[~np.isnan(vec)]
    t_otsu = filters.threshold_otsu(vec)
    # use Marching Squares algorithm to detect contours on ndwi image
    im_ndwi_buffer = np.copy(im_ndwi)
    im_ndwi_buffer[~im_ref_buffer] = np.nan
    contours = measure.find_contours(im_ndwi_buffer, t_otsu)

    # remove contours that contain NaNs (due to cloud pixels in the contour)
    contours_nonans = []
    for k in range(len(contours)):
        if np.any(np.isnan(contours[k])):
            index_nan = np.where(np.isnan(contours[k]))[0]
            contours_temp = np.delete(contours[k], index_nan, axis=0)
            if len(contours_temp) > 1:
                contours_nonans.append(contours_temp)
        else:
            contours_nonans.append(contours[k])
    contours = contours_nonans

    return contours 

def split_letters(image, num_letters=6, debug=False):
    '''
    split full captcha image into `num_letters` lettersself.
    return list of letters binary image (0: white, 255: black)
    '''
    # move left
    left = crop(image, ((0, 0), (0, image.shape[1]-SHIFT_PIXEL)), copy=True)
    image[:,:-SHIFT_PIXEL] = image[:,SHIFT_PIXEL:]
    image[:,-SHIFT_PIXEL:] = left
    # binarization
    binary = image > BINARY_THRESH
    # find contours
    contours = find_contours(binary, 0.5)
    contours = [[
        [int(floor(min(contour[:, 1]))), int(floor(min(contour[:, 0])))], # top-left point
        [int(ceil(max(contour[:, 1]))), int(ceil(max(contour[:, 0])))]  # down-right point
      ] for contour in contours]
    # keep letters order
    contours = sorted(contours, key=lambda contour: contour[0][0])
    # find letters box
    letter_boxs = []
    for contour in contours:
        if len(letter_boxs) > 0 and contour[0][0] < letter_boxs[-1][1][0] - 5:
            # skip inner contour
            continue
        # extract letter boxs by contour
        boxs = get_letter_boxs(binary, contour)
        for box in boxs:
            letter_boxs.append(box)
    # check letter outer boxs number
    if len(letter_boxs) != num_letters:
        print('ERROR: number of letters is NOT valid', len(letter_boxs))
        # debug
        if debug:
            print(letter_boxs)
            plt.imshow(binary, interpolation='nearest', cmap=plt.cm.gray)
            for [x_min, y_min], [x_max, y_max] in letter_boxs:
                plt.plot(
                    [x_min, x_max, x_max, x_min, x_min],
                    [y_min, y_min, y_max, y_max, y_min],
                    linewidth=2)
            plt.xticks([])
            plt.yticks([])
            plt.show()
        return None

    # normalize size (40x40)
    letters = []
    for [x_min, y_min], [x_max, y_max] in letter_boxs:
        letter = resize(image[y_min:y_max, x_min:x_max], LETTER_SIZE)
        letter = img_as_ubyte(letter < 0.6)
        letters.append(letter)

    return letters 

def _mask_to_polygon(mask, gdir=None):
    """Converts a mask to a single polygon.

    The mask should be a single entity with nunataks: I didnt test for more
    than one "blob".

    Parameters
    ----------
    mask: 2d array with ones and zeros
        the mask to convert
    gdir: GlacierDirectory
        for logging

    Returns
    -------
    (poly, poly_no_nunataks) Shapely polygons
    """

    regions, nregions = label(mask, structure=LABEL_STRUCT)
    if nregions > 1:
        rid = ''
        if gdir is not None:
            rid = gdir.rgi_id
        log.debug('(%s) we had to cut a blob from the catchment', rid)
        # Check the size of those
        region_sizes = [np.sum(regions == r) for r in np.arange(1, nregions+1)]
        am = np.argmax(region_sizes)
        # Check not a strange glacier
        sr = region_sizes.pop(am)
        for ss in region_sizes:
            if (ss / sr) > 0.2:
                log.info('(%s) this blob was unusually large', rid)
        mask[:] = 0
        mask[np.where(regions == (am+1))] = 1

    nlist = measure.find_contours(mask, 0.5)
    # First is the exterior, the rest are nunataks
    e_line = shpg.LinearRing(nlist[0][:, ::-1])
    i_lines = [shpg.LinearRing(ipoly[:, ::-1]) for ipoly in nlist[1:]]

    poly = shpg.Polygon(e_line, i_lines).buffer(0)
    if not poly.is_valid:
        raise GeometryError('Mask to polygon conversion error.')
    poly_no = shpg.Polygon(e_line).buffer(0)
    if not poly_no.is_valid:
        raise GeometryError('Mask to polygon conversion error.')
    return poly, poly_no 

def bathymetry(projection, x, y, z, args):
    lat, lon = get_latlon_coords(projection, x, y, z)
    if len(lat.shape) == 1:
        lat, lon = np.meshgrid(lat, lon)

    xpx = x * 256
    ypx = y * 256

    with Dataset(current_app.config['ETOPO_FILE'] % (projection, z), 'r') as dataset:
        data = dataset["z"][ypx:(ypx + 256), xpx:(xpx + 256)] * -1
        data = data[::-1, :]

    LEVELS = [100, 200, 500, 1000, 2000, 3000, 4000, 5000, 6000]

    normalized = matplotlib.colors.LogNorm(vmin=1, vmax=6000)(LEVELS)
    cmap = matplotlib.colors.LinearSegmentedColormap.from_list(
        'transparent_gray',
        [(0, 0, 0, 1), (0, 0, 0, 0.5)]
    )
    colors = cmap(normalized)

    fig = plt.figure()
    fig.set_size_inches(4, 4)
    ax = plt.Axes(fig, [0, 0, 1, 1])
    ax.set_axis_off()
    fig.add_axes(ax)

    for i, l in enumerate(LEVELS):
        contours = measure.find_contours(data, l)

        for n, contour in enumerate(contours):
            ax.plot(contour[:, 1], contour[:, 0], color=colors[i], linewidth=1)

    plt.xlim([0, 255])
    plt.ylim([0, 255])

    with contextlib.closing(BytesIO()) as buf:
        plt.savefig(
            buf,
            format='png',
            dpi=64,
            transparent=True,
        )
        plt.close(fig)
        buf.seek(0)
        im = Image.open(buf)

        buf2 = BytesIO()
        im.save(buf2, format='PNG', optimize=True)
        return buf2

    return None 

def computeShoreline(self,z,lvl=0.):
        """
        This function computes the shoreline position for a given sea-level.

        Args:
            z: mesh relative elevation to sea-level.
            lvl: water level defined in the input.

        Returns:
            - contourPts - numpy array containing the contour coordinates.
        """

        contours = measure.find_contours(z.T,level=lvl)
        contourList = []
        start = True

        # Loop through each contour
        for c in range(len(contours)):
            tmpts =  contours[c]
            tmpts[:,0] = tmpts[:,0]*self.dx+self.xi.min()
            tmpts[:,1] = tmpts[:,1]*self.dx+self.yi.min()
            closed = False
            if tmpts[0,0] == tmpts[-1,0] and tmpts[0,1] == tmpts[-1,1]:
                closed = True

            # Remove duplicate points
            unique = OrderedDict()
            for p in zip(tmpts[:,0], tmpts[:,1]):
                unique.setdefault(p[:2], p)
            pts = numpy.asarray(list(unique.values()))

            if closed:
                cpts = numpy.zeros((len(pts)+1,2), order='F')
                cpts[0:len(pts),0:2] = pts
                cpts[-1,0:2] = pts[0,0:2]

                # Get contour length
                arr = cpts
                val = (arr[:-1,:] - arr[1:,:]).ravel()
                dist = val.reshape((arr.shape[0]-1,2))
                lgth = numpy.sum(numpy.sqrt(numpy.sum(dist**2, axis=1)))
            else:
                lgth = 1.e8
                cpts = pts

            if len(cpts) > 2 and lgth > self.mlen:
                contourList.append(cpts)
                if start:
                    contourPts = cpts
                    start = False
                else:
                    contourPts = numpy.concatenate((contourPts,cpts))

        return contourPts 

def mask2poly(mask, threshold=0.5):
    """Takes a mask and returns a MultiPolygon

    Parameters
    ----------
    mask : array
        Sparse or dense array to identify polygon contours within.
    threshold : float, optional
        Threshold value used to separate points in and out of resulting
        polygons. 0.5 will partition a boolean mask, for an arbitrary value
        binary mask choose the midpoint of the low and high values.

    Output
    ------
    MultiPolygon
        Returns a MultiPolygon of all masked regions.

    """
    mask = _reformat_mask(mask)

    verts_list = []
    for z, m in enumerate(mask):
        if issparse(m):
            m = np.array(m.astype('byte').todense())

        if (m != 0).sum() == 0:
            # If the plane is empty, just skip it
            continue

        # Add an empty row and column around the mask to make sure edge masks
        # are correctly determined
        expanded_dims = (m.shape[0] + 2, m.shape[1] + 2)
        expanded_mask = np.zeros(expanded_dims, dtype=float)
        expanded_mask[1:m.shape[0] + 1, 1:m.shape[1] + 1] = m

        verts = find_contours(expanded_mask.T, threshold)

        # Subtract off 1 to shift coords back to their real space,
        # but also add 0.5 to move the coordinates back to the corners,
        # so net subtract 0.5 from every coordinate
        verts = [np.subtract(x, 0.5).tolist() for x in verts]

        v = []
        for poly in verts:
            new_poly = [point + [z] for point in poly]
            v.append(new_poly)
        verts_list.extend(v)

    return _reformat_polygons(verts_list)