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
test_x = h_test

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
Test2_x = Test2_h

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
Test3_y = Test3_h

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
y = h
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
y = h

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):
B = tail
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, top
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)
blobylist.append(pixel)
nsum = 0
nxsum = 0
nysum = 0
for pixel in blob:
nsum += pic[pixel][pixel]
nxsum += (pic[pixel][pixel] * pixel)
nysum += (pic[pixel][pixel] * pixel)

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
std = abs(popt)
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

for pixel in blob:
bloblist.append(pic[pixel][pixel])
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

Physics WebpagesRHUL WebpagesCampus Connect • Royal Holloway, University of London, Egham, Surrey TW20 0EX; Tel/Fax +44 (0)1784 434455/437520

Topic revision: r3 - 23 Nov 2017 - LunaMarieBorella1 Home   Public Web  P P P P P View Edit        Copyright © 2008-2021 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding RHUL Physics Department TWiki? Send feedback