def anaPsf(xaxis,yaxis,data):

    # centroids
    xCentroid = (xaxis*data).sum()/data.sum()
    yCentroid = (yaxis*data).sum()/data.sum()

    # radial sum -- around Centroid
    raxis = numpy.sqrt((xaxis-xCentroid)*(xaxis-xCentroid)+(yaxis-yCentroid)*(yaxis-yCentroid))

    # histogram 
    nsumbin = 1000
    npix,bin_edges = scipy.histogram(raxis,nsumbin,(0.,100.))
    rsumpix,bin_edges = scipy.histogram(raxis,nsumbin,(0.,100.),weights=data)

    # calculate ee80
    rsumpixNorm = rsumpix/rsumpix.sum()
    rcumsum = numpy.cumsum(rsumpixNorm)
    # at this point rcumsum[0] is = rsumpixNorm[0], rsumsum[1] to the sum of rsumpixnorm[0] and [1] etc.
    # so rcumsum[0] is the integral for r<bin_edges[1]   and in general rcumsum[i] is the integral of r<bin_edges[i+1]
    # thus icumsum gives the appropriate limits for each rcumsum bin
    #

    icumsum = bin_edges[1:nsumbin+1]
    ee80 = numpy.interp(0.8,rcumsum,icumsum)

    # calculate polarization (w/o seeing, CCD diffusion)
    norm = data.sum()
    qxx = ( (xaxis-xCentroid)*(xaxis-xCentroid)*data ).sum() / norm
    qyy = ( (yaxis-yCentroid)*(yaxis-yCentroid)*data ).sum() / norm
    qyx = ( (yaxis-yCentroid)*(xaxis-xCentroid)*data ).sum() / norm
    
    e1 = (qxx-qyy)/(qxx+qyy)
    e2 = 2.0*qyx/(qxx+qyy)

    return ee80,qxx,qyy,qyx,e1,e2

def generate_scipy_comparison(csvPathname):
    # this is some hack code for reading the csv and doing some percentile stuff in scipy
    from numpy import loadtxt, genfromtxt, savetxt

    dataset = loadtxt(
        open(csvPathname, 'r'),
        delimiter=',',
        dtype='float64');

    print "csv read for training, done"

    # we're going to strip just the last column for percentile work
    # used below
    NUMCLASSES = 10
    print "csv read for training, done"

    # data is last column
    # drop the output
    print dataset.shape
    from scipy import histogram
    import numpy
    if 1==1:
        print "histogram of dataset"
        print histogram(dataset,bins=NUMCLASSES)
        print numpy.mean(dataset, axis=0, dtype=numpy.float64)
        print numpy.std(dataset, axis=0, dtype=numpy.float64, ddof=0)
        print numpy.std(dataset, axis=0, dtype=numpy.float64, ddof=1)

    from scipy import stats

def genMCHistogramsOpenCL(distribution, rng, iterations=100, numBins=1000):
    # get OpenCL CPU devices
    openCLDevices = [device for device in clsim.I3CLSimOpenCLDevice.GetAllDevices() if device.cpu]
    if len(openCLDevices)==0:
        raise RuntimeError("No CPU OpenCL devices available!")
    openCLDevice = openCLDevices[0]

    openCLDevice.useNativeMath=False
    workgroupSize = 1
    workItemsPerIteration = 10240
    print("           using platform:", openCLDevice.platform)
    print("             using device:", openCLDevice.device)
    print("            workgroupSize:", workgroupSize)
    print("    workItemsPerIteration:", workItemsPerIteration)
    
    tester = clsim.I3CLSimRandomDistributionTester(device=openCLDevice,
                                                   workgroupSize=workgroupSize,
                                                   workItemsPerIteration=workItemsPerIteration,
                                                   randomService=rng,
                                                   randomDistribution=distribution)
    
    print("maxWorkgroupSizeForKernel:", tester.maxWorkgroupSize)
    
    angles = tester.GenerateRandomNumbers(iterations)
    samples = len(angles)
    
    print("generated")
    
    angles = numpy.array(angles) # convert to numpy array
    print("converted")
    
    numAng_orig, binsAng = scipy.histogram(numpy.arccos(angles)*(180./math.pi), range=(0.,180.), bins=numBins)
    print("hist1 complete")

    numCos_orig, binsCos = scipy.histogram(angles, range=(-1.,1.), bins=numBins)
    print("hist2 complete")
    
    del angles # not needed anymore
    print("deleted")
    
    numAng=[]
    for i, number in enumerate(numAng_orig):
        binWidth = math.cos(binsAng[i]*math.pi/180.) - math.cos(binsAng[i+1]*math.pi/180.)
        numAng.append(float(number)/float(samples)/binWidth)
    numAng=numpy.array(numAng)
    
    numCos=[]
    for i, number in enumerate(numCos_orig):
        numCos.append(float(number)/float(samples)/float(2./float(numBins)))
    numCos=numpy.array(numCos)
    
    binsAng = numpy.array(binsAng[:-1])+(binsAng[1]-binsAng[0])/2.
    binsCos = numpy.array(binsCos[:-1])+(binsCos[1]-binsCos[0])/2.
    
    return dict(cos=dict(num=numCos, bins=binsCos), ang=dict(num=numAng, bins=binsAng))

def getDominantColor(img_url):
    if r.exists(img_url):
        cache_result = r.hmget(img_url, ['r', 'g', 'b'])
        return cache_result
        
    NUM_CLUSTERS = 5
    im = Image.open(StringIO.StringIO(urllib2.urlopen(img_url).read()))
    img_arr = scipy.misc.fromimage(im)
    img_shape = img_arr.shape
    
    if len(img_shape) > 2:
        img_arr = img_arr.reshape(scipy.product(img_shape[:2]), img_shape[2])
    
    codes, _ = scipy.cluster.vq.kmeans(img_arr, NUM_CLUSTERS)
    
    original_codes = codes
    for low, hi in [(60, 200), (35, 230), (10, 250)]:
        codes = scipy.array([code for code in codes if not (all([c < low for c in code]) or all([c > hi for c in code]))])
        if not len(codes):
            codes = original_codes
        else:
            break

    vecs, _ = scipy.cluster.vq.vq(img_arr, codes)
    counts, bins = scipy.histogram(vecs, len(codes))

    index_max = scipy.argmax(counts)
    peak = codes[index_max]
    color = [int(c) for c in peak[:3]]
    r.hmset(img_url, {'r':color[0], 'g':color[1], 'b':color[2]})
    #r.expire(img_url, 86400)
    return color

def hist(fcms, index, savefile=None, display=True, **kwargs):
    """Plot overlay histogram.

    fcms is a list of FCMData objects/arrays
    index is channel to plot
    """
    figure = pylab.figure()
    for fcm in fcms:
        if isinstance(index, str):
            index = fcm.name_to_index(index)
        y = fcm[:, index]
        h, b = histogram(y, bins=200, **kwargs)
        b = (b[:-1] + b[1:]) / 2.0
        unused_x = pylab.linspace(min(y), max(y), 100)
        pylab.plot(b, h, label=fcm.name)
    pylab.legend()
    pylab.xlabel(fcms[0].channels[index])

    if display:
        pylab.show()

    if savefile:
        pylab.savefig(savefile)

    return figure

 def test_Intensity_1(self):
     """Test a case of distributed intensity values."""
     # Create label image with only one region
     label_image = scipy.zeros(2*2*2, dtype=scipy.int8).reshape(2,2,2)
     
     # Create original image with two equally distibuted intensity value
     original_image = scipy.zeros(2*2*2, dtype=scipy.int8)
     original_image[:4] = -1
     original_image[4:] = 1
     original_image = original_image.reshape(2,2,2)
     
     # Initialize object
     statistics = LabelImageStatistics(label_image, original_image)
     
     # Computed expected result
     i = scipy.array([-1,-1,-1,-1,1,1,1,1])
     h = scipy.histogram(i, statistics._intensity_distribution_local_histogram_width)
     hr = scipy.array(h[0]) / float(h[0].sum())
     g = stats.norm(*stats.norm.fit(i))
     r = abs(hr - g.pdf(h[1][:-1]))
     r *= h[1][-2] - h[1][0]
     r = r.sum()
     
     # Check created intensity distribution
     intensity_distributions = statistics.get_intensity_distributions()
     self.assertEqual(len(intensity_distributions), 1)
     self.assertEqual(intensity_distributions[0], i.std())
     
     intensity_distribution_histogram = statistics.get_intensity_distribution_histogram()
     self.assertEqual(intensity_distribution_histogram[0][statistics.get_intensity_distribution_histogram_width()/2], 1)
     self.assertEqual(intensity_distribution_histogram[0].max(), 1)
     self.assertEqual(intensity_distribution_histogram[0].min(), 0)
     self.assertEqual(intensity_distribution_histogram[1].mean(), i.std())

def include_image_level_features(orig_imga,fvector_l,featsel):
    # include grayscale values ?
    f_input_gray = featsel['input_gray']
    if f_input_gray is not None:
        shape = f_input_gray
        #print orig_imga.shape
        fvector_l += [sp.misc.imresize(colorconv.gray_convert(orig_imga), shape).ravel()]

    # include color histograms ?
    f_input_colorhists = featsel['input_colorhists']
    if f_input_colorhists is not None:
        nbins = f_input_colorhists
        colorhists = sp.empty((3,nbins), 'f')
        if orig_imga.ndim == 3:
            for d in xrange(3):
                h = sp.histogram(orig_imga[:,:,d].ravel(),
                                 bins=nbins,
                                 range=[0,255])
                binvals = h[0].astype('f')
                colorhists[d] = binvals
        else:
            raise ValueError, "orig_imga.ndim == 3"
            #h = sp.histogram(orig_imga[:,:].ravel(),
            #                 bins=nbins,
            #                 range=[0,255])
            #binvals = h[0].astype('f')
            #colorhists[:] = binvals

        #feat_l += [colorhists.ravel()]
        fvector_l += [colorhists.ravel()]

    return fvector_l

def cluster_colors(image_url, num_clusters=5):
    """
    Return the most clustered colors of an image.
    Use scipy's k-means clustering algorithm.
    """

    print 'Reading image...'
    response = requests.get(image_url)
    im = Image.open(StringIO(response.content))
    im = im.resize((150, 150))      # optional, to reduce time
    ar = scipy.misc.fromimage(im)
    shape = ar.shape
    ar = ar.reshape(scipy.product(shape[:2]), shape[2])
    ar = ar.astype(float)

    print 'Finding clusters...'
    # k-means clustering
    codes, dist = scipy.cluster.vq.kmeans(ar, num_clusters)
    vecs, dist = scipy.cluster.vq.vq(ar, codes)         # assign codes
    counts, bins = scipy.histogram(vecs, len(codes))    # count occurrences
    sorted_index = sorted(range(len(counts)), key=lambda index: counts[index], reverse=True)

    most_common_colors = []
    for index in sorted_index:
        peak = codes[index]
        peak = peak.astype(int)
        colour = ''.join(format(c, '02x') for c in peak)
        most_common_colors.append('#' + colour)
    return most_common_colors

def create_histogram(parameter_name, nbins=100, writeFile=True, skipfirst=0, truncate=False, smooth=False):
	"""
	Returns a histogram and some statistics about this parameter.
		
	@param writeFile: if true, write the histogram to paramname.histogram
	"""
	f = "%s-chain-0.prob.dump" % parameter_name
	values = numpy.recfromtxt(f)[skipfirst::nevery]

	statistics = {
		'min':   float(values.min()),
		'max':   float(values.max()),
		'stdev': float(values.std()),
		'mean':  float(values.mean()),
		'median':float(numpy.median(values)),
		'q1':    float(scipy.stats.scoreatpercentile(values, 25)),
		'q3':    float(scipy.stats.scoreatpercentile(values, 75)),
		'p5':    float(scipy.stats.scoreatpercentile(values, 5)),
		'p95':    float(scipy.stats.scoreatpercentile(values, 95)),
	}
	
	hist = scipy.histogram(values, bins = nbins if not smooth else nbins*10, normed=True)
	histwithborders = numpy.dstack([hist[1][0:nbins], hist[1][1:nbins+1], hist[0]])
	if writeFile:
		scipy.savetxt('%s.histogram' % parameter_name, histwithborders[0], delimiter="\t")
	return histwithborders[0], statistics

def genMCHistogramsOpenCL(distribution, range, iterations=1000, numBins=1000):
    tester = clsim.I3CLSimRandomDistributionTester(device=openCLDevice,
                                                   workgroupSize=workgroupSize,
                                                   workItemsPerIteration=workItemsPerIteration,
                                                   randomService=rng,
                                                   randomDistribution=distribution)
    
    values = tester.GenerateRandomNumbers(iterations)
    samples = len(values)
    print("generated")
    
    values = numpy.array(values)/I3Units.nanometer # convert to numpy array and convert units
    print("converted")
    
    range_width=range[1]-range[0]
    
    num_orig, bins = scipy.histogram(values, range=range, bins=numBins)
    print("hist1 complete")
    
    del values # not needed anymore
    print("deleted")
    
    num=[]
    for number in num_orig:
        num.append(float(number)/float(samples)/float(range_width/float(numBins)))
    num=numpy.array(num)
    
    bins = numpy.array(bins[:-1])+(bins[1]-bins[0])/2.
    
    return dict(num=num, bins=bins)

 def _computeEntropy(self):
     ''' Compute the entropy of the histogram of distances from the pivot element. Low entropy
     scores means the distribution is concentrated, and thus may be a good candidate for splitting.
     High entropy (at limit, a uniform distribution), may indicate that there is no good separation
     of the elements in this node.
     '''
     assert(self.Pivot != None)
     assert(self.Ds != None)
     assert(len(self.Ds) >= self.Tau)
     assert(len(self.Ds) == len(self.items))
         
     #create a list of distances not including the sample which was selected as the pivot
     #...which will have a distance of zero, within numerical errors.
     Dx = [D for D in self.Ds if D>0.01]
     
     #compute histogram using 10 bins of the Dx list
     HistInfo = scipy.histogram(Dx, bins=10)
     pk = scipy.array( HistInfo[0] )
     epsilon = 0.000001
     H = entropy(pk+epsilon)  #avoids log0 warnings
     #print "Histogram: ", HistInfo[0]
     #print "Entropy: %f"%H
     #print "Range: Min(>0)=%f, Max=%f, Mean=%f, Median=%f"%(min(Dx),max(self.Ds),scipy.mean(self.Ds),scipy.median(self.Ds))
     
     return H
     

def getPredominantColor(filename):
    im = Image.open(filename).convert('RGB')

    # Convert to numpy array
    ar = scipy.misc.fromimage(im)

    # Get dimensions
    shape = ar.shape

    # Convert to bidimensional array of width x height rows and 3 columns (RGB)
    ar = ar.reshape(scipy.product(shape[:2]), shape[2])

    # Find cluster centers and their distortions
    # codes contains the RGB value of the centers
    codes, dist = scipy.cluster.vq.kmeans(ar.astype(float), NUM_CLUSTERS)

    # Maps all the pixels in the image to their respective centers
    vecs, dist = scipy.cluster.vq.vq(ar, codes)

    # Counts the occurances of each color (NUM_CLUSTER different colors after the mapping)
    counts, bins = scipy.histogram(vecs, len(codes))

    # Find most frequent color
    index_max = scipy.argmax(counts)
    peak = codes[index_max]

    return peak.astype(int)

def getImageDescriptor(model, im, conf):
	im = standardizeImage(im)
	height, width = im.shape[:2]
	numWords = model.vocab.shape[1]
	frames, descrs = getPhowFeatures(im, conf.phowOpts)
	# quantize appearance
	if model.quantizer == 'vq':
		binsa, _ = vq(descrs.T, model.vocab.T)
	elif model.quantizer == 'kdtree':
		raise ValueError('quantizer kdtree not implemented')
	else:
		raise ValueError('quantizer {0} not known or understood'.format(model.quantizer))
	hist = []
	for n_spatial_bins_x, n_spatial_bins_y in zip(model.numSpatialX, model.numSpatialX):
		binsx, distsx = vq(frames[0, :], linspace(0, width, n_spatial_bins_x))
		binsy, distsy = vq(frames[1, :], linspace(0, height, n_spatial_bins_y))
		# binsx and binsy list to what spatial bin each feature point belongs to
		if (numpy.any(distsx < 0)) | (numpy.any(distsx > (width/n_spatial_bins_x+0.5))):
			print ("something went wrong")
			import pdb; pdb.set_trace()
		if (numpy.any(distsy < 0)) | (numpy.any(distsy > (height/n_spatial_bins_y+0.5))):
			print ("something went wrong")
			import pdb; pdb.set_trace()
		# combined quantization
		number_of_bins = n_spatial_bins_x * n_spatial_bins_y * numWords
		temp = arange(number_of_bins)
		# update using this: http://stackoverflow.com/questions/15230179/how-to-get-the-linear-index-for-a-numpy-array-sub2ind
		temp = temp.reshape([n_spatial_bins_x, n_spatial_bins_y, numWords])
		bin_comb = temp[binsx, binsy, binsa]
		hist_temp, _ = histogram(bin_comb, bins=range(number_of_bins+1), density=True)
		hist.append(hist_temp)
	
	hist = hstack(hist)
	hist = array(hist, 'float32') / sum(hist)
	return hist

def reduce_colors(image, k):
    '''Apply kmeans algorithm.
        Input:   image, number of clusters to use
        Returns: colors, 
                 counts per color, 
                 new image
    '''
    if k > 32:
        print "Setting colors to maximum allowed of 32"
        k = 32
    rows, cols, rgb = image.shape
    # reshape the image in a single row array of RGB pixels
    image_row = np.reshape(image,(rows * cols, 3))
    #HERE ADD CODE TO GET A GOOD GUESS OF COLORS AND PASS THAT AS
    #SECOND ARGUMENT TO kmeans
    #image_array_sample = shuffle(image_row, random_state=0)[:1000]
    #kguess = kmeans(image_array_sample, k)
    #colors,_ = kmeans(image_row, kguess)
    # perform the clustering
    colors,_ = kmeans(image_row, k)
    # vector quantization, assign to each pixel the index of the nearest centroid (i=1..k)
    qnt,_ = vq(image_row,colors)
    # reshape the qnt vector to the original image shape
    image_centers_id = np.reshape(qnt,(rows, cols))
    # assign the color value to each pixel
    newimage = colors[image_centers_id]
    #count number of pixels of each cluster color
    counts,bins = sp.histogram(qnt, len(colors))
    return colors, counts, newimage

def get_dominant_color(image_path):
    '''
    Parse image and return dominant color in image.

    @param image_path: Image path to parse.
    @return: Return dominant color, format as hexadecimal number. 
    '''
    # print 'reading image'
    im = Image.open(image_path)
    im = im.resize((150, 150))      # optional, to reduce time
    ar = scipy.misc.fromimage(im)
    shape = ar.shape
    ar = ar.reshape(scipy.product(shape[:2]), shape[2])
    
    # print 'finding clusters'
    NUM_CLUSTERS = 5
    codes, dist = scipy.cluster.vq.kmeans(ar, NUM_CLUSTERS)
    # print 'cluster centres:\n', codes
    
    vecs, dist = scipy.cluster.vq.vq(ar, codes)         # assign codes
    counts, bins = scipy.histogram(vecs, len(codes))    # count occurrences
    
    index_max = scipy.argmax(counts)                    # find most frequent
    peak = codes[index_max]
    colour = ''.join(chr(c) for c in peak).encode('hex')
    # print 'most frequent is %s (#%s)' % (peak, colour)
    
    return "#%s" % (colour[0:6])

def genMCHistogramsHost(distribution, hist_range, iterations=10000000, numBins=1000):
    print("generating (host)")

    values = []
    for i in range(iterations):
        values.append(distribution.SampleFromDistribution(rng, []))
    samples = len(values)
    print("generated (host)")
    
    values = numpy.array(values)/I3Units.nanometer # convert to numpy array and convert units
    print("converted (host)")
    
    range_width=hist_range[1]-hist_range[0]
    
    num_orig, bins = scipy.histogram(values, range=hist_range, bins=numBins)
    print("hist1 complete (host)")
    
    del values # not needed anymore
    print("deleted (host)")
    
    num=[]
    for number in num_orig:
        num.append(float(number)/float(samples)/float(range_width/float(numBins)))
    num=numpy.array(num)
    
    bins = numpy.array(bins[:-1])+(bins[1]-bins[0])/2.
    
    return dict(num=num, bins=bins)

def tag_images_with_color_value(NUM_CLUSTERS = 4, INPUT_FOLDER = './data/covers/'):

    isbn = list()
    cover_color = list()

    files = os.listdir(INPUT_FOLDER)
    for eachFile in files:
        print eachFile
        im = Image.open(INPUT_FOLDER + eachFile)
        im = im.resize((50, 50))                          # optional, to reduce time
        ar = scipy.misc.fromimage(im)
        shape = ar.shape
        print len(shape)

        if len(shape) == 2:
            ar = ar.reshape(scipy.product(shape[:1]), shape[1])
        else:
            ar = ar.reshape(scipy.product(shape[:2]), shape[2])

        # finding clusters
        codes, dist = scipy.cluster.vq.kmeans(ar, NUM_CLUSTERS)
        # cluster centres:\n', codes

        vecs, dist = scipy.cluster.vq.vq(ar, codes)         # assign codes
        counts, bins = scipy.histogram(vecs, len(codes))    # count occurrences

        index_max = scipy.argmax(counts)                    # find most frequent
        peak = codes[index_max]
        colour = ''.join(chr(c) for c in peak).encode('hex')

        isbn.append(eachFile[:-4])
        cover_color.append(colour)

    result = zip(isbn, cover_color)
    return result

def degree_distrib(net, deg_type="total", node_list=None, use_weights=True,
                   log=False, num_bins=30):
    '''
    Computing the degree distribution of a network.
    
    Parameters
    ----------
    net : :class:`~nngt.Graph` or subclass
        the network to analyze.
    deg_type : string, optional (default: "total")
        type of degree to consider ("in", "out", or "total").
    node_list : list or numpy.array of ints, optional (default: None)
        Restrict the distribution to a set of nodes (default: all nodes).
    use_weights : bool, optional (default: True)
        use weighted degrees (do not take the sign into account: all weights
        are positive).
    log : bool, optional (default: False)
        use log-spaced bins.
    
    Returns
    -------
    counts : :class:`numpy.array`
        number of nodes in each bin
    deg : :class:`numpy.array`
        bins
    '''
    ia_node_deg = net.get_degrees(node_list, deg_type, use_weights)
    ra_bins = sp.linspace(ia_node_deg.min(), ia_node_deg.max(), num_bins)
    if log:
        ra_bins = sp.logspace(sp.log10(sp.maximum(ia_node_deg.min(),1)),
                               sp.log10(ia_node_deg.max()), num_bins)
    counts,deg = sp.histogram(ia_node_deg, ra_bins)
    ia_indices = sp.argwhere(counts)
    return counts[ia_indices], deg[ia_indices]

def getDomIMAGEColor( imName ):
	# Reference:
	# 	http://stackoverflow.com/questions/3241929/
	# 	python-find-dominant-most-common-color-in-an-image

	# number of k-means clusters
	NUM_CLUSTERS = 4

	# Open target image
	im = imName
	im = im.resize((150, 150))      # optional, to reduce time
	ar = scipy.misc.fromimage(im)
	shape = ar.shape
	ar = ar.reshape(scipy.product(shape[:2]), shape[2])
	ar = ar.astype(float)

	# Find clusters
	codes, dist = scipy.cluster.vq.kmeans(ar, NUM_CLUSTERS)
	vecs, dist = scipy.cluster.vq.vq(ar, codes)         # assign codes
	counts, bins = scipy.histogram(vecs, len(codes))    # count occurrences

	# Find most frequent
	index_max = scipy.argmax(counts)                    
	peak = codes[index_max]
	color = ''.join(chr(int(c)) for c in peak).encode('hex')

	return (peak, color)

    def determine_dominant_color_in_image(self, image):
        NUM_CLUSTERS = 5
            
        ar = scipy.misc.fromimage(image)
        shape = ar.shape
        if len(shape) > 2:
            ar = ar.reshape(scipy.product(shape[:2]), shape[2])

        codes, _ = scipy.cluster.vq.kmeans(ar, NUM_CLUSTERS)
        # print "Before: %s" % codes
        original_codes = codes
        for low, hi in [(60, 200), (35, 230), (10, 250)]:
            codes = scipy.array([code for code in codes 
                                 if not ((code[0] < low and code[1] < low and code[2] < low) or
                                         (code[0] > hi and code[1] > hi and code[2] > hi))])
            if not len(codes): codes = original_codes
            else: break
        # print "After: %s" % codes
    
        vecs, _ = scipy.cluster.vq.vq(ar, codes)         # assign codes
        counts, bins = scipy.histogram(vecs, len(codes))    # count occurrences
        # colors = [''.join(chr(c) for c in code).encode('hex') for code in codes]
        # total = scipy.sum(counts)
        # print dict(zip(colors, [count/float(total) for count in counts]))
        index_max = scipy.argmax(counts)                    # find most frequent
        peak = codes[index_max]
        color = ''.join(chr(c) for c in peak).encode('hex')
        # print 'most frequent is %s (#%s)' % (peak, color)
        
        return color[:6]