def test_empty_selem():
    # check that min, max and mean returns zeros if structuring element is
    # empty

    image = np.zeros((5, 5), dtype=np.uint16)
    out = np.zeros_like(image)
    mask = np.ones_like(image, dtype=np.uint8)
    res = np.zeros_like(image)
    image[2, 2] = 255
    image[2, 3] = 128
    image[1, 2] = 16

    elem = np.array([[0, 0, 0], [0, 0, 0]], dtype=np.uint8)

    rank.mean(image=image, selem=elem, out=out, mask=mask,
              shift_x=0, shift_y=0)
    assert_equal(res, out)
    rank.geometric_mean(image=image, selem=elem, out=out, mask=mask,
                        shift_x=0, shift_y=0)
    assert_equal(res, out)
    rank.minimum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_equal(res, out)
    rank.maximum(image=image, selem=elem, out=out, mask=mask,
                 shift_x=0, shift_y=0)
    assert_equal(res, out)

def filter_img(img):
    selem = square(11)
    img[:, :, 0] = rank.mean(img[:, :, 0], selem=selem)
    img[:, :, 1] = rank.mean(img[:, :, 1], selem=selem)
    img[:, :, 2] = rank.mean(img[:, :, 2], selem=selem)
    #return np.array(img, dtype=float)
    return img_as_float(img)

    def test_inplace_output(self):
        # rank filters are not supposed to filter inplace

        selem = disk(20)
        image = (np.random.rand(500, 500) * 256).astype(np.uint8)
        out = image
        with pytest.raises(NotImplementedError):
            rank.mean(image, selem, out=out)

def test_random_sizes():
    # make sure the size is not a problem

    niter = 10
    elem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=np.uint8)
    for m, n in np.random.random_integers(1, 100, size=(10, 2)):
        mask = np.ones((m, n), dtype=np.uint8)

        image8 = np.ones((m, n), dtype=np.uint8)
        out8 = np.empty_like(image8)
        rank.mean(image=image8, selem=elem, mask=mask, out=out8,
                  shift_x=0, shift_y=0)
        assert_equal(image8.shape, out8.shape)
        rank.mean(image=image8, selem=elem, mask=mask, out=out8,
                  shift_x=+1, shift_y=+1)
        assert_equal(image8.shape, out8.shape)

        image16 = np.ones((m, n), dtype=np.uint16)
        out16 = np.empty_like(image8, dtype=np.uint16)
        rank.mean(image=image16, selem=elem, mask=mask, out=out16,
                  shift_x=0, shift_y=0)
        assert_equal(image16.shape, out16.shape)
        rank.mean(image=image16, selem=elem, mask=mask, out=out16,
                  shift_x=+1, shift_y=+1)
        assert_equal(image16.shape, out16.shape)

        rank.mean_percentile(image=image16, mask=mask, out=out16,
                             selem=elem, shift_x=0, shift_y=0, p0=.1, p1=.9)
        assert_equal(image16.shape, out16.shape)
        rank.mean_percentile(image=image16, mask=mask, out=out16,
                             selem=elem, shift_x=+1, shift_y=+1, p0=.1, p1=.9)
        assert_equal(image16.shape, out16.shape)

def test_selem_dtypes():

    image = np.zeros((5, 5), dtype=np.uint8)
    out = np.zeros_like(image)
    mask = np.ones_like(image, dtype=np.uint8)
    image[2, 2] = 255
    image[2, 3] = 128
    image[1, 2] = 16

    for dtype in (np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64):
        elem = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=dtype)
        rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0)
        assert_equal(image, out)
        rank.mean_percentile(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0)
        assert_equal(image, out)

def zoom(imName, img, imageclass, conf):
	imName = imName+1
	im = Image.open(img)
	if not isdir(conf.output_folder+imageclass):
		try:
			mkdir(conf.output_folder+imageclass)
		except:
			pass
	fileMatch = 0
	for file in listdir(conf.output_folder+"/"+imageclass):
		if fnmatch.fnmatch(file, str(imName)+"_zoom_"+'*.jpg'):
			fileMatch = fileMatch+1
			if fileMatch>=3:
				if conf.VERBOSE: print "exists, so breaking: #"+str(imName)+" in "+str(imageclass)
				return str(imName)
	x, y = im.size
	x1=0
	y1=0
	means=[]
	while y1<=y-480:
		while x1<=x-640:
			ims = im.crop((x1, y1, x1+640, y1+480))
			mean1 = mean(np.array(ims)[:,:,1], disk(700))
			means.append(((x1,y1),int(mean1[0][0])))
			x1 = x1+160
		x1=0
		y1 = y1+120
	zoomNRC(imName, img, imageclass, conf, conf.mean_threshold, means, im)
	return str(imName)

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 substract_mean(self, radius_disk):
     circle = disk(radius_disk)
     for i_rot in np.arange(self.stack_height):
         stack_slice = self.score_stack[:,:,i_rot]
         norm_factor = max(stack_slice.min(), stack_slice.max(), key=abs)
         stack_slice *= 1./norm_factor 
         stack_slice = np.array(rank.mean(stack_slice, selem=circle), dtype=np.float32)
         stack_slice *= norm_factor
         self.score_stack[:,:,i_rot] = self.score_stack[:,:,i_rot]-stack_slice

def test_smallest_selem16():
    # check that min, max and mean returns identity if structuring element
    # contains only central pixel

    image = np.zeros((5, 5), dtype=np.uint16)
    out = np.zeros_like(image)
    mask = np.ones_like(image, dtype=np.uint8)
    image[2, 2] = 255
    image[2, 3] = 128
    image[1, 2] = 16

    elem = np.array([[1]], dtype=np.uint8)
    rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0)
    assert_equal(image, out)
    rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0)
    assert_equal(image, out)
    rank.maximum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0)
    assert_equal(image, out)

def skmean(image):

    from skimage.filters.rank import mean

    mean_filtered = mean(image, disk(30))

    print mean_filtered.min(), mean_filtered.max()

    return mean_filtered

def test_16bit():
    image = np.zeros((21, 21), dtype=np.uint16)
    selem = np.ones((3, 3), dtype=np.uint8)

    for bitdepth in range(17):
        value = 2 ** bitdepth - 1
        image[10, 10] = value
        assert rank.minimum(image, selem)[10, 10] == 0
        assert rank.maximum(image, selem)[10, 10] == value
        assert rank.mean(image, selem)[10, 10] == int(value / selem.size)

    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 _coarsenImage(image, f):
    '''
    seems to be a more precise (but slower)
    way to down-scale an image
    '''
    from skimage.morphology import square
    from skimage.filters import rank
    from skimage.transform._warps import rescale
    selem = square(f)
    arri = rank.mean(image, selem=selem)
    return rescale(arri, 1 / f, order=0)

 def _segment_edge_areas(self, edges, disk_size, mean_threshold, min_object_size):
     """
     Takes a greyscale image (with brighter colors corresponding to edges) and returns a
     binary image where white indicates an area with high edge density and black indicates low density.
     """
     # Convert the greyscale edge information into black and white (ie binary) image
     threshold = threshold_otsu(edges)
     # Filter out the edge data below the threshold, effectively removing some noise
     raw_channel_areas = edges <= threshold
     # Smooth out the data
     channel_areas = rank.mean(raw_channel_areas, disk(disk_size)) < mean_threshold
     # Remove specks and blobs that are the result of artifacts
     clean_channel_areas = remove_small_objects(channel_areas, min_size=min_object_size)
     # Fill in any areas that are completely surrounded by the areas (hopefully) covering the channels
     return ndimage.binary_fill_holes(clean_channel_areas)

def test_16bit():
    image = np.zeros((21, 21), dtype=np.uint16)
    selem = np.ones((3, 3), dtype=np.uint8)

    for bitdepth in range(17):
        value = 2 ** bitdepth - 1
        image[10, 10] = value
        if bitdepth > 11:
            expected = ['Bitdepth of %s' % (bitdepth - 1)]
        else:
            expected = []
        with expected_warnings(expected):
            assert rank.minimum(image, selem)[10, 10] == 0
            assert rank.maximum(image, selem)[10, 10] == value
            assert rank.mean(image, selem)[10, 10] == int(value / selem.size)

    def color_adjustment(self, img, mask=None, gaussian_std=.0, gamma=1.0, contrast = 1.0, brightness = 0, mult_rgb = np.array([1.0, 1.0, 1.0]), blur_radius = 0):
        img **= gamma
        img *= contrast
        img += np.random.randn(*img.shape).astype('float32') * gaussian_std
        img += brightness
        img *= mult_rgb
        np.clip(img, 0.0, 1.0, img)
        blur_mask = None
        if mask is not None:
            blur_mask = random.choice([mask,1-mask,np.ones_like(mask)])
        if blur_radius > 0:
            selem = disk(blur_radius)
            tmp_img = img.copy()
            for i in range(img.shape[2]):
                img[:, :, i] = rank.mean(img[:, :, i], selem=selem,mask=blur_mask) / 255.0
            img[np.where(blur_mask == 0)] = tmp_img[np.where(blur_mask==0)]

        return img

    def meanFilter(self, m=3, array=numpy.empty(0)):
        """
            Mean filtering, replaces the intensity value, by the average 
            intensity of a pixels neighbours including itself. 
            m is the size of the filter, default is 3x3

            @method meanFilter
            @param m {int} The width and height of the m x m filtering matrix, 
                default is 3.
            @param array {numpy array} the array which the operation is carried 
                out on.
        """
        self.__printStatus("Mean filtering " + str(m) + "x" + str(m) + "...")
        if not array.any():
            array = self.image_array
        if array.dtype not in ["uint8", "uint16"]:
            array = numpy.uint8(array)
        mean3x3filter = rank.mean(array, square(m), mask=self.mask)
        self.image_array = mean3x3filter * self.mask
        self.__printStatus("[done]", True)
        return self

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 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

unique = np.unique(R)
centroid = ndimage.measurements.center_of_mass(gel, R, unique)
# finding the intensity and priting results:
for i in range(len(centroid)):
	print 'Spot %d x: %f y: %f intensity: %d' % (i+1, centroid[i][0], centroid[i][1], gel[int(centroid[i][0]), int(centroid[i][1])])

colored = label2rgb(R, gel, bg_label=0)
for i in range(len(centroid)):
	rr, cc = circle(centroid[i][0], centroid[i][1], 2)
	colored[rr, cc] = (1,0,0)

io.imsave(argv[2], colored)

print "Smoothing the image before running watershed..."
loc_mean = mean(gel, disk(1))
smooth_M = watershed(loc_mean, ' ')
print "Number of spots found on the smoothed image without any post process is:"
print len(np.unique(smooth_M))
R = morphology.remove_small_objects(smooth_M.astype(int), 6)
colored = label2rgb(R, gel, bg_label=0)
io.imsave('smooth_cleaned.png', colored)
print "Number of spots found on the smoothed image after image processing is:"
print len(np.unique(R))

#############################
#Answering the question 1, using different edge operator methods to detect markers:
print 'The number of spots found using roberts gradient method is:'
M_roberts = watershed(gel, 'roberts')
R_roberts = morphology.remove_small_objects(M_roberts.astype(int), 64)
print len(np.unique(R_roberts))