def focus_score(self): 
     f_score = (color.rgb2grey(self.img) - erosion(color.rgb2grey(self.img), square(4)))
     non_zero_pixel_area = self.get_nonzero_pixel_area(f_score)
     #print("focus score: " + str(np.sum(f_score) / non_zero_pixel_area))
     #plt.imshow(f_score)
     #plt.show()
     return np.sum(f_score) / non_zero_pixel_area 

def compare_images(imageA, imageB, title, show_plot=True):
    """
    computes the mean squared error and structural similarity
    """

    # index values for mean squared error
    if VERBOSE: print("comparing mean squared error...")
    m = mse(imageA, imageB)

    # convert the images to grayscale
    if VERBOSE: print("converting to greyscale...")
    imageA_grey = rgb2grey(imageA)
    imageB_grey = rgb2grey(imageB)

    # uses image copies to avoid runtime warning for ssim computation
    img1_grey = np.copy(imageA_grey)
    img2_grey = np.copy(imageB_grey)

    # index values for structural similarity
    if VERBOSE: print("comparing structural similarity...")
    s = ssim(img1_grey, img2_grey)

    if show_plot:
        if VERBOSE: print("plotting images...")
        try:
            import matplotlib.pyplot as plt
        except:
            print("Error importing pyplot from matplotlib, please install matplotlib package first...")
            sys.tracebacklimit=0
            raise Exception("Importing matplotlib failed")

        # setup the figure
        fig, ax = plt.subplots(2, 2)
        fig.suptitle("%s\nMSE: %.5f, SSIM: %.5f" % (title, m, s))

        ax[0][0].text(-10, -10, 'MSE: %.5f' %(m))

        # show first image
        ax[0][0].imshow(imageA, cmap=plt.cm.gray)
        ax[0][0].axis("off")

        # show the second image
        ax[0][1].imshow(imageB, cmap=plt.cm.gray)
        ax[0][1].axis("off")

        ax[1][0].text(-10, -10, 'SSIM: %.5f' %(s))

        # show first grey image
        ax[1][0].imshow(img1_grey, cmap=plt.cm.gray)
        ax[1][0].axis("off")

        # show the second grey image
        ax[1][1].imshow(img2_grey, cmap=plt.cm.gray)
        ax[1][1].axis("off")

        # show the images
        plt.show()

    return m, s

def read_image(filename):
    """Read an image from the disk and output data arrays."""
    image_array_rgb = misc.imread(filename, mode='RGB')
    #  image_array_grey = misc.imread(filename, flatten=True, mode='F')
    image_array_grey = color.rgb2grey(image_array_rgb)*255
    image_array_luv = color.rgb2luv(image_array_rgb)
    return image_array_rgb, image_array_grey, image_array_luv

def featurize(img_name):
    """Load an image and convert it into a dictionary of features"""
    img = plt.imread(os.path.join('stimuli', img_name + '.png'))
    height, width, _ = img.shape
    features = defaultdict(int)
    for y in range(height):
        for x in range(width):
            features['red'] += img[y][x][0]
            features['green'] += img[y][x][1]
            features['blue'] += img[y][x][2]
            features['alpha'] += img[y][x][3]

    grey = color.rgb2grey(img)
    for y in range(height):
        for x in range(width):
            for key, value in per_pixel(grey, y, x):
                features[key] += value

    # Normalize over image size
    for key, value in features.items():
        features[key] = float(value) / height / width

    features['blob'] = feature.blob_dog(grey).shape[0]
    features['corners'] = feature.corner_peaks(
        feature.corner_harris(grey)).shape[0]
    return features

def feature_extraction(raw_data):
    image = color.rgb2grey(raw_data)

    fd, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16),
                        cells_per_block=(1, 1), visualise=True)

    return hog_image

def modify(img):
    """Randomly modify an image
    
    This is a preprocessing step for training an OCR classifier. It takes
    in an image and casts it to greyscale, reshapes it, and adds some
    (1) rotations, (2) translations and (3) noise.
    
    If more efficiency is needed, we could factor out some of the initial
    nonrandom transforms.
    """
    
    block_size = np.random.uniform(20, 40)
    rotation = 5*np.random.randn()
    
    #print 'BLOCK SIZE', block_size
    #print 'ROTATION  ', rotation
    
    img = color.rgb2grey(img)
    img = transform.resize(img, output_shape=(50,30))
    img = filter.threshold_adaptive(img, block_size=block_size)
    
    # rotate the image
    img = np.logical_not(transform.rotate(np.logical_not(img), rotation))
    # translate the image
    img = shift(img)
    # add some noise to the image
    img = noise(img)
    
    img = transform.resize(img, output_shape=(25,15))
    return filter.threshold_adaptive(img, block_size=25)

def dhash(picture):
    "Compute dhash as uint64."
    img = rgb2grey(resize(picture, (9, 8)))
    h = np.zeros([8], dtype=np.uint8)
    for a in range(8):
        h[a] = TWOS[img[a] > img[a + 1]].sum()
    return (BIGS * h).sum()

def register_feature_calculators():
    return [
        lambda img: GaborFilter.compute_feats(rgb2grey(img), GaborFilter.generate_kernels(2)),
        # lambda img: GLCM.compute_feats(rgb2grey(img), [1, 5, 10, 20], [0, np.pi / 4, np.pi / 2, np.pi * 3 / 4]),
        lambda img: ColorAnalyzer.compute_feats(img, 150, 255, ColorAnalyzer.ColorChannel.Green),
        lambda img: ColorAnalyzer.compute_feats(img, 50, 150, ColorAnalyzer.ColorChannel.Hue),
    ]

def _compute_auto_correlation(image, sigma):
    """Compute auto-correlation matrix using sum of squared differences.

    Parameters
    ----------
    image : ndarray
        Input image.
    sigma : float
        Standard deviation used for the Gaussian kernel, which is used as
        weighting function for the auto-correlation matrix.

    Returns
    -------
    Axx : ndarray
        Element of the auto-correlation matrix for each pixel in input image.
    Axy : ndarray
        Element of the auto-correlation matrix for each pixel in input image.
    Ayy : ndarray
        Element of the auto-correlation matrix for each pixel in input image.

    """

    if image.ndim == 3:
        image = img_as_float(rgb2grey(image))

    imx, imy = _compute_derivatives(image)

    # structure tensore
    Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0)
    Axy = ndimage.gaussian_filter(imx * imy, sigma, mode='constant', cval=0)
    Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode='constant', cval=0)

    return Axx, Axy, Ayy

def processOneImage(inputPath, outputPath):
    image = io.imread(inputPath)
    greyImage = rgb2grey(image)
    threshold = threshold_otsu(greyImage)
    imgout = closing(greyImage > threshold, square(1))
    imgout = crop(imgout)
    imgout = transform.resize(imgout, (max(imgout.shape), max(imgout.shape)))
    io.imsave(outputPath, imgout)

def analyze_image(args):
    """Analyze all wells from all trays in one image."""
    filename, config = args
    LOGGER.debug(filename)
    rows = config["rows"]
    columns = config["columns"]
    well_names = config["well_names"]

    name = splitext(basename(filename))[0]
    if config["parse_dates"]:
        try:
            index = convert_to_datetime(fix_date(name))
        except ValueError as err:
            return {"error": str(err), "filename": filename}
    else:
        index = name

    try:
        image = rgb2grey(imread(filename))
    except OSError as err:
        return {"error": str(err), "filename": filename}

    plate_images = cut_image(image)

    data = dict()

    for i, plate_name in zip(config["plate_indexes"], config["plate_names"]):
        plate = data[plate_name] = dict()
        plate[config["index_name"]] = index
        plate_image = plate_images[i]
        if i // 3 == 0:
            calibration_plate = config["left_image"]
            positions = config["left_positions"]
        else:
            calibration_plate = config["right_image"]
            positions = config["right_positions"]

        try:
            edge_image = canny(plate_image, CANNY_SIGMA)
            offset = align_plates(edge_image, calibration_plate)

            # Add the offset to get the well centers in the analyzed plate.
            well_centers = generate_well_centers(
                np.array(positions) + offset, config["plate_size"], rows,
                columns)
            assert len(well_centers) == rows * columns
            # Add a minimal value to avoid zero division.
            plate_image /= (1 - plate_image + float_info.epsilon)

            well_intensities = [find_well_intensity(plate_image, center)
                                for center in well_centers]

            for well, intensity in zip(well_names, well_intensities):
                plate[well] = intensity
        except (AttributeError, IndexError) as err:
            return {"error": str(err), "filename": filename}

    return data

def image():
    """Load a single image from p1 brain directory.
    output: a single image as a numpy array"""


    inputDir = '{}'.format(all.__path__[0])
    img = load_image('p1-D3-01b.jpg',inputDir)
    img = color.rgb2grey(img)
    return img 

def apply_watermark(filename):
    img = io.imread(filename) # Image in gray scale
    img = color.rgb2grey(img)
    if img.dtype.name != 'uint8':
        img = img * 255
        img = img.astype(numpy.uint8)
    image = img.copy()
    blocks = []
    width, height = image.shape
    hor_block = width / 4
    ver_block = height / 4
    block_counter = 0
    for x in range(0, hor_block):
        for y in range(0, ver_block):
            x_coor = x * 4
            y_coor = y * 4
            block = image[x_coor: x_coor + 4, y_coor: y_coor + 4]
            blocks.append(block)
            block_counter += 1
    n = block_counter
    k = Functions.get_biggest_prime(n)

    for index in range(0, n):
        block_B = blocks[index]
        block_A = (blocks[Functions.mapping(index + 1, k, n) - 1]).copy()
        for x in range(0, 4):
            for y in range(0, 4):
                block_B[x, y] = Functions.removeLSB(block_B[x, y])
        avg_B = Functions.average(block_B)
        for i in range(0, 2):
            for j in range(0, 2):
                i_coor = i * 2
                j_coor = j * 2
                blockBS = block_B[i_coor: i_coor+2, j_coor: j_coor+2]
                average = Functions.average(blockBS)
                v = 0
                if average >= avg_B:
                    v = 1
                p = 1
                if Functions.ones_in_sixMSB(average) % 2 == 0:
                    p = 0
                subblock_a = block_A[i_coor: i_coor+2, j_coor: j_coor+2].copy()
                avg_as = Functions.average(subblock_a)
                r = Functions.split_binary_sixMSB(avg_as)
                if v == 1:
                    v = 2
                if p == 1:
                    p = 2
                if r[2] == 1:
                    r[2] = 2
                if r[4] == 1:
                    r[4] = 2
                blockBS[0][0] = (blockBS[0][0] + v + r[0])
                blockBS[0][1] = (blockBS[0][1] + p + r[1])
                blockBS[1][0] = (blockBS[1][0] + r[2] + r[3])
                blockBS[1][1] = (blockBS[1][1] + r[4] + r[5])
    return image

def normalize(image, subtractMin):
    """
    @params { array-like } image Skimage type acceptable
    @return { np.ndarray }
    """
    if subtractMin:
        return color.rgb2grey(image)
    else:
        return np.divide(a, np.max(a))

def stainspace_to_2d_array(ihc_xyz, channel):
    #rescale = rescale_intensity(ihc_xyz[:, :, channel], out_range=(0,1))
    #stain_array = np.dstack((np.zeros_like(rescale), rescale, rescale))

    #try to not reverse engineer rescale right now
    stain_array = ihc_xyz[:, :, channel]
    #plt.imshow(stain_array)
    gray_array = rgb2grey(stain_array)
    #plt.imshow(gray_array)
    return gray_array

def load_th_image(image_name):
    image = data.load(image_name)
    grey = rgb2grey(image)

    th = threshold_otsu(grey)

    bin = grey >= th

    bin[50:-50,50:-50] = 0
    return bin

def predict(image_url, fileName=None):
    svmObject = getSVMInstance()
    image_url = image_url.replace(' ', '')[:-1] # to remove \n
    fileName = "temp.png"
    getImageAndSave(image_url)
    image = imread(fileName)
    image = skimage.transform.resize(image, (100,100))
    image = rgb2grey(image)
    fd = hog(image, orientations=9, pixels_per_cell=(16, 16), cells_per_block=(1, 1), visualise=False,normalise=True)
    prediction = svmObject.predict(fd)
    return prediction

def preprocess(image, height=50, block_size=50):
    """Turn to greyscale, scale to a height, and then threshold to binary
    """

    image = color.rgb2grey(image)
    size_factor = float(height) / image.shape[0]
    new_size = [int(e * size_factor) for e in image.shape]
    image = transform.resize(image, new_size)
    image = filter.threshold_adaptive(image, block_size=30)

    return image

 def extract_histogram(self, bins=10):
     """Extract grey intensity histogram."""
     assert len(self.images) > 0, 'No images loaded! Did you call ' \
                                  'load_images() ?'
     histograms = []
     for image in self.images:
         grey = skicol.rgb2grey(image)
         hist_values, bins = np.histogram(grey, range=(0, 1), bins=bins)
         histograms.append(hist_values)
     histograms = np.array(histograms)
     histograms = histograms.astype('float')
     return histograms

def watershed_counter(path_image):
    # load the image and convert it to a floating point data type
    rgb_image = img_as_float(io.imread(path_image))
    image = rgb2grey(rgb_image)

    bin_image = image > threshold_otsu(image)
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5))
    ax1.imshow(image, cmap=plt.cm.gray)
    ax1.set_title('Original Image')
    ax1.axis('off')

    ax2.hist(image)
    ax2.set_title('Otsu Thresholded Histogram')
    ax2.axvline(threshold_otsu(image), color='r')

    ax3.imshow(bin_image, cmap=plt.cm.gray)
    ax3.set_title('Thresholded Image')
    ax3.axis('off')

    # Now we want to separate the two objects in image
    # Generate the markers as local maxima of the distance to the background

    distance = ndi.distance_transform_edt(bin_image)
    local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), labels=bin_image)
    markers = ndi.label(local_maxi)[0]
    labels = watershed(distance, markers, mask=bin_image)

    regions = regionprops(labels)
    regions = [r for r in regions if r.area > 50]
    num = len(regions)

    fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7))
    ax0, ax1, ax2, ax3 = axes

    ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
    ax0.set_title('Overlapping objects')
    ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
    ax1.set_title('Distances')
    ax2.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
    ax2.set_title(str(num) + ' Total Objects')
    ax3.imshow(rgb_image, cmap=plt.cm.gray, interpolation='nearest')
    ax3.contour(labels, [0.5], linewidths=1.2, colors='y')
    ax3.axis('off')

    for ax in axes:
        ax.axis('off')

    fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
                        right=1)

    print num

    plt.show()