Commit caad58fb authored by Markus Gaug's avatar Markus Gaug
Browse files

added changes suggested in merge request discussion

parent 623670aa
Pipeline #6012 passed with stages
in 2 minutes and 51 seconds
......@@ -23,50 +23,30 @@ CONFIGURATION:
6) Write the conversion into the parameter "conv"
7) Run the program
"""
from scipy import optimize as opt
import astropy
from scipy.optimize import leastsq
from scipy import interpolate, arange, zeros
from astropy.io import fits
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import copy
from scipy import *
from matplotlib import cm
import math
from pylab import figure, axes, pie, title, show
import sys
#
# user configuration:
#
# image file in FITS format (green):
file='/Users/markusgaug/CTA/LIDAR/Measurements/20170410_Laser/IMG_0767-G.fits'
#file='./_DSC0530-G.fits'
#file='/Users/markusgaug/CTA/LIDAR/Measurements/20170410_Laser/IMG_0767-G.fits'
# image cutout (from gimp) coordinates
rangey1=600
rangey2=1650
rangex1=1000
rangex2=1800
#rangex1=2000
#rangex2=2300
#rangey1=760
#rangey2=1100
# conversion from points to mm (please use GIMP and the reference grid on the picture!
conv = 1./8.31 # conversion from points to mm
#conv = 0.844921 # 1./8.31 # conversion from points to mm
# --------------------------------------------------------------------------------
#file='/Users/markusgaug/CTA/LIDAR/Measurements/20170410_Laser/IMG_0826-G.fits'
#rangey1=1010
#rangey2=1100
#rangex1=1950
#rangex2=2050
# --------------------------------------------------------------------------------
beamdiameterdef = 1./math.exp(1.)**2 # definition of beam diameter
securityfac = 9 # times the fitted half widths beyond which pixels are counted as background
idx = rangex2-rangex1
idy = rangey2-rangey1
#rangex1=900
#rangex2=1850
#rangey1=1250
#rangey2=2700
#conv = 0.0833 # conversion from points to mm
# define model function and pass independant variables x and y as a list
def gaussian(height, center_x, center_y, width_x, width_y, theta, offset):
"""Returns a gaussian function with the given parameters"""
width_x = float(width_x)
......@@ -78,24 +58,31 @@ def gaussian(height, center_x, center_y, width_x, width_y, theta, offset):
# print (height, center_x, center_y, width_x, width_y, theta, offset)
if (height<0):
return lambda x, y: 99999999.
# raise ValueError("height = %f must be positive!" % height)
if (center_x < 0.):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("center_x = %f must be positive!" % center_x)
if (center_y < 0.):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("center_y = %f must be positive!" % center_y)
if (width_x < 0.):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("width_x = %f must be positive!" % width_x)
if (width_y < 0.):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("width_y = %f must be positive!" % width_y)
if (offset < 0.):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("offset = %f must be positive!" % offset)
if (theta > 3.*np.pi):
return lambda x, y: 99999999.
return lambda x, y: 99999999.
# raise ValueError("theta = %f must be smaller than 3pi!" % theta)
return lambda x, y: offset + height * np.exp( - (a* (center_x - x) ** 2 + c * (center_y - y) ** 2 + 2*b * (center_x - x)*(center_y-y)) / 2.)
#-(((center_x - x) / width_x) ** 2 + ((center_y - y) / width_y) ** 2) / 2)
def moments(data):
"""Returns (height, x, y, width_x, width_y)
"""Returns (height, x, y, width_x, width_y, theta, offset)
the gaussian parameters of a 2D distribution by calculating its
moments """
total = data.sum()
......@@ -107,9 +94,9 @@ def moments(data):
row = data[int(x), :]
width_y = np.sqrt(np.abs((np.arange(row.size) - x) ** 2 * row).sum() / row.sum())
height = data.max()
offset = total / (450*450)
offset = total / (col.size*row.size)
# print height, x, y, width_x, width_y, -45*np.pi/180, offset
print('from moments: height=', height, 'x=', x, 'y=', y, 'width_x=', width_x, 'width_y=', width_y, 'theta=', -45*np.pi/180, 'offset=', offset)
return height, x, y, 2.*width_x, width_y, -45.*np.pi/180, offset
......@@ -119,181 +106,260 @@ def fitgaussian(data):
params = moments(data)
errorfunction = lambda p: np.ravel(gaussian(*p)(*np.indices(data.shape)) -
data)
p, success = opt.leastsq(errorfunction, params)
p, success = leastsq(errorfunction, params)
return p
hdulist = fits.open(file)
print (hdulist.info)
hdu = hdulist[0]
print (hdu)
#scidata_red = hdu.data[0]
#scidata_green = hdu.data[1]
# scidata_blue = hdu.data[2]
# print (scidata_red[0])
# print (scidata_red.shape)
# scidata_sl = scidata_green[1225:1425,1040:1235]
colors = ('red', 'green', 'blue')
titles = ('full', 'background-subtracted', 'blue')
f, axarr = plt.subplots(1, 3)
x = 0.
y = 0.
width_x = 0.
width_y = 0.
for i in range(0, 1):
def run_script(file,rangey1,rangey2,rangex1,rangex2):
print('ranges: ',rangey1,rangey2,rangex1,rangex2)
outfile = str(file).replace(".fits","")
outfile1 = outfile+"_1.pdf"
print ("outfile=",outfile1)
beamdiameterdef = 1./math.e**2 # definition of beam diameter
securityfac = 5 # times the fitted half widths beyond which pixels are counted as background
fac = 0.025 # required precision in units ellipse half width and length
upx = 100 # upper range of plot (in mm)
idx = rangex2-rangex1
idy = rangey2-rangey1
dim = int(0.5*(idx+idy)*4)
hdulist = fits.open(file)
print (hdulist.info)
hdu = hdulist[0]
print (hdu)
print (hdu.shape)
colors = ('red', 'red', 'red', 'black')
titles = ('linear, fitted', 'log-scale', 'background-subtracted', 'beam diameter and off-region')
f, axarr = plt.subplots(2, 2)
x = 0.
y = 0.
width_x = 0.
width_y = 0.
scidata = hdu.data
scidata_sl = scidata[rangey1:rangey2,rangex1:rangex2]
# print scidata_sl
axarr[i].imshow(scidata_sl, origin='lower', cmap='gnuplot')
axarr[i].set_title(titles[i], color=colors[i])
# print scidata_sl
axarr[0, 0].imshow(scidata_sl, origin='lower', cmap='gnuplot')
axarr[0, 0].set_title(titles[0], color=colors[0], fontsize=10)
params = fitgaussian(scidata_sl)
print(' height, x, y, width_x, width_y, theta, offset')
print('Fitted parameters: ', [ "{:f} ".format(x) for x in params] )
fit = gaussian(*params)
axarr[i].contour(fit(*np.indices(scidata_sl.shape))) #, cmap=plt.cm.copper)
axarr[0, 0].contour(fit(*np.indices(scidata_sl.shape))) #, cmap=plt.cm.copper)
(height, x, y, width_x, width_y, theta, offset) = params
axarr[i].text(0., -0.3, '$\sigma_x$={:3.1f}, $\sigma_y$={:3.1f}, $alpha$={:.0f}$\circ$'.format(width_y, width_x, theta*57.3,prec='.1'),
verticalalignment='bottom', horizontalalignment='left',
transform=axarr[i].transAxes,
color=colors[i], fontsize=15)
num = 0
offset = 0.
scidata = hdu.data
scidata_sl = scidata[rangey1:rangey2,rangex1:rangex2]
scidata_sl2 = np.array(scidata_sl) + 0.
#print (scidata_sl2.shape)
#print (scidata_sl2.T.shape)
sum = 0.
numsum = 0
a = (np.cos(theta) ** 2) / ( securityfac* width_x ** 2) + (np.sin(theta) ** 2) / ( securityfac* width_y ** 2)
b = -(np.sin(2 * theta)) / (2 * securityfac* width_x ** 2) + (np.sin(2 * theta)) / (2 * securityfac* width_y ** 2)
c = (np.sin(theta) ** 2) / ( securityfac* width_x ** 2) + (np.cos(theta) ** 2) / ( securityfac* width_y ** 2)
for ix in range(0, idy):
for iy in range(0, idx):
# print (ix, x, width_x, iy, y, width_y, (((ix-x)/width_x)**2+((iy-y)/width_y)**2))
# if ((((ix - x) / width_x * 2) ** 2 + ((iy - y) / width_y * 2) ** 2) > 1):
if (( a* (ix - x) ** 2 + c * (iy - y) ** 2 + 2*b* (ix - x)*(iy - y)) > 1):
offset += scidata_sl2[ix, iy]
# print ('ix=',ix,'iy=',iy,'offset=',offset,'scidata=',scidata_sl2[ix,iy])
num += 1
else:
sum += scidata_sl2[ix, iy]
numsum += 1
print ('offset=', offset, ' num=', num)
if (offset != 0):
axarr[0, 0].text(0.,0.9, r'$\sigma_x$={:3.1f}, $\sigma_y$={:3.1f}, $\alpha$={:.0f}$^\circ$'.format(width_y, width_x, theta*57.3,prec='.1'),
verticalalignment='bottom', horizontalalignment='left',
transform=axarr[0, 0].transAxes,
color='white', fontsize=8)
num = 0
offset = 0.
scidata = hdu.data
scidata_sl = scidata[rangey1:rangey2,rangex1:rangex2]
scidata_sl2 = np.array(scidata_sl) + 0.
scidata_log = np.array(scidata_sl) + 0.
summ = 0.
numsum = 0
# for rotated ellipsee equations, see e.g. https://math.stackexchange.com/questions/426150/what-is-the-general-equation-of-the-ellipse-that-is-not-in-the-origin-and-rotate
a = (np.cos(theta) ** 2) / ( (securityfac* width_x) ** 2) + (np.sin(theta) ** 2) / ( (securityfac* width_y) ** 2)
b = -(np.sin(2 * theta)) / (2 * (securityfac* width_x) ** 2) + (np.sin(2 * theta)) / (2 * (securityfac* width_y) ** 2)
c = (np.sin(theta) ** 2) / ( (securityfac* width_x) ** 2) + (np.cos(theta) ** 2) / ( (securityfac* width_y) ** 2)
for ix in range(0, idy):
for iy in range(0, idx):
# print (ix, x, width_x, iy, y, width_y, (((ix-x)/width_x)**2+((iy-y)/width_y)**2))
# if ((((ix - x) / width_x * 2) ** 2 + ((iy - y) / width_y * 2) ** 2) > 1):
if (( a* (ix - x) ** 2 + c * (iy - y) ** 2 + 2*b* (ix - x)*(iy - y)) > 1):
offset += scidata_sl2[ix, iy]
# print ('ix=',ix,'iy=',iy,'offset=',offset,'scidata=',scidata_sl2[ix,iy])
num += 1
else:
summ += scidata_sl2[ix, iy]
numsum += 1
print ('sum=', offset, ' num=', num)
offset /= num
print (offset)
sum -= numsum * offset
dim = 60
integrals = np.zeros((dim, 1), dtype=float)
w_x = np.zeros((dim, 1), dtype=float)
w_y = np.zeros((dim, 1), dtype=float)
wxx = np.zeros((dim, 1), dtype=float)
wyy = np.zeros((dim, 1), dtype=float)
aa = np.zeros((dim, 1), dtype=float)
bb = np.zeros((dim, 1), dtype=float)
cc = np.zeros((dim, 1), dtype=float)
fac = 0.05
for ii in range(0, dim):
w_x[ii] = (ii + 1) * width_x * fac
w_y[ii] = (ii + 1) * width_y * fac
print ('i=',ii, 'wx=',w_x[ii],' wy=',w_y[ii])
wxx[ii] = 2. * w_x[ii] * conv # convert half-axes to full axes
wyy[ii] = 2. * w_y[ii] * conv # convert half-axes to full axes
aa[ii] = (np.cos(theta) ** 2) / ( w_x[ii] ** 2) + (np.sin(theta) ** 2) / ( w_y[ii] ** 2)
bb[ii] = -(np.sin(2 * theta)) / (2 * w_x[ii] ** 2) + (np.sin(2 * theta)) / (2 * w_y[ii] ** 2)
cc[ii] = (np.sin(theta) ** 2) / ( w_x[ii] ** 2) + (np.cos(theta) ** 2) / ( w_y[ii] ** 2)
scidata_test = np.zeros(scidata_sl2.shape , dtype=float)
for ix in range(0, idy):
for iy in range(0, idx):
scidata_sl2[ix, iy] -= offset
scidata_sl2[ix, iy] /= sum
for ii in range(0, dim):
# if ((((ix - x) / w_x[ii] * 2) ** 2 + ((iy - y) / w_y[ii] * 2) ** 2) < 1):
if (( aa[ii]* (ix - x) ** 2 + cc[ii] * (iy - y) ** 2 + 2*bb[ii]* (ix - x)*(iy - y)) > 1):
integrals[ii] += scidata_sl2[ix, iy]
if (ii==25):
scidata_test[ix,iy] += 1.
# print (scidata_sl[0])
for ii in range(0, dim):
print ("wx=%04.3f" % w_x[ii], "wy=%04.3f" % w_y[ii], " int=%.04f" % integrals[ii])
print ('offset=', offset)
summ -= numsum * offset
integrals = np.zeros((dim, 1), dtype=float)
irradiance = np.zeros((dim, 1), dtype=float)
# w_x = np.zeros((dim, 1), dtype=float)
# w_y = np.zeros((dim, 1), dtype=float)
# wxx = np.zeros((dim, 1), dtype=float)
# wyy = np.zeros((dim, 1), dtype=float)
# aa = np.zeros((dim, 1), dtype=float)
# bb = np.zeros((dim, 1), dtype=float)
# cc = np.zeros((dim, 1), dtype=float)
sfac = int(securityfac / fac) - 1
ii = np.arange(dim)
w_x = (ii + 1) * width_x * fac
w_y = (ii + 1) * width_y * fac
wxx = 2. * w_x * conv # convert half-axes to full axes
wyy = 2. * w_y * conv # convert half-axes to full axes
# print(list(zip(ii, w_x, w_y, wxx, wyy)))
# for rotated ellipsee equations, see e.g.
# https://math.stackexchange.com/questions/426150/what-is-the-general-equation-of-the-ellipse-that-is-not-in-the-origin-and-rotate
aa = (np.cos(theta) ** 2) / ( w_x[ii] ** 2) + (np.sin(theta) ** 2) / ( w_y[ii] ** 2)
bb = -(np.sin(2 * theta)) / (2 * w_x[ii] ** 2) + (np.sin(2 * theta)) / (2 * w_y[ii] ** 2)
cc = (np.sin(theta) ** 2) / ( w_x[ii] ** 2) + (np.cos(theta) ** 2) / ( w_y[ii] ** 2)
scidata_test = np.zeros(scidata_sl2.shape , dtype=float)
scidata_test2 = np.zeros(scidata_sl2.shape , dtype=float)
for ix in range(idy):
for iy in range(idx):
scidata_log[ix, iy] = np.log(scidata_sl2[ix,iy])
scidata_sl2[ix, iy] -= offset
scidata_sl2[ix, iy] /= summ
for ii in range(dim):
# if ((((ix - x) / w_x[ii] * 2) ** 2 + ((iy - y) / w_y[ii] * 2) ** 2) < 1):
if (( aa[ii]* (ix - x) ** 2 + cc[ii] * (iy - y) ** 2 + 2*bb[ii]* (ix - x)*(iy - y)) > 1):
integrals[ii] += scidata_sl2[ix, iy]
if (ii==sfac):
scidata_test[ix,iy] += 2.
for ii in range(dim-1,0,-1):
if (integrals[ii] > beamdiameterdef):
for ix in range(idy):
for iy in range(idx):
if (( aa[ii]* (ix - x) ** 2 + cc[ii] * (iy - y) ** 2 + 2*bb[ii]* (ix - x)*(iy - y)) > 1):
scidata_test[ix,iy] += 1.
break
irradiance[0] = (1. - integrals[0]) / (wxx[0]*wyy[0])
for ii in range(1, dim-1):
irradiance[ii] = (integrals[ii-1] - integrals[ii]) / (wxx[ii]*wyy[ii]-wxx[ii-1]*wyy[ii-1])
for ii in range(0, dim):
irradiance[dim-ii-1] /= irradiance[0]
irradiance[dim-1] = 0;
for ii in range(0, dim):
print ("wx=%04.3f" % w_x[ii], "wy=%04.3f" % w_y[ii], " int=%.04f" % integrals[ii], " irr=%.04f" % irradiance[ii])
axarr[0, 1].imshow(scidata_log, origin='lower', cmap='gnuplot')
axarr[0, 1].set_title(titles[1], color=colors[1], fontsize=10)
axarr[1, 0].imshow(scidata_sl2, origin='lower', cmap='gnuplot')
axarr[1, 0].set_title(titles[2], color=colors[2], fontsize=10)
axarr[1, 1].imshow(scidata_test, origin='lower',cmap='gnuplot')
axarr[1, 1].set_title(titles[3], color=colors[3], fontsize=10)
f.tight_layout()
f.savefig(outfile+"_1.pdf")
# Interpolation object
fiy = interpolate.interp1d(integrals[:,0], wyy, assume_sorted = False)
fix = interpolate.interp1d(integrals[:,0], wxx, assume_sorted = False)
fry = interpolate.interp1d(irradiance[:,0],wyy, assume_sorted = False)
frx = interpolate.interp1d(irradiance[:,0],wxx, assume_sorted = False)
print ("ee=", beamdiameterdef, " interp=", fiy(beamdiameterdef))
diamy = float(fiy(beamdiameterdef))
diamx = float(fix(beamdiameterdef))
diary = float(fry(beamdiameterdef))
diarx = float(frx(beamdiameterdef))
print (diamx, diamy)
g,f2 = plt.subplots(1,1)
f2.plot(wyy,integrals,'ro',color='green', label=r'major axis: d={:3.1f}mm'.format(diamy, prec='.1'))
f2.plot(wxx,integrals,'ro',color='red', label=r'minor axis: d={:3.1f}mm'.format(diamx, prec='.1'))
f2.legend(loc="upper right", prop={'size':9})
plt.xlim([-10,upx])
f2.set_xlabel('diameter (mm)')
f2.set_ylabel('intensity coverage')
f2.grid()
f2.hlines(beamdiameterdef,0.,upx, linestyles='dashed')
g.tight_layout()
g.savefig(outfile+"_2.pdf")
idx2 = idx
idy2 = idy
if (idy > idx):
idx2 = idy
if (idx > idy):
idy2 = idx
fig = plt.figure()
ax1 = fig.add_subplot(121, projection='3d')
scidata = hdu.data
scidata_sl = scidata[rangey1:rangey2,rangex1:rangex2]
axarr[1].imshow(scidata_sl2, origin='lower', cmap='gnuplot')
axarr[1].set_title(titles[1], color=colors[1])
X = np.arange(idx2)
Y = np.arange(idy2)
XX, YY = np.meshgrid(X, Y)
ZZ = zeros((idx2, idy2), dtype=float)
for i in range(idx):
for j in range(idy):
ZZ[i, j] = scidata_sl[j][i]
for i in range(idx2):
for j in range(idy2):
if (ZZ[i, j] == 0):
ZZ[i, j] = ZZ[idx-1, idy-1]
axarr[2].imshow(scidata_test, origin='lower',cmap='gnuplot')
surf1 = ax1.plot_surface(XX, YY, ZZ, linewidth=0, antialiased=False, cmap=cm.coolwarm)
ax1.set_title(titles[0], color=colors[0])
ax2 = fig.add_subplot(122, projection='3d')
XL = arange(0, idx2, 1) ######## this is Y (50)
YL = arange(0, idy2, 1) ######## this is X (45)
XXL, YYL = np.meshgrid(XL, YL)
ZZL = zeros((idx2, idy2), dtype=float)
for i in range(idx):
for j in range(idy):
ZZL[i, j] = np.log(scidata_sl[j][i])
for i in range(idx2):
for j in range(idy2):
if (ZZL[i, j] == 0):
ZZL[i, j] = ZZL[idx-1, idy-1]
g,f2 = plt.subplots(1,1)
axes = plt.gca()
f2.plot(wyy,integrals,'ro',color=colors[1], label="major axis")
f2.plot(wxx,integrals,'ro',color=colors[0], label="minor axis")
f2.legend(loc="upper right", prop={'size':30})
axes.set_xlim([0,800*conv])
axes.set_xlabel('diameter (mm)')
axes.set_ylabel('coverage')
f2.grid()
surf2 = ax2.plot_surface(XXL, YYL, ZZL, linewidth=0, antialiased=False, cmap=cm.coolwarm)
ax2.set_title(titles[1], color=colors[1])
# draw vertical line from (70,100) to (70, 250)
f2.hlines(beamdiameterdef,0.,600*conv, linestyles='dashed')
# plt.show()
#################################################################
#################################################################
fig.savefig(outfile+"_3.pdf")
idx2 = idx
idy2 = idy
if (idy > idx):
idx2 = idy
if (idx > idy):
idy2 = idx
fig = plt.figure()
ax1 = fig.add_subplot(121, projection='3d')
if __name__ == '__main__':
if len(sys.argv) > 1: # there is at least one argument
print('Number of arguments:', len(sys.argv), 'arguments.')
print('Argument List:', str(sys.argv))
scidata = hdu.data
scidata_sl = scidata[rangey1:rangey2,rangex1:rangex2]
# print (scidata_sl[0])
file = sys.argv[1]
rangey1 = int(float(sys.argv[2]))
rangey2 = int(float(sys.argv[3]))
rangex1 = int(float(sys.argv[4]))
rangex2 = int(float(sys.argv[5]))
conv = float(sys.argv[6])
X = arange(0, idx2, 1) ######## this is Y (50)
Y = arange(0, idy2, 1) ######## this is X (45)
XX, YY = np.meshgrid(X, Y)
ZZ = zeros((idx2, idy2), dtype=float)
for i in range(idx):
for j in range(idy):
ZZ[i, j] = scidata_sl[j][i]
run_script(file,rangey1,rangey2,rangex1,rangex2)
surf1 = ax1.plot_surface(XX, YY, ZZ, linewidth=0, antialiased=False, cmap=cm.coolwarm)
ax1.set_title(titles[0], color=colors[0])
plt.show()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment