def plot_img_and_hist(img, axes, bins=256):
    """Plot an image along with its histogram and cumulative histogram.

    """
    img = img_as_float(img)
    ax_img, ax_hist = axes
    ax_cdf = ax_hist.twinx()

    # Display image
    ax_img.imshow(img, cmap=plt.cm.gray)
    ax_img.set_axis_off()
    ax_img.set_adjustable('box-forced')

    # Display histogram
    ax_hist.hist(img.ravel(), bins=bins, histtype='step', color='black')
    ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
    ax_hist.set_xlabel('Pixel intensity')
    ax_hist.set_xlim(0, 1)
    ax_hist.set_yticks([])

    # Display cumulative distribution
    img_cdf, bins = exposure.cumulative_distribution(img, bins)
    ax_cdf.plot(bins, img_cdf, 'r')
    ax_cdf.set_yticks([])

    return ax_img, ax_hist, ax_cdf

def match_hist(image, cdf, bin_centers, nbins=256):
    '''Modify pixels of input image so that its histogram matches target image histogram, specified by:
    cdf, bin_centers = cumulative_distribution(target_image)

    Parameters
    ----------
    image : array
        Image to be transformed.
    cdf : 1D array
        Values of cumulative distribution function of the target histogram.
    bin_centers ; 1D array
        Centers of bins of the target histogram.
    nbins : int, optional
        Number of bins for image histogram.

    Returns
    -------
    out : float array
        Image array after histogram matching.

    References
    ----------
    [1] Matlab implementation histoMatch(MTX, N, X) by Simoncelli, 7/96.
    '''
    image = img_as_float(image)
    old_cdf, old_bin = exposure.cumulative_distribution(image, nbins) # Unlike [1], we didn't add small positive number to the histogram
    new_bin = np.interp(old_cdf, cdf, bin_centers)
    out = np.interp(image.ravel(), old_bin, new_bin)
    return out.reshape(image.shape)

def test_equalize_ubyte():
    with expected_warnings(['precision loss']):
        img = skimage.img_as_ubyte(test_img)
    img_eq = exposure.equalize_hist(img)

    cdf, bin_edges = exposure.cumulative_distribution(img_eq)
    check_cdf_slope(cdf)

def plithist(im, nbins=256):
    f = plt.figure(figsize=(6, 6))
    plt.hist(im.ravel(), bins=nbins, histtype='step', color='black')
    img_cdf, bins = exposure.cumulative_distribution(im, nbins)
    plt.plot(bins, img_cdf, 'r')
    plt.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
    plt.xlabel('Pixel intensity')
    plt.xlim(0, 1)
    plt.yticks([])

def test_equalize_masked():
    img = skimage.img_as_float(test_img)
    mask = np.zeros(test_img.shape)
    mask[50:150, 50:250] = 1
    img_mask_eq = exposure.equalize_hist(img, mask=mask)
    img_eq = exposure.equalize_hist(img)

    cdf, bin_edges = exposure.cumulative_distribution(img_mask_eq)
    check_cdf_slope(cdf)

    assert not (img_eq == img_mask_eq).all()

def plot_hist(img, bins=256):
    """Plot histogram and cumulative histogram for image"""
    img_cdf, bins = exposure.cumulative_distribution(img, bins)
    plt.hist(img.ravel(), bins=bins)
    plt.ylabel('Number of pixels')
    plt.xlabel('Pixel intensiy')

    ax_cdf = plt.twinx()
    ax_cdf.plot(bins, img_cdf, 'r')
    xmin, xmax = dtype_range[img.dtype.type]
    plt.xlim(xmin, xmax)

    ax_cdf.set_ylabel('Fraction of total intensity')

def plot_img_and_hist(img, axes, bins=256):
    img = img_as_float(img)
    ax_img, ax_hist = axes
    ax_cdf = ax_hist.twinx()
    ax_img.imshow(img, cmap=plt.cm.gray)
    ax_img.set_axis_off()
    ax_hist.hist(img.ravel(), bins=bins, histtype='step', color='black')
    ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
    ax_hist.set_xlabel('Pixel intensity')
    ax_hist.set_xlim(0, 1)
    ax_hist.set_yticks([])
    img_cdf, bins = exposure.cumulative_distribution(img, bins)
    ax_cdf.plot(bins, img_cdf, 'r')
    ax_cdf.set_yticks([])
    return ax_img, ax_hist, ax_cdf

def RFshow(data,ax=None,cmap='RdYlBu',cmap_norm='equalize',hs=True,
              zf=10,azdeg=45,altdeg=45,dx=1,dy=1,fraction=1.5,blend_mode='alpha',
              alpha=0.7,contours=False,levels=32,colorbar=True,cb_contours=False,
              cb_ticks='linear',nSigma=1,fRange=None,FFTbins=512,**kwargs):

    # equalize the colormap:
    cdf, bins = exposure.cumulative_distribution(data[~np.isnan(data)].flatten(),nbins=256)
    my_cmap = equalizeColormap(cmap,bins,cdf)

    # convert input data to masked array
    data = np.ma.masked_array(data, np.isnan(data))

    fig = ax.get_figure()

    im = ax.imshow(data,cmap=my_cmap,**kwargs)

    # time axis:
    ax.xaxis_date()
    xxstart, xxend = ax.get_xlim()
    # we want 11 ticks:
    ax.xaxis.set_ticks(np.linspace(xxstart, xxend, 11, endpoint=True))
    #ax.set_xticklabels(np.linspace(xxstart, xxend, 11, endpoint=True))
    date_format = mdates.DateFormatter('%H:%M:%S')
    ax.xaxis.set_major_formatter(date_format)
    # fig.autofmt_xdate()
    for tick in ax.get_xticklabels():
    	tick.set_rotation(90)

    # frequencies axis:
    # we want 11 ticks:
    plt.yticks(np.linspace(fRange[0], fRange[1], 11, endpoint=True))
    yformatter = FuncFormatter(fmtMHz)
    ax.yaxis.set_major_formatter(yformatter)

    plt.xlabel('UTC time', labelpad=10)
    plt.ylabel('FFT frequencies in MHz (over %d bins)' % (FFTbins) )

    # colorbar:
    newTicks = stats_boundaries(data,nSigma,nSigma)
    cb1 = fig.colorbar(im,ticks=newTicks)

    cb1.ax.set_xlabel('relative\npower', labelpad=10)

    cb1.update_normal(im)

def histogramEqualize(image, maxIntensity):
  width, height = image.shape

  # Get cdf from image.
  cdf, binCenters = exposure.cumulative_distribution(image, maxIntensity)
  binCenters = binCenters.tolist()
  
  # Intensity transformation: Each pixel becomes the cumulative probability
  # that its intensity will show up, multiplied by the intended maximum 
  # intensity.
  for i in range(width):
    for j in range(height):
        
        try:
          probability = cdf[binCenters.index(image[i][j])]
        except:
          probability = 1
        
        image[i][j] = int(probability * maxIntensity)
  
  return image, binCenters, cdf

def plot_img_and_hist(img, axes, bins=256):

    ax_img, ax_hist = axes
    ax_cdf = ax_hist.twinx()

    # Display image
    ax_img.imshow(img, cmap='gray')
    ax_img.set_axis_off()

    # Display histogram
    ax_hist.hist(img.ravel() * 255, bins=bins)
    ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
    ax_hist.set_xlabel('Pixel intensity')


    ax_hist.set_xlim(0, 255)
    # Display cumulative distribution
    img_cdf, bins = exposure.cumulative_distribution(img*255, bins)
    ax_cdf.plot(bins, img_cdf, 'r')

    return ax_img, ax_hist, ax_cdf

def EMD(saliency_map1, saliency_map2, sub_sample=1/32.0):
    '''
    Earth Mover's Distance measures the distance between two probability distributions
    by how much transformation one distribution would need to undergo to match another
    (EMD=0 for identical distributions).

    Parameters
    ----------
    saliency_map1 : real-valued matrix
        If the two maps are different in shape, saliency_map1 will be resized to match saliency_map2.
    saliency_map2 : real-valued matrix

    Returns
    -------
    EMD : float, positive
    '''
    map2 = np.array(saliency_map2, copy=False)
    # Reduce image size for efficiency of calculation
    map2 = resize(map2, np.round(np.array(map2.shape)*sub_sample), order=3, mode='nearest')
    map1 = resize(saliency_map1, map2.shape, order=3, mode='nearest')
    # Histogram match the images so they have the same mass
    map1 = match_hist(map1, *exposure.cumulative_distribution(map2))
    # Normalize the two maps to sum up to 1,
    # so that the score is independent of the starting amount of mass / spread of fixations of the fixation map
    map1 = normalize(map1, method='sum')
    map2 = normalize(map2, method='sum')
    # Compute EMD with OpenCV
    # - http://docs.opencv.org/modules/imgproc/doc/histograms.html#emd
    # - http://stackoverflow.com/questions/5101004/python-code-for-earth-movers-distance
    # - http://stackoverflow.com/questions/12535715/set-type-for-fromarray-in-opencv-for-python
    r, c = map2.shape
    x, y = np.meshgrid(range(c), range(r))
    signature1 = cv.CreateMat(r*c, 3, cv.CV_32FC1)
    signature2 = cv.CreateMat(r*c, 3, cv.CV_32FC1)
    cv.Convert(cv.fromarray(np.c_[map1.ravel(), x.ravel(), y.ravel()]), signature1)
    cv.Convert(cv.fromarray(np.c_[map2.ravel(), x.ravel(), y.ravel()]), signature2)
    return cv.CalcEMD2(signature2, signature1, cv.CV_DIST_L2)

def plot_img_and_hist(img, axes, bins=256):
    """Plot an image along with its histogram and cumulative histogram.

    """
    ax_img, ax_hist = axes
    ax_cdf = ax_hist.twinx()

    # Display image
    ax_img.imshow(img, cmap=plt.cm.gray)
    ax_img.set_axis_off()

    # Display histogram
    ax_hist.hist(img.ravel(), bins=bins)
    ax_hist.ticklabel_format(axis="y", style="scientific", scilimits=(0, 0))
    ax_hist.set_xlabel("Pixel intensity")

    xmin, xmax = dtype_range[img.dtype.type]
    ax_hist.set_xlim(xmin, xmax)

    # Display cumulative distribution
    img_cdf, bins = exposure.cumulative_distribution(img, bins)
    ax_cdf.plot(bins, img_cdf, "r")

    return ax_img, ax_hist, ax_cdf

def imshow_hs(data,ax=None,cmap='geosoft',cmap_norm='equalize',hs=True,
              zf=10,azdeg=45,altdeg=45,dx=1,dy=1,fraction=1.5,blend_mode='alpha',
              alpha=0.7,contours=False,levels=32,colorbar=True,cb_contours=False,
              cb_ticks='linear',nSigma=1,**kwargs):
    '''
    Display an array with optional hillshading and contours. The colormap can be
    normalised by equalisation or by clipping extremes (autolevels).
    
    Parameters
    ----------
    data : 2D array
        Grid to plot. Arrays with NaNs and masked arrays are supported.
    ax : matplotlib axes instance
        This indicates where to draw the figure. Create new figure if absent.
    cmap : string
        Name of the colormap to use to display the array. The default 'geosoft' is
        the blue to pink clra colormap from Geosoft Oasis Montaj.
    cmap_norm : string
        Type of normalisation of the colormap. 
        Possible values are:
            'equalize' (or 'equalization')
                Increases contrast by distributing intensities across all the 
                possible colours. With this option, it is not the data that is normalised 
                but the colormap, based on the data. 
            'auto' (or 'autolevels')
                Stretches the histogram of the colormap so that dark colours become
                darker and the bright colours become brighter. Two extra parameters control 
                the amount of clipping at the extremes: minPercent (default to 10%) and
                maxPercent (default to 90%)
    hs : boolean
        If True, the array is displayed in colours over a grey hillshaded version
        of the data.
    zf : number
        Vertical exaggeration (Z factor) for hillshading.
    azdeg : number
        The azimuth (0-360, degrees clockwise from North) of the light source.
    altdeg : number
        The altitude (0-90, degrees up from horizontal) of the light source.
    dx : number, optional
        cell size in the x direction
    dy : number, optional
        cell size in the y direction
    fraction : number
        Increases or decreases the contrast of the hillshade. 
    blend_mode :  {'alpha', 'hsv', 'overlay', 'soft'} 
        The type of blending used to combine the colormapped data values with the 
        illumination intensity. Default is 'alpha' and the effect is controlled
        by the alpha parameter.
    alpha : float
        Controls the transparency of the data overlaid over the hillshade.
        1.0 is fully opaque while 0.0 is fully transparent.
    contours : Boolean
        If True, adds contours to the map. The number of calculated contours is 
        defined by:
            levels : integer
                Number of contour levels.
    colorbar : Boolean
        If True, draw a colorbar on the right-hand side of the map. The colorbar
        shows the distribution of colors, as modified by the normalization algorithm.
    cb_ticks : string
        If left as default ('linear') the ticks and labels on the colorbar are 
        spaced linearly in the standard way. Otherwise (any other keyword, for example
        'stats'), the mean and two ticks at + and - nSigma*(standard deviation) 
        are shown instead.
            nSigma : integer (default is 1)
                Size of the interval to show between ticks on the colorbar. 
    cb_contours : Boolean
        Add lines corresponding to contours on the colorbar.
    kwargs : other optional arguments
        Can be used to pass other arguments to imshow, such as 'origin' and 'extent'.
        
    Notes
    -----
    This function exploits the hillshading capabilities implemented in
    matplotlib.colors.LightSource. It adds additional blending mode (alpha compositing,
    see https://en.wikipedia.org/wiki/Alpha_compositing) and normalising functions
    for the data (equalization).

    '''
    # modify colormap if required
    if cmap_norm in ['equalize','equalization']:
        # histogram equalization
        cdf, bins = exposure.cumulative_distribution(data[~np.isnan(data)].flatten(),nbins=256)
        my_cmap = equalizeColormap(cmap,bins,cdf)
    elif cmap_norm in ['auto','autolevels']:
        # autolevels
        minP = kwargs.pop('minPercent',10) # also removes the key from the dictionary
        maxP = kwargs.pop('maxPercent',90)
        my_cmap = normalizeColormap(cmap,norm='autolevels',minPercent=minP,maxPercent=maxP)
    elif cmap in plt.colormaps():
        # colormap defined as string (recognised name)
        my_cmap = plt.get_cmap(cmap) 
    else:
        # colormap is one of the extra ones added by the colors module 
        my_cmap = load_cmap(cmap)  # raises error if not recognised
        
    # create figure or retrieve the one already there
    if ax:
        fig = ax.get_figure()
    else:
#.........这里部分代码省略.........

def grayhist(img, *args, **histkwargs):
    """Plot an image along with its histogram and cumulative histogram.

    ADAPTED FROM SCIKIT IMAGE GALLERY
    http://scikit-image.org/docs/dev/auto_examples/plot_local_equalize.html

    Parameters
    ----------
    bins : (Number bins, defaults to 256)

    cdf : bool(False) or str(color)
        Plot cumulative distribution function over histogram.
    If cdf = color, interpreted as line color eg (cdf = 'r') 
    plots a red line for CDF.   

    lw / ls : CDF Line styles

    xlim : set (xs, xf) or "auto" 
        Return cropped histogram between x-limits.  If "auto", min and max
    brigntess of image are used.  

    Returns
    -------
    tuple : (n, bins, patches) or ([n0, n1, ...], bins, [patches0, patches1,...])

    Notes
    -----
    Unlike standard histogram, this returns axes rather than the
    histogram parameters.  Because this method changes api for xlim,
    IE user can prescribe xlimits through call signature, it is easier to just
    crop the image instead of changing the plot limits to account for the
    various cases.  Therefore, it would return output for cropped image
    histogram, which could lead to confusion.

    See matplotlib hist API for all plt.hist() parameters.
    http://matplotlib.org/api/pyplot_api.html
    """

    if img.ndim == 3:
        img = rgb2uint(img, warnmsg = True)

    # Histogram plotting kwargs
    bins = histkwargs.pop('bins', 256) #used several places
    cdf = histkwargs.pop('cdf', False)
    title = histkwargs.pop('title', None)
    histkwargs.setdefault('color', 'black')
    histkwargs.setdefault('alpha', 0.5)
    histkwargs.setdefault('orientation', 'vertical')

    # CDF line plotting kwargs
    lw = histkwargs.pop('lw', 2)    
    ls = histkwargs.pop('ls', '-')

    xlim = histkwargs.pop('xlim', None)

    # Set the range based on scikit image dtype range 
    # (not quite right for rgb)
    xmin, xmax = pp_dtype_range(img)

    if xlim:
        # ALSO SET VLIM FROM AUTO!
        if xlim =='auto':
            xlim = img.min(), img.max()

        rmin, rmax = xlim
        if rmin < xmin or rmax > xmax:
            raise UtilsError("Range %s out of bounds (%s, %s)" %
                             (xlim, xmin, xmax))
        else:
            xmin, xmax = xlim    

    raveled_img = img[(img >= xmin) & (img <= xmax)]

    if histkwargs['orientation'] == 'horizontal':
        raise UtilsError("horizontal orientation not supported.")

    axes, kwargs = _parse_ax(*args, **histkwargs)    

    # Matplotlib
    if not axes:
        fig, axes = plt.subplots()

    # Display histogram
    histout = axes.hist(raveled_img, bins=bins, **histkwargs)
    axes.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
    axes.set_xlabel('Pixel intensity')

    # Display cumulative distribution
    if cdf:
        if cdf is not True:
            lcolor = cdf
        else:
            lcolor = 'r'
        ax_cdf = axes.twinx()
        img_cdf, bins = exposure.cumulative_distribution(img, bins)
        ax_cdf.plot(bins, img_cdf, color=lcolor, lw=lw, ls=ls)

    axes.set_xlim(xmin, xmax) #is necessary
    if title:
        axes.set_title(title)
#.........这里部分代码省略.........

ax3.imshow(matched)
ax3.set_title('Matched')

plt.tight_layout()
plt.show()


######################################################################
# To illustrate the effect of the histogram matching, we plot for each
# RGB channel, the histogram and the cumulative histogram. Clearly,
# the matched image has the same cumulative histogram as the reference
# image for each channel.

fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(8, 8))


for i, img in enumerate((image, reference, matched)):
    for c, c_color in enumerate(('red', 'green', 'blue')):
        img_hist, bins = exposure.histogram(img[..., c], source_range='dtype')
        axes[c, i].plot(bins, img_hist / img_hist.max())
        img_cdf, bins = exposure.cumulative_distribution(img[..., c])
        axes[c, i].plot(bins, img_cdf)
        axes[c, 0].set_ylabel(c_color)

axes[0, 0].set_title('Source')
axes[0, 1].set_title('Reference')
axes[0, 2].set_title('Matched')

plt.tight_layout()
plt.show()

def test_equalize_float():
    img = skimage.img_as_float(test_img)
    img_eq = exposure.equalize_hist(img)

    cdf, bin_edges = exposure.cumulative_distribution(img_eq)
    check_cdf_slope(cdf)

from skimage.data import camera, moon
from skimage.exposure import cumulative_distribution, equalize_hist
from skimage.io import imshow, show
from numpy import vectorize

im_camera = camera()
im_moon = moon()
dist_camera, bins = cumulative_distribution(im_camera)
dist_moon, bins = cumulative_distribution(im_moon)


# def rechercher_transformation(dist_source, dist_target):
#     nb_source = len(dist_source)
#     nb_target = len(dist_target)
#     transformation = [nb_source - 1] * nb_source
#     i, j = 0, 0
#     while j < nb_target:
#         while i < nb_source and dist_source[i] < dist_target[j]:
#             transformation[i] = j
#             i += 1
#         j += 1
#     return transformation


def rechercher_transformation(dist_source, dist_target):
    nb_source = len(dist_source)
    nb_target = len(dist_target)
    transformation = [nb_source - 1] * nb_source
    i, j = 0, 0
    ind0, val0 = 0, 0
    while j < nb_target:

def plot_cdf(image, ax=None):
    img_cdf, bins = exposure.cumulative_distribution(image)
    ax.plot(bins, img_cdf, 'r')
    ax.set_ylabel("Fraction of pixels below intensity")

def egalisation(image):
    dist_image, bins = cumulative_distribution(image)
    nb = len(dist_image)
    transfo = rechercher_transformation(dist_image, [k / nb for k in range(nb)])
    return appliquer_transformation(image, transfo)

def match(im_source, im_target):
    dist_source, bins = cumulative_distribution(im_source)
    dist_target, bins = cumulative_distribution(im_target)
    transfo = rechercher_transformation(dist_source, dist_target)
    return appliquer_transformation(im_source, transfo)