# -*- coding: utf-8 -*-
"""
"""

import cv2
import numpy as np
from matplotlib import pyplot as plt
from preproces_max import *
from bisect import bisect
"""
From :: https://code.google.com/p/agpy/source/browse/trunk/agpy/radialprofile.py?r=317img.shape

"""
def azimuthalAverage(image, center=None, stddev=False, returnradii=False, return_nr=False, 
        binsize=0.5, weights=None, steps=False, interpnan=False, left=None, right=None):
    """
    Calculate the azimuthally averaged radial profile.

    image - The 2D image
    center - The [x,y] pixel coordinates used as the center. The default is 
             None, which then uses the center of the image (including 
             fractional pixels).
    stddev - if specified, return the azimuthal standard deviation instead of the average
    returnradii - if specified, return (radii_array,radial_profile)
    return_nr   - if specified, return number of pixels per radius *and* radius
    binsize - size of the averaging bin.  Can lead to strange results if
        non-binsize factors are used to specify the center and the binsize is
        too large
    weights - can do a weighted average instead of a simple average if this keyword parameter
        is set.  weights.shape must = image.shape.  weighted stddev is undefined, so don't
        set weights and stddev.
    steps - if specified, will return a double-length bin array and radial
        profile so you can plot a step-form radial profile (which more accurately
        represents what's going on)
    interpnan - Interpolate over NAN values, i.e. bins where there is no data?
        left,right - passed to interpnan; they set the extrapolated values

    If a bin contains NO DATA, it will have a NAN value because of the
    divide-by-sum-of-weights component.  I think this is a useful way to denote
    lack of data, but users let me know if an alternative is prefered...
    
    """
    # Calculate the indices from the image
    y, x = np.indices(image.shape)

    if center is None:
        center = np.array([(x.max()-x.min())/2.0, (y.max()-y.min())/2.0])

    r = np.hypot(x - center[0], y - center[1])

    if weights is None:
        weights = np.ones(image.shape)
    elif stddev:
        raise ValueError("Weighted standard deviation is not defined.")

    # the 'bins' as initially defined are lower/upper bounds for each bin
    # so that values will be in [lower,upper)  
    nbins = (np.round(r.max() / binsize)+1)
    maxbin = nbins * binsize
    bins = np.linspace(0,maxbin,nbins+1)
    # but we're probably more interested in the bin centers than their left or right sides...
    bin_centers = (bins[1:]+bins[:-1])/2.0

    # Find out which radial bin each point in the map belongs to
    whichbin = np.digitize(r.flat,bins)

    # how many per bin (i.e., histogram)?
    # there are never any in bin 0, because the lowest index returned by digitize is 1
    nr = np.bincount(whichbin)[1:]

    # recall that bins are from 1 to nbins (which is expressed in array terms by arange(nbins)+1 or xrange(1,nbins+1) )
    # radial_prof.shape = bin_centers.shape
    if stddev:
        radial_prof = np.array([image.flat[whichbin==b].std() for b in xrange(1,nbins+1)])
    else:
        radial_prof = np.array([(image*weights).flat[whichbin==b].sum() / weights.flat[whichbin==b].sum() for b in xrange(1,nbins+1)])

    #import pdb; pdb.set_trace()

    if interpnan:
        radial_prof = np.interp(bin_centers,bin_centers[radial_prof==radial_prof],radial_prof[radial_prof==radial_prof],left=left,right=right)

    if steps:
        xarr = np.array(zip(bins[:-1],bins[1:])).ravel() 
        yarr = np.array(zip(radial_prof,radial_prof)).ravel() 
        return xarr,yarr
    elif returnradii: 
        return bin_centers,radial_prof
    elif return_nr:
        return nr,bin_centers,radial_prof
    else:
        return radial_prof

def get_center(img, per):
    tt = int((per/100.0)*max_img(img))
    print tt
    ret,thresh = cv2.threshold(img,tt,255,cv2.THRESH_BINARY | cv2.THRESH_OTSU)
    contours,hierarchy = cv2.findContours(thresh, 1, 2)
    #contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
    #Max Area contour
    cnt,mm=0,0
    for c in contours:
        if mm < cv2.contourArea(c):
            mm = cv2.contourArea(c)
            cnt = c
    
    ellipse = cv2.fitEllipse(cnt)
    return ((ellipse[0][0]+ellipse[1][0])/2,(ellipse[0][1]+ellipse[1][1])/2)


def radial_profile(data, center):
    y, x = np.indices((data.shape))
    r = np.sqrt((x - center[0])**2 + (y - center[1])**2)
    r = r.astype(np.int)

    tbin = np.bincount(r.ravel(), data.ravel())
    nr = np.bincount(r.ravel())
    radialprofile = tbin / nr
    return radialprofile
    

img_num = 104176
per = 60#%

img = cv2.imread('images_test_rev1/'+str(img_num)+'.jpg')
img = img[100:300, 100:300,:]
img = cv2.resize(img, (100, 100))

plt.subplot(121), plt.imshow(img)
plt.tight_layout()
plt.axis('off')
plt.title('Image')

img = cv2.imread('images_test_rev1/'+str(img_num)+'.jpg', 0)
img = img[100:300, 100:300]
img = cv2.resize(img, (100, 100))


rad = radial_profile(img,get_center(img, per))
ccc = np.cumsum(rad)
#plt.subplot(1,3,2), plt.plot(rad)
plt.subplot(122), plt.plot(ccc)
b50 = bisect(ccc, ccc[-1]/2)
plt.plot([b50,b50],[0,ccc[b50]],  'k-', lw=2, label='50%')
plt.plot([0,b50],[ccc[b50],ccc[b50]], 'k-', lw=2)
plt.legend(bbox_to_anchor=(0, 0, 1, 0))
#plt.plot([b50,ccc[b50]], [0, ccc[b50]])

plt.title('Cumulative light Intensity')

plt.ylabel('Cumulative Intensity')
plt.xlabel('Radius')
plt.tight_layout()