small analysis program for images taken with the LASER spot and a reference grid

beam_diameters_from_Laser_image.py

+#!/usr/bin/python
+
+"""
+A small analysis program for images taken with the LASER spot and a reference grid. 
+The image needs to have been transformed to FITS format (preferably the green color part)
+
+The program read the fits-file, and fits an elliptic Gaussian to the LASER spot
+It then subtracts the background (calculated from beyond an ellipse of "securityfac" times the 
+fitted half lengths and widths and normalizes the content. 
+
+The background-subtracted image is then integrated in concentric ellipses around the fitted center, 
+and the image integral contained outside the ellipses calculated. 
+
+The result is plot of half width/length vs. image integral, from which one can read off 
+the beam diameters at the reference value of 1/e^2 (indicated by a horizontal line)
+
+CONFIGURATION:
+1) Open the image first with Gimp and select a suitable cutout (containing the full Laser image and some constant background)
+2) Write the full image name in the parameter "file"
+3) Use the cursor to pick the edge points of the cutout
+4) Write these edge coordinates into the variables "rangex1, rangex2, rangey1 and rangey2"
+5) Use the reference grid to calculate the conversion from points to mm (always using Gimp)
+6) Write the conversion into the parameter "conv"
+7) Run the program
+"""
+from scipy import optimize as opt
+import astropy
+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
+
+#
+# user configuration: 
+#
+# image file in FITS format (green):
+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
+# conversion from points to mm (please use GIMP and the reference grid on the picture!
+conv   = 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)
+    width_y = float(width_y)
+    a = (np.cos(theta) ** 2) / (    width_x ** 2) + (np.sin(theta) ** 2) / (    width_y ** 2)  
+    b = -(np.sin(2 * theta)) / (2 * width_x ** 2) + (np.sin(2 * theta))  / (2 * width_y ** 2)   
+    c = (np.sin(theta) ** 2) / (    width_x ** 2) + (np.cos(theta) ** 2) / (    width_y ** 2)  
+
+#    print (height, center_x, center_y, width_x, width_y, theta, offset)
+    if (height<0):
+        return lambda x, y: 99999999.
+    if (center_x < 0.):
+        return lambda x, y: 99999999.
+    if (center_y < 0.):
+        return lambda x, y: 99999999.
+    if (width_x < 0.):
+        return lambda x, y: 99999999.
+    if (width_y < 0.):
+        return lambda x, y: 99999999.
+    if (offset < 0.):
+        return lambda x, y: 99999999.
+    if (theta > 3.*np.pi):
+        return lambda x, y: 99999999.
+    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)
+    the gaussian parameters of a 2D distribution by calculating its
+    moments """
+    total = data.sum()
+    X, Y = np.indices(data.shape)
+    x = (X * data).sum() / total
+    y = (Y * data).sum() / total
+    col = data[:, int(y)]
+    width_x = np.sqrt(np.abs((np.arange(col.size) - y) ** 2 * col).sum() / col.sum())
+    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)
+
+#    print height, x, y, width_x, width_y, -45*np.pi/180, offset
+    return height, x, y, 2.*width_x, width_y, -45.*np.pi/180, offset
+
+
+def fitgaussian(data):
+    """Returns (height, x, y, width_x, width_y)
+    the gaussian parameters of a 2D distribution found by a fit"""
+    params = moments(data)
+    errorfunction = lambda p: np.ravel(gaussian(*p)(*np.indices(data.shape)) -
+                                       data)
+    p, success = opt.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):  
+    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])
+
+    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)
+    (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): 
+    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])
+
+
+axarr[1].imshow(scidata_sl2, origin='lower', cmap='gnuplot')
+axarr[1].set_title(titles[1], color=colors[1])
+
+axarr[2].imshow(scidata_test, origin='lower',cmap='gnuplot')
+
+
+
+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()
+
+# draw vertical line from (70,100) to (70, 250)
+f2.hlines(beamdiameterdef,0.,600*conv, linestyles='dashed')
+
+#################################################################
+#################################################################
+
+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]
+#    print (scidata_sl[0])
+
+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]
+
+surf1 = ax1.plot_surface(XX, YY, ZZ, linewidth=0, antialiased=False, cmap=cm.coolwarm)
+ax1.set_title(titles[0], color=colors[0])
+
+plt.show()