def handle(self, *args, **options):
        better_thans = BetterThan.objects.all() #.filter(pk__lte=50)

        ds = SupervisedDataSet(204960, 1)
        for better_than in better_thans:
            bt = imread(better_than.better_than.image.file)
            wt = imread(better_than.worse_than.image.file)
            better_than.better_than.image.file.close()
            better_than.worse_than.image.file.close()

            bt = filters.sobel(bt)
            wt = filters.sobel(wt)

            bt_input_array = np.reshape(bt, (bt.shape[0] * bt.shape[1]))
            wt_input_array = np.reshape(wt, (wt.shape[0] * wt.shape[1]))
            input_1 = np.append(bt_input_array, wt_input_array)
            input_2 = np.append(wt_input_array, bt_input_array)
            ds.addSample(np.append(bt_input_array, wt_input_array), [-1])
            ds.addSample(np.append(wt_input_array, bt_input_array), [1])
        
        net = buildNetwork(204960, 2, 1)

        train_ds, test_ds = ds.splitWithProportion(options['train_test_split'])
        _, test_ds = ds.splitWithProportion(options['test_split'])

        trainer = BackpropTrainer(net, ds)

        avgerr = trainer.testOnData(dataset=test_ds)
        print 'untrained avgerr: {0}'.format(avgerr)

        trainer.train()

        avgerr = trainer.testOnData(dataset=test_ds)
        print 'trained avgerr: {0}'.format(avgerr)

def nrOfEdgePixels(rgbimage, intensityImage):
    redEdges = sobel(rgbimage[:,:,0])
    grayEdges = sobel(intensityImage)
    t = redEdges - grayEdges
    t[t < 0.05] = 0
    t[t >= 0.05] = 1
    return convolve2d(t, np.ones((17,17)), mode="same")

def getSubImages(img, pixels, size):
    subImages = []
    originals = []
    for i in range(len(img)):
        subImageRow = []
        originalRow = []
        for j in range(len(img[i])):
            if i % pixels == 0 and j % pixels == 0 and i+size-1 < len(img) and j+size-1 < len(img[i]):
                subImage = []
                for k in range(i, i+size, int(size/20)):
                    line = []
                    for l in range(j, j+size, int(size/20)):
                        line.append(img[k][l])
                    subImage.append(line)
                originalRow.append(subImage)
                if preprocess == preprocessing.SOBEL:
                    subImage = denoise_bilateral(subImage, sigma_range=0.1, sigma_spatial=15)
                    subImage = sobel(subImage)
                elif preprocess == preprocessing.HOG:
                    subImage = useHoG(subImage)
                else:
                    subImage = denoise_bilateral(subImage, sigma_range=0.1, sigma_spatial=15)
                    subImage = sobel(subImage)
                    subImage = useHoG(subImage)
                subImageRow.append(subImage)
        if len(subImageRow) > 0:
            subImages.append(subImageRow)
            originals.append(originalRow)
    return subImages, originals

def sobel(data, sliceId=2):
    edges = np.zeros(data.shape)
    if sliceId == 2:
        for idx in range(data.shape[2]):
            edges[:, :, idx] = skifil.sobel(data[:, :, idx])
    elif sliceId == 0:
        for idx in range(data.shape[0]):
            edges[idx, :, :] = skifil.sobel(data[idx, :, :])
    return edges

def color_edge():
        image=data.astronaut()
        r=np.abs(filters.sobel(image[:,:,0]))
        r=np.uint8(r/r.max()*255)
        io.imsave('astronautedger.png',r)
        g=np.abs(filters.sobel(image[:,:,1]))
        g=np.uint8(g/g.max()*255)
        io.imsave('astronautedgeg.png',g)
        b=np.abs(filters.sobel(image[:,:,2]))
        b=np.uint8(b/b.max()*255)
        io.imsave('astronautedgeb.png',b)

def filter_bank(img, coeff_resolution):
    """
    Calculates the responses of an image to M filters.
    Returns 2-d array of the vectorial responses
    """

    h, w = img.shape

    im = np.reshape(img, (h*w, 1))

    e1 = np.reshape(entropy(img, disk(coeff_resolution*5)), (h*w, 1))
    e2 = np.reshape(entropy(img, disk(coeff_resolution*8)), (h*w, 1))
    e3 = np.reshape(entropy(img, disk(coeff_resolution*10)), (h*w, 1))

    g1 = np.reshape(gradient(img, disk(1)), (h*w, 1))
    g2 = np.reshape(gradient(img, disk(coeff_resolution*3)), (h*w, 1))
    g3 = np.reshape(gradient(img, disk(coeff_resolution*5)), (h*w, 1))

    m1 = np.reshape(ndi.maximum_filter(256-img, size=coeff_resolution*2, mode='constant'), (h*w, 1))
    m2 = np.reshape(ndi.maximum_filter(256-img, size=coeff_resolution*4, mode='constant'), (h*w, 1))
    m3 = np.reshape(ndi.maximum_filter(256-img, size=coeff_resolution*7, mode='constant'), (h*w, 1))

    #c = np.reshape(canny(img), (h*w, 1))
    s = np.reshape(sobel(img), (h*w, 1))

    return np.column_stack((im, e1, e2, e3, g1, g2, g3, m1, m2, m3, s))

def pestFeatureExtraction(filename):
	selem = disk(8)
	image = data.imread(filename,as_grey=True)
	thresh = threshold_otsu(image)
	elevation_map = sobel(image)
	markers = np.zeros_like(image)

	if ((image<thresh).sum() > (image>thresh).sum()):
		markers[image < thresh] = 1
		markers[image > thresh] = 2
	else:
		markers[image < thresh] = 2
		markers[image > thresh] = 1

	segmentation = morphology.watershed(elevation_map, markers)
	segmentation = dilation(segmentation-1, selem)
	segmentation = ndimage.binary_fill_holes(segmentation)

	segmentation = np.logical_not(segmentation)
	image[segmentation]=0;

	hist = np.histogram(image.ravel(),256,[0,1])

	hist = list(hist[0])
	hist[:] = [float(x) / (sum(hist) - hist[0]) for x in hist]
	hist.pop(0)

	features = np.empty( (1, len(hist)), 'float' )
	
	a = np.array(list(hist))
	f = a.astype('float')
	features[0,:]=f[:]

	return features

def prepare(img):
    """
    Pre-process the image before translation detection, here we transform to black and white and use edge-detection.
    :param img: An image (as numpy array)
    :return: The preprocessed image (as numpy array)
    """
    return sobel(rgb2gray(img))

def testSkimage():
    img = Image.open('../img/1.png')
    img = np.array(img)
    imggray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    # (thresh, imgbw) = cv2.threshold(imggray, 128, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)

    # canny detector
    # from skimage.feature import canny
    # edges = canny(imggray/ 255.)
    from scipy import ndimage as ndi
    # fill_imgbw = ndi.binary_fill_holes(edges)
    # label_objects, nb_labels = ndi.label(fill_imgbw)
    # sizes = np.bincount(label_objects.ravel())
    # mask_sizes = sizes > 20
    # mask_sizes[0] = 0
    # cleaned_imgbw = mask_sizes[label_objects]

    markers = np.zeros_like(imggray)
    markers[imggray < 120] = 1
    markers[imggray > 150] = 2

    from skimage.filters import sobel
    elevation_map = sobel(imggray)
    from skimage.morphology import watershed
    segmentation = watershed(elevation_map, markers)

    # from skimage.color import label2rgb
    # segmentation = ndi.binary_fill_holes(segmentation - 10)
    # labeled_coins, _ = ndi.label(segmentation)
    # image_label_overlay = label2rgb(labeled_coins, image=imggray)
    plt.imshow(segmentation, cmap='gray')
    plt.show()
    return

def _segment_watershed(image):
	elevation_map = sobel(image)
	markers = np.zeros(image.shape) # initialize markers as zero array 

	
	# determine thresholds for markers
	sorted_pixels = np.sort(image, axis=None)
	max_int = np.mean(sorted_pixels[-10:])
	min_int = np.mean(sorted_pixels[:10])
	#max_int = np.max(orig_image)
	#min_int = np.min(orig_image)
	
	alpha_min = 0.01
	alpha_max = 0.4
	thresh_background = (1-alpha_min)*min_int	+	alpha_min*max_int
	thresh_spots = 		(1-alpha_max)*min_int	+	alpha_max*max_int
	
	markers[image < thresh_background] = 1 # mark background
	markers[image > thresh_spots] = 2 # mark background
	
	segmentation = watershed(elevation_map, markers)
	segmentation = segmentation-1
	segmentation = ndi.binary_fill_holes(segmentation)	# fill holes
	
	return segmentation

    def _calc_crispness(self, grey_array):
        """Calculate three measures of the crispness of an channel.

        PARAMETERS
        ----------
        grey_array : 2D numpy array
            Raw data for the grey channel.

        PRODUCES
        --------
        crispnesses : list
            Three measures of the crispness in the grey channel of types:

            - ``sobel``, ``canny``, and ``laplace``
        """
        grey_array = grey_array/255
        sobel_var = filters.sobel(grey_array).var()
        canny_array = feature.canny(grey_array, sigma=1).var()
        canny_ratio = np.sum(canny_array == True)/float(
                                                    len(canny_array.flatten()))
        laplace_var = filters.laplace(grey_array, ksize=3).var()
        self.feature_data.extend([sobel_var, canny_ratio, laplace_var])
        if self.columns_out:
            self.column_names.extend(['crisp_sobel', 'crisp_canny',
                                      'crisp_laplace'])

def filter(data,filtType,par):

    if   filtType == "sobel":       filt_data = sobel(data)
    elif filtType == "roberts":     filt_data = roberts(data)
    elif filtType == "canny":       filt_data = canny(data)
    elif filtType == "lowpass_avg":
        from scipy import ndimage
        p=int(par)
        kernel = np.ones((p,p),np.float32)/(p*p)
        filt_data = ndimage.convolve(data, kernel)
    elif filtType == "highpass_avg":
        from scipy import ndimage
        p=int(par)
        kernel = np.ones((p,p),np.float32)/(p*p)
        lp_data = ndimage.convolve(data, kernel)
        filt_data = data - lp_data
    elif filtType == "lowpass_gaussian":
        filt_data = gaussian(data, sigma=float(par))
    elif filtType == "highpass_gaussian":
        lp_data   = gaussian(data, sigma=float(par))
        filt_data = data - lp_data

    #elif filtType ==  "gradient":
       
    return filt_data

 def transform(self,X):
     imgs = []
     for x in X:
         if x.ndim == 3:
             x =self.rgb2gray(x)
         imgs.append(sobel(x).ravel())
     return np.vstack(imgs)

def op_vs_ip(subid, image_types, imagepaths, op_direc, overlays):
	
	
	img_data_group=[]
	img_shape_group=[]
	ol_data_group=[]
	ol_shape_group=[]
	for i, path in enumerate(imagepaths):	

		axial_slice, cor_slice, sag_slice, img_aspect_axial, img_aspect_cor, img_aspect_sag = pull_midslices(path)
		if os.path.isfile(overlays[i]):
			axial_slice_ol, cor_slice_ol, sag_slice_ol, img_aspect_axial_ol, img_aspect_cor_ol, img_aspect_sag_ol = pull_midslices(overlays[i])
			ol_data_group.append([axial_slice_ol, cor_slice_ol, sag_slice_ol])
			ol_shape_group.append([img_aspect_axial_ol, img_aspect_cor_ol, img_aspect_sag_ol])
		else:
			ol_data_group.append(['null','null','null'])
			ol_shape_group.append(['null','null','null'])
		## Append to Matrices
		img_data_group.append([axial_slice, cor_slice, sag_slice])
		img_shape_group.append([img_aspect_axial,img_aspect_cor,img_aspect_sag])
		


	my_cmap=plt.cm.gray


	fig, axarr = plt.subplots(ncols=np.shape(img_shape_group)[1], nrows=np.shape(img_shape_group)[0], figsize=(np.shape(img_shape_group)[0]*5,np.shape(img_shape_group)[1]*5))
	plt.suptitle(subid+' File Comparison', fontsize=20)	
	
	titlearray=['Axial', 'Coronal', 'Saggital']
	
	for x in range(0,np.shape(img_shape_group)[0]):
		for y in range(0,np.shape(img_shape_group)[1]):
			im = axarr[x, y].imshow(img_data_group[x][y], cmap=my_cmap, aspect=img_shape_group[x][y])
			axarr[x, y].set_xlabel('(Right) Radiological Convention (Left)', fontsize=10)
			axarr[x, y].set_title(image_types[x]+' '+titlearray[y])
			#divider = make_axes_locatable(axarr[x, y])
			#cax_ = divider.append_axes("right", size="5%", pad=0.05)
			#cbar = plt.colorbar(im, cax=cax_, ticks=MultipleLocator(round(np.max(img_data_group[x][y])/5, 1)))
			axarr[x, y].xaxis.set_visible(False)
			axarr[x, y].yaxis.set_visible(False)




			if os.path.isfile(overlays[x]):
				if x == 1:
					thresh=0.25
				if x == 2:
					thresh=0.4
				sl=np.array(ol_data_group[x][y]).astype(np.float64)
				sl=filters.sobel(sl)
				sl=preprocessing.binarize(sl, np.max(sl)*thresh)
				sl[sl < 1] = 'Nan'
				axarr[x, y].imshow(sl, cmap='autumn', aspect=ol_shape_group[x][y])

	#plt.show()
	plt.tight_layout()
	plt.autoscale()
	plt.savefig(op_direc)

def main():
    """Load image, apply sobel (to get x/y gradients), plot the results."""
    img = data.camera()

    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    sobel_x = np.rot90(sobel_y)  # rotates counter-clockwise

    # apply x/y sobel filter to get x/y gradients
    img_sx = signal.correlate(img, sobel_x, mode="same")
    img_sy = signal.correlate(img, sobel_y, mode="same")

    # combine x/y gradients to gradient magnitude
    # scikit-image's implementation divides by sqrt(2), not sure why
    img_s = np.sqrt(img_sx ** 2 + img_sy ** 2) / np.sqrt(2)

    # create binarized image
    threshold = np.average(img_s)
    img_s_bin = np.zeros(img_s.shape)
    img_s_bin[img_s > threshold] = 1

    # generate ground truth (scikit-image method)
    ground_truth = skifilters.sobel(data.camera())

    # plot
    util.plot_images_grayscale(
        [img, img_sx, img_sy, img_s, img_s_bin, ground_truth],
        [
            "Image",
            "Sobel (x)",
            "Sobel (y)",
            "Sobel (magnitude)",
            "Sobel (magnitude, binarized)",
            "Sobel (Ground Truth)",
        ],
    )

def number_nucleus(image):

    elevation_map = sobel(image)
    markers = np.zeros_like(image)
    markers[image < 250] = 1
    markers[image > 2000] = 2

    segmentation = watershed(elevation_map, markers)
    label_img = label(segmentation)
    prop = regionprops(label_img)

    width, height = plt.rcParams['figure.figsize']
    plt.rcParams['image.cmap'] = 'gray'

    image_label_overlay = label2rgb(label_img, image=image)

    fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(15, 8))
    ax1.imshow(image_label_overlay)
    ax2.imshow(image, cmap=plt.cm.gray, interpolation='nearest')

    # create list of region with are < 1000
    image_labeled = [region for region in prop if region.area > 5000]


    return len(image_labeled)

def main():
    from skimage import data, io, filters
    testfolder='/Users/davidgreenfield/Downloads/pics_boots/'
    testimage='B00A0GVP8A.jpg'
    image = io.imread(testfolder+testimage,flatten=True) # or any NumPy array!
    edges = filters.sobel(image)
    io.imshow(edges)
    io.show()

def sobel(frame):
    """
    return the sobel importance of an rgb image
    """
    frame = grayscale(frame)
    frame = filters.sobel(frame)
    # print(frame.max())
    return normalize(frame)

def getContourImage(imageFile):
    planeImage = data.imread(imageFile,True)
    imageArray = np.asarray(planeImage)
    averageColor = np.mean(imageArray)
    imageArray = getBlackAndWhiteImage(imageArray,averageColor*0.85)
    imageArray = filters.sobel(imageArray)
    imageArray = morphology.dilation(imageArray,morphology.disk(3))
    return imageArray

def test_sobel_vertical():
    """Sobel on a vertical edge should be a vertical line."""
    i, j = np.mgrid[-5:6, -5:6]
    image = (j >= 0).astype(float)
    result = filters.sobel(image) * np.sqrt(2)
    j[np.abs(i) == 5] = 10000
    assert (np.all(result[j == 0] == 1))
    assert (np.all(result[np.abs(j) > 1] == 0))