def test_median_default_value():
    a = np.zeros((3, 3), dtype=np.uint8)
    a[1] = 1
    full_selem = np.ones((3, 3), dtype=np.uint8)
    assert_equal(rank.median(a), rank.median(a, full_selem))
    assert rank.median(a)[1, 1] == 0
    assert rank.median(a, disk(1))[1, 1] == 1

def test_percentile_median():
    # check that percentile p0 = 0.5 is identical to local median
    img = data.camera()
    img16 = img.astype(np.uint16)
    selem = disk(15)
    # check for 8bit
    img_p0 = rank.percentile(img, selem=selem, p0=.5)
    img_max = rank.median(img, selem=selem)
    assert_equal(img_p0, img_max)
    # check for 16bit
    img_p0 = rank.percentile(img16, selem=selem, p0=.5)
    img_max = rank.median(img16, selem=selem)
    assert_equal(img_p0, img_max)

def watershed(image):
    hsv_image = color.rgb2hsv(image)

    low_res_image = rescale(hsv_image[:, :, 0], SCALE)
    local_mean = mean(low_res_image, disk(50))
    local_minimum_flat = np.argmin(local_mean)
    local_minimum = np.multiply(np.unravel_index(local_minimum_flat, low_res_image.shape), round(1 / SCALE))

    certain_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_bone_pixels[
    local_minimum[0] - INITIAL_WINDOW_SIZE/2:local_minimum[0]+INITIAL_WINDOW_SIZE/2,
    local_minimum[1] - INITIAL_WINDOW_SIZE/2:local_minimum[1]+INITIAL_WINDOW_SIZE/2
    ] = True

    certain_non_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_non_bone_pixels[0:BORDER_SIZE, :] = True
    certain_non_bone_pixels[-BORDER_SIZE:-1, :] = True
    certain_non_bone_pixels[:, 0:BORDER_SIZE] = True
    certain_non_bone_pixels[:, -BORDER_SIZE:-1] = True

    smoothed_hsv = median(hsv_image[:, :, 0], disk(50))
    threshold = MU * np.median(smoothed_hsv[certain_bone_pixels])

    possible_bones = np.zeros_like(hsv_image[:, :, 0])
    possible_bones[smoothed_hsv < threshold] = 1

    markers = np.zeros_like(possible_bones)
    markers[certain_bone_pixels] = 1
    markers[certain_non_bone_pixels] = 2

    labels = morphology.watershed(-possible_bones, markers)

    return labels

def image_prep(img):
	# apply filters for better contrast and noise removal
	img_gamma = adjust_gamma(img, 0.7)
	img_median = median(img_gamma, disk(1))

	# apply threshold
	val =  threshold_otsu(img_median)
	img_otsu = img_median > val

	# label image regions
	label_image = label(img_otsu)

	candidates = []
	for region in regionprops(label_image):
		minr, minc, maxr, maxc = region.bbox
		if (maxr - minr > maxc - minc):
			candidates.append(region)

	# find numbers
	areas = []
	for candidate in candidates:
		areas.append(candidate.area)
	areas.sort()
	areas.reverse()

	n = 1
	v = []
	for candidate in candidates:
		if (candidate.area == areas[0] or candidate.area == areas[1] or candidate.area == areas[2]):
			v.append(candidate.image)
			imsave('num%d.png' % n, candidate.image)
			n += 1
	return v

def _repair_bad_pixels_generic(raw, coords, method='median'):
    if median is None:
        raise RuntimeError('scikit-image is required for repair_bad_pixels if the Bayer pattern is not 2x2')
    
    color_masks = _colormasks(raw)
        
    rawimg = raw.raw_image_visible
    r = 5
    kernel = np.ones((r,r))
    for color_mask in color_masks:       
        mask = np.zeros_like(color_mask)
        mask[coords[:,0],coords[:,1]] = True
        mask &= color_mask
        
        # interpolate all bad pixels belonging to this color
        if method == 'mean':
            # With mean filtering we have to ignore the bad pixels as they
            # would influence the mean too much.
            # FIXME could lead to invalid values if bad pixels are clustered
            #       such that no valid pixels are left in a block
            raise NotImplementedError
        elif method == 'median':
            # bad pixels won't influence the median in most cases and just using
            # the color mask prevents bad pixel clusters from producing
            # bad interpolated values (NaNs)
            smooth = median(rawimg, kernel, mask=color_mask)
        else:
            raise ValueError
        
        rawimg[mask] = smooth[mask]   

def _find_bad_pixel_candidates_generic(raw, isCandidateFn):
    if median is None:
        raise RuntimeError('scikit-image is required if the Bayer pattern is not 2x2')
    
    color_masks = _colormasks(raw)   
    rawimg = raw.raw_image_visible    
    coords = []
    
    # median filtering for each channel    
    r = 5
    kernel = np.ones((r,r))
    for mask in color_masks:
        # skimage's median is quite slow, it uses an O(r) filtering algorithm.
        # There exist O(log(r)) and O(1) algorithms, see https://nomis80.org/ctmf.pdf.
        # Also, we only need the median values for the masked pixels.
        # Currently, they are calculated for all pixels for each color.
        med = median(rawimg, kernel, mask=mask)
        
        # detect possible bad pixels
        candidates = isCandidateFn(rawimg, med)
        candidates &= mask
        
        y,x = np.nonzero(candidates)
        # note: the following is much faster than np.transpose((y,x))
        candidates = np.empty((len(y),2), dtype=y.dtype)
        candidates[:,0] = y
        candidates[:,1] = x
        
        coords.append(candidates)
        
    return coords

def preproc_2d_img(img):
    #silence Possible precision loss when converting from float64 to uint8 warning
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        picture = np.rot90(img)
        picture /= np.max(picture)
        # apply median filter to make stripes more vissible
        picture = median(picture, disk(1))
        return picture

    def smooth(self):
        # TODO: there is non nan in the ff img, or?
        mask = self.flatField == 0
        from skimage.filters.rank import median, mean
        from skimage.morphology import disk

        ff = mean(median(self.flatField, disk(5), mask=~mask),
                  disk(13), mask=~mask)

        return ff.astype(float) / ff.max(), mask

def label_from_thresh(frame, thresh, parameters):

    smooth = parameters['smooth']
    min_distance = np.int(parameters['nuc_distance'])
    image = rank.median(frame, disk(smooth))
    image = rank.enhance_contrast(image, disk(smooth))
    im_max = image.max()
    if im_max < thresh:
        return np.zeros(image.shape, dtype=np.int32)
    else:
        image = image > thresh
        distance = ndimage.distance_transform_edt(image)
        local_maxi = peak_local_max(distance,
                                    footprint=disk(min_distance),
                                    #min_distance=min_distance,
                                    indices=False, labels=image)
        markers = ndimage.label(local_maxi)[0]
        return watershed(-distance, markers, mask=image)

def compactness(_img):

    ## Lookup table for the empirical distribution of median(img)/img
    ## in the case of a random image (white noise), for varying
    ## proportions of white pixels in the image (0.1, 0.2, ..., 1.0).
    ## The distributions are approximated by Gaussians, and the
    ## corresponding means and standard deviations are stored.

    prop = np.linspace(0.1, 1.0, 10)
    emp_distrib_mean = np.array([4.4216484763437432e-06, 0.0011018116582350559,
                                 0.042247116747488218, 0.34893587605251208,
                                 1.0046008733628913, 1.4397675817057451,
                                 1.4115741958770296, 1.2497935146232551,
                                 1.1111058415275834, 1.0])
    emp_distrib_std = np.array([2.7360459073474441e-05, 0.00051125394394966434,
                                0.0038856377648490894, 0.012029872915543046,
                                0.013957075037020938, 0.0057246251730834283,
                                0.0028750796874699143, 0.0023709207886137384,
                                0.0015018959493632007, 0.0])
    if _img.ndim != 2:
        raise ValueError('The input image must be a 2D binary image.')

    _img = (_img != 0).astype(np.uint8)
    _med = median(_img, disk(3))

    # "proportion of white pixels" in the image:
    swp = np.sum(_img, dtype=np.float64)
    pwp = swp / _img.size

    # compactness coefficient
    cf = np.sum(_med, dtype=np.float64) / swp

    # standardize using the "closest" Gaussian from the list of empirical
    # distributions:
    k = np.argmin(np.abs(prop - pwp))

    # this should make the coeff more or less normally distributed N(0,1)
    cf = (cf - emp_distrib_mean[k]) / emp_distrib_std[k]

    return cf

def test_line_profile_dynamic():
    """Test a line profile updating after an image transform"""
    image = data.coins()[:-50, :]  # shave some off to make the line lower
    image = skimage.img_as_float(image)
    viewer = ImageViewer(image)

    lp = LineProfile(limits='dtype')
    viewer += lp

    line = lp.get_profiles()[-1][0]
    assert line.size == 129
    assert_almost_equal(np.std(viewer.image), 0.208, 3)
    assert_almost_equal(np.std(line), 0.229, 3)
    assert_almost_equal(np.max(line) - np.min(line), 0.725, 1)

    viewer.image = skimage.img_as_float(median(image,
                                               selem=disk(radius=3)))

    line = lp.get_profiles()[-1][0]
    assert_almost_equal(np.std(viewer.image), 0.198, 3)
    assert_almost_equal(np.std(line), 0.220, 3)
    assert_almost_equal(np.max(line) - np.min(line), 0.639, 1)

def check_all():
    np.random.seed(0)
    image = np.random.rand(25, 25)
    selem = morphology.disk(1)
    refs = np.load(os.path.join(skimage.data_dir, "rank_filter_tests.npz"))

    assert_equal(refs["autolevel"], rank.autolevel(image, selem))
    assert_equal(refs["autolevel_percentile"], rank.autolevel_percentile(image, selem))
    assert_equal(refs["bottomhat"], rank.bottomhat(image, selem))
    assert_equal(refs["equalize"], rank.equalize(image, selem))
    assert_equal(refs["gradient"], rank.gradient(image, selem))
    assert_equal(refs["gradient_percentile"], rank.gradient_percentile(image, selem))
    assert_equal(refs["maximum"], rank.maximum(image, selem))
    assert_equal(refs["mean"], rank.mean(image, selem))
    assert_equal(refs["mean_percentile"], rank.mean_percentile(image, selem))
    assert_equal(refs["mean_bilateral"], rank.mean_bilateral(image, selem))
    assert_equal(refs["subtract_mean"], rank.subtract_mean(image, selem))
    assert_equal(refs["subtract_mean_percentile"], rank.subtract_mean_percentile(image, selem))
    assert_equal(refs["median"], rank.median(image, selem))
    assert_equal(refs["minimum"], rank.minimum(image, selem))
    assert_equal(refs["modal"], rank.modal(image, selem))
    assert_equal(refs["enhance_contrast"], rank.enhance_contrast(image, selem))
    assert_equal(refs["enhance_contrast_percentile"], rank.enhance_contrast_percentile(image, selem))
    assert_equal(refs["pop"], rank.pop(image, selem))
    assert_equal(refs["pop_percentile"], rank.pop_percentile(image, selem))
    assert_equal(refs["pop_bilateral"], rank.pop_bilateral(image, selem))
    assert_equal(refs["sum"], rank.sum(image, selem))
    assert_equal(refs["sum_bilateral"], rank.sum_bilateral(image, selem))
    assert_equal(refs["sum_percentile"], rank.sum_percentile(image, selem))
    assert_equal(refs["threshold"], rank.threshold(image, selem))
    assert_equal(refs["threshold_percentile"], rank.threshold_percentile(image, selem))
    assert_equal(refs["tophat"], rank.tophat(image, selem))
    assert_equal(refs["noise_filter"], rank.noise_filter(image, selem))
    assert_equal(refs["entropy"], rank.entropy(image, selem))
    assert_equal(refs["otsu"], rank.otsu(image, selem))
    assert_equal(refs["percentile"], rank.percentile(image, selem))
    assert_equal(refs["windowed_histogram"], rank.windowed_histogram(image, selem))

def generate():
    np.random.seed(None)
    ohe = OneHotEncoder(sparse=False)

    hmap = np.array(DS.diamond_square((200, 200), -1, 1, 0.35))
    + np.array(DS.diamond_square((200, 200), -1, 1, 0.55))
    + np.array(DS.diamond_square((200, 200), -1, 1, 0.75))

    hmap_flatten = np.array(hmap).ravel()[:, None]
    kmeans = KMeans(n_clusters=n_colors, random_state=0).fit(hmap_flatten)
    labels_hmap = kmeans.predict(hmap_flatten)[:, None]

    # Back to rectangular
    labels_hmap = labels_hmap.reshape([hmap.shape[0], hmap.shape[1]])
    labels_hmap = median(labels_hmap.astype(np.uint8), disk(5))
    labels_hmap = resize(labels_hmap, dims, order=0, preserve_range=True)

    labels_hmap = ohe.fit_transform(labels_hmap.ravel()[:, None])

    # Reshape
    hmap_masks = labels_hmap.reshape([dims[0], dims[1], n_colors])
    hmap_masks = hmap_masks.transpose([2, 0, 1])

    return hmap_masks

    def background_subtraction(self, img, method='avg'):
        #width, height = img.shape
        
        if method=='avg': 
            # vigra
            #kernel = vigra.filters.averagingKernel(radius)
            #bgsub = img - vigra.filters.convolve(self.ut.to_float(img), kernel)
            
            # with skimage           
            se = disk(self.settings.background_subtraction['radius'])
            bgsub = img.astype(np.dtype('float')) - rank.mean(img, se)
            bgsub[bgsub < 0] = 0
            bgsub = bgsub.astype(img.dtype)

        if method=='med': 
            # vigra
            #kernel = vigra.filters.averagingKernel(radius)
            #bgsub = img - vigra.filters.convolve(self.ut.to_float(img), kernel)
            
            # with skimage           
            se = disk(self.settings.background_subtraction['radius'])
            bgsub = img.astype(np.dtype('float')) - rank.median(img, se)
            bgsub[bgsub < 0] = 0
            bgsub = bgsub.astype(img.dtype)
            
        elif method=='constant_median':
            # vigra
            #bgsub = img - np.median(np.array(img))
            
            # with skimage            
            bgsub = img - np.median(img)

            bgsub[bgsub < 0] = 0
            bgsub = bgsub.astype(img.dtype)
            
        return bgsub

def cr_med(image, selem):
    return median(image=image, selem=selem)

def median_filter(image, radius):
    return median(image, selem=disk(radius))

from skimage.morphology import disk

noise = np.random.random(noisy_image.shape)
noisy_image = img_as_ubyte(data.camera())
noisy_image[noise > 0.99] = 255
noisy_image[noise < 0.01] = 0

fig, ax = plt.subplots(2, 2, figsize=(10, 7), sharex=True, sharey=True)
ax1, ax2, ax3, ax4 = ax.ravel()

ax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)
ax1.set_title('Noisy image')
ax1.axis('off')
ax1.set_adjustable('box-forced')

ax2.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax2.set_title('Median $r=1$')
ax2.axis('off')
ax2.set_adjustable('box-forced')


ax3.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax3.set_title('Median $r=5$')
ax3.axis('off')
ax3.set_adjustable('box-forced')


ax4.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax4.set_title('Median $r=20$')
ax4.axis('off')
ax4.set_adjustable('box-forced')

def median_each(image, selem):
    return median(image, selem)

def median_hsv(image, selem):
    return median(image, selem)

 def transform(self, images):
     from skimage.filters.rank import median
     from skimage.morphology import disk
     return np.array([median(image, disk(self._radius))
                      for image in images])