<jupyter_start><jupyter_text>`beam_size()` Algorithm**<NAME>****June 2020, Version 5**The ISO 11146 standard recommends masking the image with a rotated rectangular mask around the center of the image.  This notebook explains how that was implemented and then shows some results for artifically generated images of non-circular Gaussian beams. As noise increases, the first-order parameters (beam center) are robust, but the second-order parameters (diameters) are shown to be much more sensitive to image noise.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import laserbeamsize as lbs

def side_by_side_plot(h,v,xc,yc,dx,dy,phi,noise=0,offset=0):
    test = lbs.beam_test_image(h,v,xc,yc,dx,dy,phi,noise=noise)
    xc_found, yc_found, dx_found, dy_found, phi_found = lbs.beam_size(test)

    plt.subplots(1,2,figsize=(12,5))
    plt.subplot(1,2,1)
    plt.imshow(test,cmap='gray')
    plt.plot(xc,yc,'ob',markersize=2)
    plt.title('Original (%d,%d), dx=%.0f, dy=%.0f, $\phi$=%.0f°' % (xc,yc,dx,dy,np.degrees(phi)))

    plt.subplot(1,2,2)
    plt.imshow(test,cmap='gray')
    xp,yp = lbs.rotated_rect_arrays(xc_found, yc_found, dx_found, dy_found, phi_found)
    plt.plot(xp,yp,':y')
    plt.plot(xc_found,yc_found,'ob',markersize=2)
    plt.title('Found (%d,%d), dx=%.0f, dy=%.0f, $\phi$=%.0f°' % (xc_found,yc_found,dx_found,dy_found,np.degrees(phi_found)))<jupyter_output><empty_output><jupyter_text>Integration AreaISO 11146-3 states::> All integrations ... are performed on a rectangular integration area which is centred to the beam centroid, defined by the spatial first order moments, orientated parallel to the principal axes of the power density distribution, and sized three times the beam widths $d_{\sigma x}$ and $d_{\sigma y}$...This turned out to be surprisingly fiddly (most likely because masked numpy arrays did't work as I expected).  In the pictures that follow, the dotted rectangle shows the final integration area.  Algorithm OverviewThis shows the effect of background removal, then rectangle masking, then locating the first<jupyter_code>xc=300
yc=200
dx=100
dy=50
phi=np.radians(30)
h=600
v=400
max_value=1023
noise = 25

# generate test image
test = lbs.beam_test_image(h,v,xc,yc,dx,dy,phi,noise=noise)

plt.subplots(3,2,figsize=(12,12))
plt.subplot(3,2,1)
plt.imshow(test)
plt.colorbar()
plt.xticks([])
plt.yticks([])
plt.title('Original')

# corners
corners = lbs.corner_mask(test)
plt.subplot(3,2,2)
plt.imshow(corners)
plt.colorbar()
plt.xticks([])
plt.yticks([])
plt.title('Corners Used')

# remove background
back, sigma = lbs.corner_background(test)
offset = int(back+3*sigma)
zeroed = lbs.subtract_threshold(test, offset)

plt.subplot(3,2,3)
plt.imshow(zeroed)
plt.colorbar()
plt.xticks([])
plt.yticks([])
plt.title('Background Removed')

# first guess at beam parameters
xc, yc, dx, dy, phi = lbs.basic_beam_size(zeroed)

plt.subplot(3,2,4)
plt.imshow(zeroed)
plt.colorbar()
plt.xticks([])
plt.yticks([])
xp,yp = lbs.ellipse_arrays(xc, yc, dx, dy, phi)
plt.plot(xp,yp,':y')
plt.title('Ellipse estimate using entire image')

mask = lbs.rotated_rect_mask(zeroed, xc, yc, dx, dy, phi)

plt.subplot(3,2,5)
plt.imshow(mask)
plt.colorbar()
plt.xticks([])
plt.yticks([])
xp,yp = lbs.ellipse_arrays(xc, yc, dx, dy, phi)
plt.plot(xp,yp,':b')
plt.title('Integration Region')

masked_image = np.copy(zeroed)
masked_image[mask < 1] = 0       # zero all values outside mask

plt.subplot(3,2,6)
plt.imshow(masked_image)
plt.colorbar()
plt.xticks([])
plt.yticks([])
xp,yp = lbs.rotated_rect_arrays(xc, yc, dx, dy, phi)
plt.plot(xp,yp,':y')
xp,yp = lbs.ellipse_arrays(xc, yc, dx, dy, phi)
plt.plot(xp,yp,':y')
plt.title('First result using only integration region')

plt.show()<jupyter_output><empty_output><jupyter_text>Integration Area Tests Test 1 Centered, Horizontal, away from edges<jupyter_code>xc=320
yc=200
dx=100
dy=40
phi=np.radians(0)
h=600
v=400

side_by_side_plot(h,v,xc,yc,dx,dy,phi)
plt.show()<jupyter_output><empty_output><jupyter_text>Test 2 Corner, Horizontal<jupyter_code>xc=550
yc=50

side_by_side_plot(h,v,xc,yc,dx,dy,phi)
plt.show()<jupyter_output><empty_output><jupyter_text>Test 3 Center, tilted 30°<jupyter_code>xc=320
yc=200
dx=100
dy=10
phi=np.radians(30)

side_by_side_plot(h,v,xc,yc,dx,dy,phi)
plt.show()<jupyter_output><empty_output><jupyter_text>Test 4 Corner, tilted 30°<jupyter_code>xc=120
yc=350

side_by_side_plot(h,v,xc,yc,dx,dy,phi)
plt.show()<jupyter_output><empty_output><jupyter_text>Test 5 Center, tilted 70°<jupyter_code>xc=320
yc=200
dx=100
dy=10
phi=np.radians(70)

side_by_side_plot(h,v,xc,yc,dx,dy,phi)
plt.show()<jupyter_output><empty_output><jupyter_text>Tests with Poisson noise Test 1. Simple, noise-free rotated elliptical beamIn this and all rest of the test functions, the maximum value in the test array is 256.No gaussian noise, works fine!<jupyter_code>xc=400
yc=150
dx=150
dy=40
phi=np.radians(20)
h=640
v=480

side_by_side_plot(h,v,xc,yc,dx,dy,phi,noise=0)
plt.show()<jupyter_output><empty_output><jupyter_text>Test 2: Rotated elliptical beam with 4% Poisson noiseHere the added noise is enough to significantly affect the estimate of beam widths<jupyter_code>xc=200
yc=200
dx=120
dy=30
phi=np.radians(20)
h=400
v=400
max_value = 255
noise = 0.04 * max_value

test = lbs.beam_test_image(h, v, xc, yc, dx, dy, phi, noise=noise, max_value=max_value)
title = 'Original (%d,%d), dx=%d, dy=%d' % (xc,yc,dx,dy)

lbs.visual_report(test, title)
#plt.savefig('plot2.pdf', bbox_inches = 'tight')
plt.show()
import imageio
beam = imageio.imread("t-hene.pgm")
lbs.visual_report(beam)
plt.savefig('hene.png', bbox_inches = 'tight')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14681978798586573, 'fraction_numerical': 0.03657243816254417, 'mean_word_length': 6.702040816326531, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '205888', 'license': 'MIT', 'n_tokens_mistral': 2480, 'n_tokens_neox': 2213, 'n_words': 494, 'path': 'docs/02-Test-Beam.ipynb', 'repo_name': 'SReich-EMI/laserbeamsize', 'size': '5665'}