from astropy.io import fits as pyfits
import numpy as np
import matplotlib.pyplot as plt 

image = pyfits.getdata("FinderTest2.fits")

black=0
white=2**16

def Test(image):
    s=Stats(image)
    m=PixelMaxCoord(image)
    c=Count(image)
    print s
    print m
    print c
    
def Stats(image):
    imgsum = np.sum(image)
    imgmean = np.mean(image)
    imgmedian = np.median(image)
    imgstd = np.std(image)
    imgmaxval = image.max()
    return [imgsum,imgmean,imgmedian,imgstd,imgmaxval]
 
def PixelMaxCoord(image):       
    image[0,0] #upper left corner
    xsize,ysize=image.shape
    bx,by=0,0
    for x in range(xsize):
        for y in range(ysize):
            value=image[x,y]
            if value>image.max():
                bx,by=x,y
    return (x,y) #returns coordinates of brightest pixel

def Count(image):
    xsize,ysize=image.shape
    thresh=1*(image>300) #set threshold so be 300 (out of possible 2**16), if value is above this it is assigned 1, if below, assigned 0
    #plt.imshow(thresh)
    #plt.show()
    clustList = [] #list of clusters of pixels (stars)
    for x in range(xsize):
        for y in range(ysize):
            #print x,y,thresh[x,y]
            clust=[] #list of pixels in star
            if thresh[x,y]==1: #if pixel value is 1 then fill function called to start finding cluster coordinates
                Fill(thresh,xsize,ysize,x,y,clust) #call fill to begin counting
                clustList.append(clust) #add each cluster found to list of all clusters
                StarCount=len(clustList) #number of stars in file
    #return clustList   
    n1=0
    ExpTime=0.5 #(s), varies
    Gain=2.39 #check this against telescope and twiki
    FluxList=[] #list of fluxes for each star
    while n1<StarCount:
         Counts=0
         n2=0
         while n2<len(clustList[n1]): #while within the list of coordinates for star
             Coords=clustList[n1] #list of coordinates in star
             Coords2=Coords[n2] 
             Counts=Counts+image[Coords2[0],Coords2[1]] #pixel values (counts) in star
             n2=n2+1
         CountsMean=Counts/len(clustList[n1]) #mean value of counts of each pixel in the star
         Flux=(Gain*CountsMean)/ExpTime
         FluxList.append(Flux) #add each flux to list of fluxes for all stars
         n1=n1+1
    return StarCount,FluxList #returns number of stars and the flux of each
              
def Fill(image,xsize,ysize,x,y,clust):
    #print x,y,image[x,y]
    if image[x,y] == 0:
       return #if pixel value is zero then ignore
    else: 
        image[x,y]=0 #if pixel value isn't zero then we set it to zero, proof in the image that all stars found when image is blank
        clust.append([x,y])
    if x>0: Fill(image,xsize,ysize,x-1,y,clust)
    if x<(xsize-1): Fill(image,xsize,ysize,x+1,y,clust)
    if y>0: Fill(image,xsize,ysize,x,y-1,clust)
    if y<(ysize-1): Fill(image,xsize,ysize,x,y+1,clust)