def cluster_process(labels, original, activations):
    rbase = np.zeros(labels.shape)
    rubase = np.zeros(labels.shape)
    rubase[range(0,20),:] = 1
    rubase[:,range(0,20)] = 1
    rubase[range(-20,-1),:] = 1
    rubase[:,range(-20,-1)] = 1
    for i in range(1, int(np.max(labels))+1):
        base = np.zeros(labels.shape)
        base[labels==i] = 1
        li = len(base.nonzero()[0])
        if li>0:
            hull = convex_hull_image(base)
            lh =len(hull.nonzero()[0])
            sel_org = base*original
            sel_act = base*activations
            cond = (li > 4000 and float(lh) / float(li) < 1.07 and perimeter(base)**2.0 / li < 30) or np.max(base * rubase) > 0.5
            # print li>4000 and float(lh)/float(li)<1.07, perimeter(base)**2.0/li<30, np.max(base*rubase)>0.5, np.min(original[base>0])
            hard_array =[li > 4000, float(lh) / float(li) < 1.07]
            optional_array = [perimeter(base)**2.0/li < 25,
                              np.percentile(sel_org[sel_org>0], 5) > 0.2,
                              np.percentile(sel_act, 90) - np.percentile(sel_act, 90)]
            print hard_array, optional_array
            if debug and li>1000:
                rs(base,'subspread cluster')
            if cond:
                rbase = rbase + base
    rbase[rubase.astype(np.bool)] = 1
    return dilation(rbase, selem)

def get_bg_mask(img):
    
    #if img.ndim == 3:
    #    bg_mask = img.any(axis=-1)
    #    bg_mask = np.invert(bg_mask) # consistent with np.ma, True if masked

    #    # make multichannel (is it really this hard?)
    #    bg_mask = np.repeat(bg_mask[:,:,np.newaxis], 3, axis=2) 
    #
    #else:
    #    bg_mask = (img != 0)
    #    bg_mask = np.invert(bg_mask) # see above

    #bound = segmentation.find_boundaries(bg_mask, mode='inner', background=1)
    #bg_mask[bound] = 1
    #min_size = img.shape[0] * img.shape[1] // 4 
    #holes = morphology.remove_small_holes(bg_mask, min_size=min_size)
    #bg_mask[holes] = 1
    
    #bg_mask = segmentation.find_boundaries(img)
    #bg_mask = morphology.remove_small_objects(bg_mask)
    #bg_mask = morphology.remove_small_holes(bg_mask)

    bg_mask = morphology.convex_hull_image(img)
    bg_mask = np.zeros_like(img)
    return bg_mask

        def onclick(self,event):
            if event.inaxes==ax:
                if event.button==1:
                    if np.sum(MASKs[self.ind])==0:
                        
                        cx, cy = circle_perimeter(int(event.ydata), int(event.xdata), 3)
                        MASKs[self.ind][cx,cy] = (220, 80, 20, 1)
                        original_image = np.copy(MASKs[self.ind][:,:,3])
                        chull = convex_hull_image(original_image)
                        MASKs[self.ind][chull,:] = (255, 0, 0, .4)


                        ax.cla()
                        data =  areaFile.attrs['ROI_patches'][:,:,self.ind]
                        ax.imshow(data,cmap='binary_r')
                        ax.imshow(MASKs[self.ind])

                elif event.button==3:
                    MASKs[self.ind] = np.zeros(MASKs[self.ind].shape)

                    
                    data = areaFile.attrs['ROI_patches'][:,:,self.ind]

                    ax.cla()
                    ax.imshow(data,cmap='binary_r')
                    ax.imshow(MASKs[self.ind])

    def internal_filter(self, lungs, coronal, body):
        #lungs = ss.get_lungs()
        lungs_sum = np.sum(lungs, axis = 0)>=5
        lungs_hull = convex_hull_image(lungs_sum)
    
        def dist_lungs(lungs_hull, body):
            """ Metoda pro zjisteni, jak blizko jsme stredu (nejvzdalenejsimu mistu od povrchu) """
            blank_body = np.ones(body.shape)
            lungs_edge_dist = scipy.ndimage.morphology.distance_transform_edt(lungs_hull)
            lungs_edge_dist_maximum = np.max(lungs_edge_dist)
            lungs_edge_dist_focus = lungs_edge_dist > 0.2 * lungs_edge_dist_maximum               # 0.1 puvodne
            blank_body[:] = lungs_edge_dist_focus
            ld = scipy.ndimage.morphology.distance_transform_edt(blank_body)
            # resized_to_orig = resize_to_shape(ld, self.ss.orig_shape)
            # return resized_to_orig
            return ld

        dist_lungs_with_body = dist_lungs(lungs_hull, body)
        lungs_mask = (dist_lungs_with_body > 0 ) & (coronal > 0)
        #print_2D_gray(lungs_mask[100])
    
        #print_it_all(ss, data3dr_tmp, lungs_mask*3, "001")
    
    
        def improve_lungs(dist_lungs_hulls, final_area_filter):
            """ Upravi masku lungs pro specialni ucely """
            lungs_hulls = dist_lungs_hulls > 0
            new_filter = np.zeros(final_area_filter.shape)
            for i, area_slice in enumerate(final_area_filter):
                convex = convex_hull_image(area_slice)
                new_filter[i] = lungs_hulls[i] & convex
            return new_filter
        
        return lungs_mask

def shift(img):
    """Shift a binary image randomly within the frame
    
    Uses a convex hull calculation to make sure it doesn't translate
    the image out of the frame.
    """
    hull = morphology.convex_hull_image(1-img)

    horizontal = np.where(np.sum(hull, axis=0) > 0)[0]
    vertical = np.where(np.sum(hull, axis=1) > 0)[0]

    max_left = -np.min(horizontal)
    max_right = img.shape[1] - np.max(horizontal)
    max_down = -np.min(vertical)
    max_up = img.shape[0] - np.max(vertical)
    
    shift_x = np.random.randint(max_left, max_right)
    shift_y = np.random.randint(max_down, max_up)

    #print "SHIFT", shift_x, shift_y
    
    def shift(xy):
        xy[:, 0] -= shift_x
        xy[:, 1] -= shift_y
        return xy
        
    return np.logical_not(transform.warp(np.logical_not(img), shift))

 def deep_struct_filter(self, data3dr_tmp, body):
     """ Najde vsechny hluboke (v ose z) objekty a jejich okoli 5 pixelu """
     
     min_bone_thr = 220
     bone = (data3dr_tmp>min_bone_thr) & body
     bone_sum = np.sum(bone, axis = 0) #> 150
 
     bone_mask = bone_sum>=1
     bone_hull = convex_hull_image(bone_mask)#convex_hull_image(bone_mask)
 
     def weak_dist_bone(self, bone_hull, body):
         """ Metoda pro zjisteni, jak blizko jsme stredu (nejvzdalenejsimu mistu od povrchu) """
         blank_body = np.ones(body.shape)
         bone_edge_dist = scipy.ndimage.morphology.distance_transform_edt(bone_hull)
         bone_edge_dist_maximum = np.max(bone_edge_dist)
         bone_edge_dist_focus = bone_edge_dist > 0*bone_edge_dist_maximum
         blank_body[:] = bone_edge_dist_focus
         ld = scipy.ndimage.morphology.distance_transform_edt(blank_body)
         return resize_to_shape(ld, self.ss.orig_shape)
 
     dh = weak_dist_bone(bone_hull, body)
 
     blur = scipy.ndimage.filters.gaussian_filter(copy.copy(data3dr_tmp), sigma=[17, 3, 3]) > 100
     db = scipy.ndimage.morphology.distance_transform_edt(1-blur)
     dbf = db >= 6
 
     struct_filter = ((dbf==False) & (dh>15))==False
     return struct_filter

def add_auto_masks_area(areaFile,addMasks=False):
    from skimage.morphology import binary_dilation, binary_erosion, disk
    from skimage import exposure
    from skimage.transform import hough_circle
    from skimage.morphology import convex_hull_image
    from skimage.feature import canny
    from skimage.draw import circle_perimeter

    import h5py


    MASKs = np.zeros(areaFile.attrs['ROI_patches'].shape)
    if 'ROI_masks' in (areaFile.attrs.iterkeys()):
        print 'Masks have already been created'
        awns = raw_input('Would you like to redo them from scratch: (answer y/n): ')
    else:
        awns = 'y'

    if awns=='y':
        MASKs = np.zeros(areaFile.attrs['ROI_patches'].shape)
        for i in range(areaFile.attrs['ROI_patches'].shape[2]):
            patch = areaFile.attrs['ROI_patches'][:,:,i]
            
            tt0 = exposure.equalize_hist(patch)
            tt = 255*tt0/np.max(tt0)
            thresh = 1*tt>0.3*255
            thresh2 = 1*tt<0.1*255

            tt[thresh] = 255
            tt[thresh2] = 0



            edges = canny(tt, sigma=2, low_threshold=20, high_threshold=30)
            try_radii = np.arange(3,5)
            res = hough_circle(edges, try_radii)
            ridx, r, c = np.unravel_index(np.argmax(res), res.shape)
            
            r, c, try_radii[ridx]
            image = np.zeros([20,20,4])
            cx, cy = circle_perimeter(c, r, try_radii[ridx]+2)
            
            try:
                image[cy, cx] = (220, 80, 20, 1)
                original_image = np.copy(image[:,:,3])
                chull = convex_hull_image(original_image)
                image[chull,:] = (255, 0, 0, .4)
                
            except:
                pass
            
            
            MASKs[:,:,i] = (1*image[:,:,-1]>0)
        if addMasks:
            areaFile.attrs['ROI_masks'] = MASKs
    else:
        pass    

    return np.array(MASKs)

 def improve_lungs(dist_lungs_hulls, final_area_filter):
     """ Upravi masku lungs pro specialni ucely """
     lungs_hulls = dist_lungs_hulls > 0
     new_filter = np.zeros(final_area_filter.shape)
     for i, area_slice in enumerate(final_area_filter):
         convex = convex_hull_image(area_slice)
         new_filter[i] = lungs_hulls[i] & convex
     return new_filter

def find_convex_hull_rectangle(grey_image, mask=None, canny_threshold=50):
    edge_map = cv2.Canny(grey_image, canny_threshold, 3 * canny_threshold, apertureSize=3)
    if mask is not None:
        edge_map = edge_map*mask
    hull = convex_hull_image(edge_map)
    (_, contours, _) = cv2.findContours(hull.astype(np.uint8), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    assert len(contours) == 1
    rect = cv2.minAreaRect(contours[0])
    max_box = cv2.boxPoints(rect)
    return max_box.astype(np.int)

def test_pathological_qhull_example():
    image = np.array(
                [[0, 0, 0, 0, 1, 0, 0],
                 [0, 0, 1, 1, 1, 1, 1],
                 [1, 1, 1, 0, 0, 0, 0]], dtype=bool)
    expected = np.array(
                [[0, 0, 0, 1, 1, 1, 0],
                 [0, 1, 1, 1, 1, 1, 1],
                 [1, 1, 1, 1, 0, 0, 0]], dtype=bool)
    assert_array_equal(convex_hull_image(image), expected)

def wrapper_regions(bestregions, opening_param = 3, mshape = ((0,1,0),(1,1,1),(0,1,0)) ):

    zdim, xdim, ydim = bestregions.shape

    wregions = np.zeros_like(bestregions)

    for sidx in range(zdim):
        if np.sum(bestregions[sidx]) > 0:
            wregions[sidx] = convex_hull_image(bestregions[sidx])

    return wregions

def fill_poly(poly_y, poly_x, shape):
    bbox = np.zeros((4), dtype=np.int32)
    bbox[0] = np.min(poly_y)
    bbox[1] = np.min(poly_x)
    bbox[2] = np.max(poly_y)
    bbox[3] = np.max(poly_x)
    
    mask = np.zeros(shape, dtype = np.bool_)
    mask[poly_y.astype(np.int), poly_x.astype(np.int)] = True
    mask = morph.convex_hull_image(mask).astype(np.int8)
    return mask, bbox

 def improve_bone_hull(bone_area_filter):
     basic_loc_filter = body & (coronal > 0)
     mask = np.zeros(bone_area_filter.shape)
     
     for i in range(len(bone_area_filter)):
         try:
             down_mask = bone_area_filter[i] & basic_loc_filter[i]
             mask[i] = convex_hull_image(down_mask) | bone_area_filter[i]
             mask[i] = scipy.ndimage.morphology.binary_fill_holes(mask[i])
         except:
             pass
     mask = mask > 0
     return mask

def process_mask(mask):
    convex_mask = np.copy(mask)
    for i_layer in range(convex_mask.shape[0]):
        mask1  = np.ascontiguousarray(mask[i_layer])
        if np.sum(mask1)>0:
            mask2 = convex_hull_image(mask1)
            if np.sum(mask2)>2*np.sum(mask1):
                mask2 = mask1
        else:
            mask2 = mask1
        convex_mask[i_layer] = mask2
    struct = generate_binary_structure(3,1)  
    dilatedMask = binary_dilation(convex_mask,structure=struct,iterations=10) 
    return dilatedMask

def test_basic():
    image = np.array(
        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 0, 1, 0, 0, 0],
         [0, 0, 1, 0, 0, 0, 1, 0, 0],
         [0, 1, 0, 0, 0, 0, 0, 1, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)

    expected = np.array(
        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 1, 1, 0, 0, 0],
         [0, 0, 1, 1, 1, 1, 1, 0, 0],
         [0, 1, 1, 1, 1, 1, 1, 1, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)

    assert_array_equal(convex_hull_image(image), expected)

    # Test that an error is raised on passing a 3D image:
    image3d = np.empty((5, 5, 5))
    with pytest.raises(ValueError):
        convex_hull_image(image3d)

def test_qhull_offset_example():
    nonzeros = (([1367, 1368, 1368, 1368, 1369, 1369, 1369, 1369, 1369, 1370,
                  1370, 1370, 1370, 1370, 1370, 1370, 1371, 1371, 1371, 1371,
                  1371, 1371, 1371, 1371, 1371, 1372, 1372, 1372, 1372, 1372,
                  1372, 1372, 1372, 1372, 1373, 1373, 1373, 1373, 1373, 1373,
                  1373, 1373, 1373, 1374, 1374, 1374, 1374, 1374, 1374, 1374,
                  1375, 1375, 1375, 1375, 1375, 1376, 1376, 1376, 1377]),
                ([151, 150, 151, 152, 149, 150, 151, 152, 153, 148, 149, 150,
                 151, 152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 154,
                 155, 146, 147, 148, 149, 150, 151, 152, 153, 154, 146, 147,
                 148, 149, 150, 151, 152, 153, 154, 147, 148, 149, 150, 151,
                 152, 153, 148, 149, 150, 151, 152, 149, 150, 151, 150]))
    image = np.zeros((1392, 1040), dtype=bool)
    image[nonzeros] = True
    expected = image.copy()
    assert_array_equal(convex_hull_image(image), expected)

def create_mask(frame):
    """"Create a big mask that encompasses all the cells"""
    
    # detect ridges
    ridges = enhance_ridges(frame)

    # threshold ridge image
    thresh = filters.threshold_otsu(ridges)
    thresh_factor = 1.1
    prominent_ridges = ridges > thresh_factor*thresh
    prominent_ridges = morphology.remove_small_objects(prominent_ridges, min_size=128)

    # the mask contains the prominent ridges
    mask = morphology.convex_hull_image(prominent_ridges)
    mask = morphology.binary_erosion(mask, disk(10))
    return mask

def test_basic():
    image = np.array(
        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 0, 1, 0, 0, 0],
         [0, 0, 1, 0, 0, 0, 1, 0, 0],
         [0, 1, 0, 0, 0, 0, 0, 1, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)

    expected = np.array(
        [[0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 1, 1, 0, 0, 0],
         [0, 0, 1, 1, 1, 1, 1, 0, 0],
         [0, 1, 1, 1, 1, 1, 1, 1, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)

    assert_array_equal(convex_hull_image(image), expected)

 def get_bone_hull(data3dr_tmp, body):
     # jeste by bylo dobre u patere vynechat konvexni obal
     bone = (data3dr_tmp>min_bone_thr) & body
     mask = np.zeros(data3dr_tmp.shape)
     
     step = 32
     for i in range(len(data3dr_tmp)):
         mask[i] = np.sum(
             bone[
             int(max(0, i - step / 2)):
             int(min(len(data3dr_tmp), i + step / 2))
             ], axis = 0
         ) >= 1
         try:
             mask[i] = convex_hull_image(mask[i])
         except:
             pass
     
     mask = mask>0
     return mask

def mask_readoutstreaks(image):
    '''logarithmic image to edge detect faint features

    This requires the `scikit-image <http://scikit-image.org/>`_ package.
    '''

    from skimage import filters as skfilter
    from skimage.morphology import convex_hull_image

    logimage = np.log10(np.clip(image, 1, 1e5)) / 5
    # Mask overexposed area + sobel edge detect
    mask = (skfilter.sobel(logimage) > 0.1) | (image > 0.6 * np.max(image))
    # pick out the feature that contain the center
    # I hope that this always is bit enough
    mask, lnum = ndimage.label(mask, structure=np.ones((3, 3), dtype=bool))

    i = mask[ndimage.center_of_mass(image)]
    mask = (mask == i)
    # fill any holes in that region
    return convex_hull_image(mask)