def color_feature(blur, hbins=15, sbins=15):
    hsv = color.rgb2hsv(blur)

    # cal hist
    h_hist = exposure.histogram(hsv[:, :, 0], nbins=hbins)
    s_hist = exposure.histogram(hsv[:, :, 1], nbins=sbins)

    return np.append(normalize(h_hist[0]), normalize(s_hist[0]))

def getColorVector(im, nbin):
    h1, v1 = exposure.histogram(im[:,:,0], nbin)
    h2, v2 = exposure.histogram(im[:,:,1], nbin)
    h3, v3 = exposure.histogram(im[:,:,2], nbin)
    h1 = h1 / (h1.sum() * 1.0)
    h2 = h2 / (h2.sum() * 1.0)
    h3 = h3 / (h3.sum() * 1.0)
    return np.append(h1,[h2,h3])

def plotImageWithHistogram(im, size, alpha=0.3, interpolation='nearest'):
    r"""Plot an image alongside its histogram

    Parameters
    ----------

    im : a numpy array
        The input image

    size : number
        The size (in in) for the single figure. The plot will be
        that high and twice that wide.

    alpha : float
        The transparency. A value of 0.3 is great for an RGB image,
        a value of 0.8 is a bit better for a grayscale image.

    Returns
    -------
    (ax__image, ax_hist) 
        The matplotlib axes for the figure.

    Examples
    --------

    from jmToolsPy3 import plotImageWithHistogram
    import numpy as np
    from skimage import data

    img1 = data.camera()
    axImg1, axHis1 = plotImageWithHistogram(img1, 5, alpha=0.8)

    img2 = data.lena()
    axImg2, axHis2 = plotImageWithHistogram(img2, 5, alpha=0.3)
    """
    from skimage import exposure
    from matplotlib import pyplot as plt
    fig, (ax_image, ax_hist) = plt.subplots(ncols=2, figsize=(2*size, size))
    ax_image.imshow(im, cmap=plt.cm.gray, interpolation=interplolation)
    if im.ndim == 2:
        hist, bin_centers = exposure.histogram(im)
        ax_hist.fill_between(bin_centers, hist, alpha=alpha, color='gray')        
    elif im.ndim == 3:
        for channel, channel_color in zip(iter_channels(im), 'rgb'):
            hist, bin_centers = exposure.histogram(channel)
            ax_hist.fill_between(bin_centers, hist, alpha=alpha, color=channel_color)
    ax_hist.set_ylabel('# pixels')
    ax_hist.set_xlabel('intensity')
    ax_hist.set_yticklabels("")
    ax_image.set_axis_off()
    # match_axes_height(ax_image, ax_hist)
    dst = ax_hist.get_position()
    src = ax_image.get_position()
    ax_hist.set_position([dst.xmin, src.ymin, dst.width, src.height])
    return ax_image, ax_hist

def test_normalize():
    im = np.array([0, 255, 255], dtype=np.uint8)
    frequencies, bin_centers = exposure.histogram(im, source_range='dtype',
                                                  normalize=False)
    expected = np.zeros(256)
    expected[0] = 1
    expected[-1] = 2
    assert_equal(frequencies, expected)
    frequencies, bin_centers = exposure.histogram(im, source_range='dtype',
                                                  normalize=True)
    expected /= 3.
    assert_equal(frequencies, expected)

 def threshHist(imgIn):
     
     imgIn1 = imgIn.ravel()
     
     # Histogram analysis to find maximum peak and associated gl
     pxCnt, gL = exposure.histogram( imgIn )
     indMaxBG = int(np.arange( len(imgIn1) )[pxCnt == pxCnt.max()]) #int()
     BGlevel = gL[indMaxBG]
     
     # Nearest min below this max is threshold 
     d1 = np.zeros( np.shape(pxCnt) )
       
     for i in range( 2 , len(pxCnt) - 1):
         # derivative approximation
         d1[i] = pxCnt[ i + 1 ] - pxCnt[ i ]
     
     i = 1
     p = 0
     while ( d1[ indMaxBG - i ] > 0): ### - i!!!
         p = indMaxBG - i
         i = i + 1
     
     t = gL[ p ]
    
     imgOut = imgProc.applyThresh( imgIn, t ) 
     return imgOut

def threshold_yen(image, nbins=256, shift=None):
    """Return threshold value based on Yen's method.

    Parameters
    ----------
    image : array
        Input image.
    nbins : int, optional
        Number of bins used to calculate histogram. This value is ignored for
        integer arrays.
    shift : int, optional
        Shift threshold value by percent up (positive) or down (negative).

    Returns
    -------
    threshold : float
        Upper threshold value. All pixels intensities that less or equal of
        this value assumed as foreground.

    References
    ----------
    .. [1] Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion
           for Automatic Multilevel Thresholding" IEEE Trans. on Image
           Processing, 4(3): 370-378
    .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
           Techniques and Quantitative Performance Evaluation" Journal of
           Electronic Imaging, 13(1): 146-165,
           http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf
    .. [3] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold

    Examples
    --------
    >>> from skimage.data import camera
    >>> image = camera()
    >>> thresh = threshold_yen(image)
    >>> binary = image <= thresh
    """
    hist, bin_centers = histogram(image, nbins)
    # On blank images (e.g. filled with 0) with int dtype, `histogram()`
    # returns `bin_centers` containing only one value. Speed up with it.
    if bin_centers.size == 1:
        return bin_centers[0]

    # Calculate probability mass function
    pmf = hist.astype(np.float32) / hist.sum()
    P1 = np.cumsum(pmf)  # Cumulative normalized histogram
    P1_sq = np.cumsum(pmf ** 2)
    # Get cumsum calculated from end of squared array:
    P2_sq = np.cumsum(pmf[::-1] ** 2)[::-1]
    # P2_sq indexes is shifted +1. I assume, with P1[:-1] it's help avoid '-inf'
    # in crit. ImageJ Yen implementation replaces those values by zero.
    crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) *
                  (P1[:-1] * (1.0 - P1[:-1])) ** 2)

    threshold = bin_centers[crit.argmax()]
    ptp = (bin_centers[-1] - bin_centers[0]) / 100.  # Peek to peek range
    if shift:
        threshold += ptp * shift
    # print("Threshold value shift", (threshold - bin_centers[0]) / float(ptp))
    return threshold

def statxture(pixels):
    """computes a variety of texture stats from
    the image histogram.
    See Digital Image Processing Using MATLAB, ch. 11"""
    average_gray_level = np.mean(pixels)
    average_contrast = np.std(pixels)

    H = histogram(pixels)[0]
    H = H / (1. * len(pixels))
    L = len(H)

    d = (L - 1.)**2

    normvar = np.var(pixels) / d
    smoothness = 1. - 1. / (1. + normvar)

    third_moment = moment(pixels,3) / d

    uniformity = np.sum(H**2)

    eps = np.finfo(float).eps
    entropy = 0. - np.sum(H * np.log2(H + eps)) 
    
    return average_gray_level, average_contrast, smoothness, \
        third_moment, uniformity, entropy

def test_all_negative_image():
    im = np.array([-128, -1], dtype=np.int8)
    frequencies, bin_centers = exposure.histogram(im)
    assert_array_equal(bin_centers, np.arange(-128, 0))
    assert frequencies[0] == 1
    assert frequencies[-1] == 1
    assert_array_equal(frequencies[1:-1], 0)

def analyse_histogram(data, roi=None, debug=False, dens_min=20, dens_max=255, minT=0.95, maxT=1.05):
    if roi == None:
        #roi = np.ones(data.shape, dtype=np.bool)
        roi = np.logical_and(data >= dens_min, data <= dens_max)

    voxels = data[np.nonzero(roi)]
    hist, bins = skiexp.histogram(voxels)
    max_peakIdx = hist.argmax()

    minT = minT * hist[max_peakIdx]
    maxT = maxT * hist[max_peakIdx]
    histTIdxs = (hist >= minT) * (hist <= maxT)
    histTIdxs = np.nonzero(histTIdxs)[0]
    histTIdxs = histTIdxs.astype(np.int)

    class1TMin = bins[histTIdxs[0]]
    class1TMax = bins[histTIdxs[-1]]

    liver = data * (roi > 0)
    class1 = np.where( (liver >= class1TMin) * (liver <= class1TMax), 1, 0)

    if debug:
        plt.figure()
        plt.plot(bins, hist)
        plt.hold(True)

        plt.plot(bins[max_peakIdx], hist[max_peakIdx], 'ro')
        plt.plot(bins[histTIdxs], hist[histTIdxs], 'r')
        plt.plot(bins[histTIdxs[0]], hist[histTIdxs[0]], 'rx')
        plt.plot(bins[histTIdxs[-1]], hist[histTIdxs[-1]], 'rx')
        plt.title('Histogram of liver density and its class1 = maximal peak (red dot) +-5% of its density (red line).')
        plt.show()

    return class1

def test_peak_float_out_of_range_dtype():
    im = np.array([10, 100], dtype=np.float16)
    nbins = 10
    frequencies, bin_centers = exposure.histogram(im, nbins=nbins, source_range='dtype')
    assert_almost_equal(np.min(bin_centers), -0.9, 3)
    assert_almost_equal(np.max(bin_centers), 0.9, 3)
    assert_equal(len(bin_centers), 10)

def _plot_histogram(ax, image, alpha=0.3, **kwargs):
    # Use skimage's histogram function which has nice defaults for
    # integer and float images.
    hist, bin_centers = exposure.histogram(image)
    ax.fill_between(bin_centers, hist, alpha=alpha, **kwargs)
    ax.set_xlabel('intensity')
    ax.set_ylabel('# pixels')

def histogram(image):

    # Iterate throught a.) Colors, b.) Channels to create lines
    for color, channel in zip('rgb', np.rollaxis(image, axis=-1)):
        counts, bin_centers = exposure.histogram(channel)
        a = plt.fill_between(bin_centers, counts, color=color, alpha=0.4)
    return a;

    def set_data(self, data):#, mask):
        self.data = data
        # self.mask = mask
        # if self.data is not None:
        self.hist, self.bins = skiexp.histogram(self.data, nbins=1000)
            # if mask is not None:
            #     self.data_m = self.data[np.nonzero(self.mask)]
            # else:
            #     self.data_m = self.data

        if self.params and self.params.has_key('data_min'):
            self.data_min = self.params['data_min']
        # elif self.data is not None:
        self.data_min = self.data.min()
        # else:
        #     self.data_min = 0
        if self.params and self.params.has_key('datam_max'):
            self.data_max = self.params['data_max']
        # elif self.data is not None:
        self.data_max = self.data.max()
        # else:
        #     self.data_max = 0

        self.ui.hypo_mean_SL.setMinimum(self.data_min)
        self.ui.heal_mean_SL.setMinimum(self.data_min)
        self.ui.hyper_mean_SL.setMinimum(self.data_min)

        self.ui.hypo_mean_SL.setMaximum(self.data_max)
        self.ui.heal_mean_SL.setMaximum(self.data_max)
        self.ui.hyper_mean_SL.setMaximum(self.data_max)

        self.update_figures()

def test_negative_overflow():
    im = np.array([-1, 127], dtype=np.int8)
    frequencies, bin_centers = exposure.histogram(im)
    assert_array_equal(bin_centers, np.arange(-1, 128))
    assert frequencies[0] == 1
    assert frequencies[-1] == 1
    assert_array_equal(frequencies[1:-1], 0)

def testSkimage2():
    img = Image.open('../img/1.png')
    # img = Image.open('../img/2.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)
    camera = imggray
    # camera = data.camera()
    val = filters.threshold_otsu(camera)

    hist, bins_center = exposure.histogram(camera)

    plt.figure(figsize=(9, 4))
    plt.subplot(131)
    plt.imshow(camera, cmap='gray', interpolation='nearest')
    plt.axis('off')
    plt.subplot(132)
    plt.imshow(camera < val, cmap='gray', interpolation='nearest')
    plt.axis('off')
    plt.subplot(133)
    plt.plot(bins_center, hist, lw=2)
    plt.axvline(val, color='k', ls='--')

    plt.tight_layout()
    plt.show()
    return

def preprocImageV(fimg):
    nbin=64
    dn = 24
    img=skimage.img_as_float(io.imread(fimg,as_grey=True))
    siz=img.shape
    sizMax=np.max(siz)
    siz2=(sizMax,sizMax)
    msk=np.zeros(siz)
    cv2.circle(msk, (siz[1]/2, siz[0]/2), (np.min(siz)/2)-3, (1,1,1), -1) ##cv2.cv.CV_FILLED
    s0 = 0.3440
    h1,xx = exposure.histogram(img[msk==1], nbin)
    s1 = np.std(img[msk==1])
    h  = h1.copy()
    h[np.argwhere(h > 0)[-1]] = 0
    h[np.argwhere(h > 0)[-1]] = 0
    st = np.percentile(img[msk==1], 1)
    p1 = np.argwhere(np.fabs(xx - st)==np.min(abs(xx - st)))[0][0]
    h[p1:np.min((nbin, p1 + np.round(dn * s1 / s0)))]=0
    p2 = np.argwhere(h==np.max(h))[0][0]
    max1=xx[p1]
    max2=xx[p2]
    ret= 0.8*(img-max1)/(max2-max1)
    ret[ret<0]=0
    ret[ret>1]=1
    ret=np.uint8(255*ret)
    ret2=np.zeros(siz2, np.uint8)
    r0=np.floor((sizMax-siz[0])/2)
    c0=np.floor((sizMax-siz[1])/2)
    ret2[r0:r0+siz[0], c0:c0+siz[1]]=ret
    # cv2.imshow("ret", ret)
    # cv2.imshow("ret2", ret2)
    # cv2.waitKey(0)
    return ret2

def test_peak_int_range_dtype():
    im = np.array([10, 100], dtype=np.int8)
    frequencies, bin_centers = exposure.histogram(im, source_range='dtype')
    assert_array_equal(bin_centers, np.arange(-128, 128))
    assert_equal(frequencies[128+10], 1)
    assert_equal(frequencies[128+100], 1)
    assert_equal(frequencies[128+101], 0)
    assert_equal(frequencies.shape, (256,))

    def estimate_healthy_pdf(self, data, mask, params):
        perc = params['perc']
        k_std_l = params['k_std_h']
        simple_estim = params['healthy_simple_estim']
        show_me = params['show_healthy_pdf_estim']

        # data = tools.smoothing_tv(data.astype(np.uint8), tv_weight)
        ints = data[np.nonzero(mask)]
        hist, bins = skiexp.histogram(ints, nbins=256)
        if simple_estim:
            mu, sigma = scista.norm.fit(ints)
        else:
            ints = data[np.nonzero(mask)]

            n_pts = mask.sum()
            perc_in = n_pts * perc / 100

            peak_idx = np.argmax(hist)
            n_in = hist[peak_idx]
            win_width = 0

            while n_in < perc_in:
                win_width += 1
                n_in = hist[peak_idx - win_width:peak_idx + win_width].sum()

            idx_start = bins[peak_idx - win_width]
            idx_end = bins[peak_idx + win_width]
            inners_m = np.logical_and(ints > idx_start, ints < idx_end)
            inners = ints[np.nonzero(inners_m)]

            # liver pdf -------------
            mu = bins[peak_idx] + params['hack_healthy_mu']
            sigma = k_std_l * np.std(inners) + params['hack_healthy_sigma']

        mu = int(mu)
        sigma = int(sigma)
        rv = scista.norm(mu, sigma)

        if show_me:
            plt.figure()
            plt.subplot(211)
            plt.plot(bins, hist)
            plt.title('histogram with max peak')
            plt.hold(True)
            plt.plot([mu, mu], [0, hist.max()], 'g')
            ax = plt.axis()
            plt.axis([0, 256, ax[2], ax[3]])
            # plt.subplot(212), plt.plot(bins, rv_l.pdf(bins), 'g')
            x = np.arange(0, 256, 0.1)
            plt.subplot(212), plt.plot(x, rv.pdf(x), 'g')
            plt.hold(True)
            plt.plot(mu, rv.pdf(mu), 'go')
            ax = plt.axis()
            plt.axis([0, 256, ax[2], ax[3]])
            plt.title('estimated normal pdf of healthy parenchym')
            # plt.show()

        return rv

def features_extractor(trimmed):

    from color_quantinization import quantinization
    image_gray = quantinization(trimmed)
    # image_gray = color.rgb2hsv(image_gray)

    featuresh = []

    hh = exposure.histogram(image_gray[:, :, 0])
    featuresh.extend(hh[0])

    hs = exposure.histogram(image_gray[:, :, 1])
    featuresh.extend(hs[0])

    hv = exposure.histogram(image_gray[:, :, 2])
    featuresh.extend(hv[0])

    return featuresh

def test_dask_histogram():
    pytest.importorskip('dask', reason="dask python library is not installed")
    import dask.array as da
    dask_array = da.from_array(np.array([[0, 1], [1, 2]]), chunks=(1, 2))
    output_hist, output_bins = exposure.histogram(dask_array)
    expected_bins = [0, 1, 2]
    expected_hist = [1, 2, 1]
    assert np.allclose(expected_bins, output_bins)
    assert np.allclose(expected_hist, output_hist)