from astropy.io import fits
 
from scipy.optimize import curve_fit
import glob
from photutils import CircularAperture
BLACK = 0


def OpenImageHist(filename):
        """Opens a FIT image and puts pixel ADU values into an array"""
        image = fits.getdata(filename)
        A = image.flatten()
        return A
        


def OpenDirectory(Directory):
    filenames = glob.glob(Directory+"*.fit")
    return filenames

def Histogram(Array, bins): 
    """Plots histogram from flattened array of data, automatic range for gaussian distribution in array."""
    h_test = plt.hist(Array, bins, (min(Array), max(Array))) #Histogram #1 over entire array
    plt.clf()
    test_y = h_test[0]
    test_x = h_test[1]
    
    PeakPosition = np.argmax(test_y)    #Find peak of histogram
    Peak = test_x[PeakPosition]          #Find position of  peak along x axis
    BottomofRange = Peak-300             #Tailor new range around peak, very vague estimate, could probably use FWHM instead? Or a percentage of peak
    TopofRange = Peak+300
    EstRange = (BottomofRange, TopofRange)#Estimate 1 of range
    
    Test2_h = plt.hist(Array, bins, EstRange) #Histogram #2 has slightly better range, still not great
    plt.clf()
    Test2_y = Test2_h[0]                 
    Test2_x = Test2_h[1]
    
    Counter = 0                           #loop over the y values (Counts in histogram) to find the first non-zero value.                                                  
    for i in Test2_y:
        if i == 0:
            Counter+=1
        if i >0:
            break
    Test2RangeLower = Test2_x[Counter]      #The next test range is between all values that are non-zero, 
                                            #this should be changed to a number dependent on the values in the histogram,
                                            #because they might not all be non-zero.
      
    Counter = 0                             #loop over reversed list, to cut off any non-zero values.
    for i in reversed(Test2_y):
        if i == 0:
            Counter+=1
        if i >0:
            break
    Test2RangeUpper = Test2_x[len(Test2_y) - Counter]
    
      
    Test2Range = (Test2RangeLower, Test2RangeUpper)
    
    Test3_h = plt.hist(Array, bins, Test2Range)
    plt.clf()
    Test3_x = Test3_h[1]
    Test3_y = Test3_h[0]
    
    Counter = 0 
    for i in Test3_y:
        if i<2000:                      #Instead of value here could have a percentage of peak count?
            Counter +=1
        if i>=2000:
            break
    NewRangeLower = int(round(Test3_x[Counter]))
    Counter = 0
    for i in reversed(Test3_y):
        if i<2000:
            Counter +=1
        if i>=2000:
            break
    NewRangeHigher = int(round(Test3_x[len(Test3_y) - Counter]))
    FinalRange = (NewRangeLower, NewRangeHigher)
    Difference = NewRangeHigher-NewRangeLower 
    bins = (Difference/2)-7
    h = plt.hist(Array, bins, FinalRange)
    xedge = h[1]
    y = h[0]        
    xcentre = (xedge + np.roll(xedge, -1))/2.0
    x = xcentre[0:len(y)]
    plt.clf()
    return x, y, NewRangeLower, NewRangeHigher
    
def GaussFunc(x,a,x0,sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
    
def FitAGaussian(Array, bins): #Not Working
    """Fits a Gaussian Distribution over histogram"""
    
    h = Histogram(Array, bins)
    x = h[0]
    y = h[1]
    
    n = len(x)
    mean = sum(x)/n 
    
    popt, pcov = curve_fit(GaussFunc, x, y, (1, mean, 1),absolute_sigma=True)
 
    return popt, pcov


def FindExposureTime(path):
    head, tail = ntpath.split(path)
    B = tail
    #B = ntpath.basename(head)
    ExposureTime = float(B[:3])
    return ExposureTime



def Bins(array):
    binwidth = 10
    bins = np.arange(0, max(array) + binwidth, binwidth)
    return bins

def OpenImage(FileName):
    """Opens a FIT image and puts pixel ADU values into an array"""
    image = fits.getdata(FileName)
    return image

def OpenDirectory(Directory):
    filenames = glob.glob(Directory+"*.fit")
    return filenames


def count(pic):
    xsize, ysize = np.shape(pic)
    temp = pic
    clust = []
    
    for x in range(xsize):
        for y in range(ysize):
            if temp[x,y] >= Threshold and temp[x,y] != WHITE:
                clust.append((fill(temp, xsize, ysize, x, y)))
                

    plt.imshow(temp, cmap='gray_r', origin='lower',interpolation='none', vmin=WHITE - 1,vmax=WHITE)         
    
    return clust

def fill(pic, xsize, ysize, x_init, y_init):
    thisClustList = []
    queue = [(x_init, y_init)]

    while len(queue) != 0:
        top = queue.pop()
        x, y = top[0], top[1]
        pval = pic[x, y]
        if pval >= Threshold and pval != WHITE:
            thisClustList.append([x, y])
            pic[x, y] = WHITE
            
            if x > 0: queue.append((x-1, y))
            if x < (xsize - 1): queue.append((x+1, y))
            if y > 0: queue.append((x, y-1))
            if y < (ysize - 1): queue.append((x, y+1))


    return thisClustList

def blob_centre(blob, pic):
    blobxlist = []
    blobylist = []
    for pixel in blob:
         blobxlist.append(pixel[1])
         blobylist.append(pixel[0])
    nsum = 0
    nxsum = 0
    nysum = 0
    for pixel in blob:
        nsum += pic[pixel[0]][pixel[1]]
        nxsum += (pic[pixel[0]][pixel[1]] * pixel[1])
        nysum += (pic[pixel[0]][pixel[1]] * pixel[0])

    xcen = nxsum/nsum
    ycen = nysum/nsum
    return xcen, ycen

#path = "C:\\Users\\Blake\\Dropbox\\Msci_project\\SimulatedStars\\50StarsWithStotFluctuation.fits"
path = "C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei_10.fit" 
#
picHist = OpenImageHist(path)
picFill = OpenImage(path)
pic = OpenImage(path)
WHITE = max(picHist) + 1
bins = Bins(picHist)
popt, pcov = FitAGaussian(picHist, bins)
mean = popt[1]
std = abs(popt[2])
Threshold = mean+std*6
#Threshold = picHist.mean() + picHist.std()

plt.figure(1)

plt.imshow(pic, cmap='gray', origin='lower',interpolation='none', vmin=np.median(pic),vmax=np.median(pic)*2.2)
plt.colorbar()
plt.figure(2)
clist = count(picFill)



print len(clist)


bloblist = []
# i = 0 
# for blob in clist:
#         print i
#         print len(blob)
#         i+=1
        
blob = clist[179]

for pixel in blob:
     bloblist.append(pic[pixel[0]][pixel[1]])
plt.figure(3)
print len(bloblist)
plt.hist(bloblist, bins = np.arange(0, max(bloblist) + 10,10))
plt.axis([200, 4000, 0, 120])
# maxb = max(bloblist)
# bloblist.index(maxb)
# centreCo = blob[]
# apertures = CircularAperture(centreCo, r=50.)
# apertures.plot(color='blue', lw=1.5, alpha=0.5)
plt.show()


<\verbatim>



-- Public.LunaBorella - 27 Oct 2017