SIM Pattern Leakage Test

The purpose of this notebook is to test how much leakage occurs through the other mask holes relative to the desired peaks for different situations.

One thing to note is that the SLM effects the electric field but the sample responds to the intensity this means that during this analysis the result of the FFTs must be squared.

Remembering our previous results estimating the modulation contrast due to small variations in the two SIM beams, it doesn't take a lot to create an interference pattern. If we want to keep the spurious intereference patterns to less than 5% of the main one than we need the main peaks to have about 1,000 times more energy than the spurious ones.


In [1]:
#%%px --local
%pylab inline
#nice plotting
#import seaborn as sns
import numexpr as ne
#for minimizing the difference between the desired frequency and the calculated one
from scipy.optimize import minimize
#Need to be able to do polynomial fitting
from sklearn.preprocessing import PolynomialFeatures
from sklearn import linear_model
#PIL allows us to write binary images, though we need a cludge, see 'Writing Binary Files.ipynb'
from PIL import Image
import os
import zipfile

from matplotlib.patches import Circle, Wedge, Polygon

from scipy.ndimage import gaussian_filter

import seaborn as sns

from skimage.draw import circle, polygon

set_cmap('gnuplot2')

from matplotlib.colors import LogNorm


Populating the interactive namespace from numpy and matplotlib
<matplotlib.figure.Figure at 0x1153550b8>

In [2]:
#newest one
def pattern_gen(angle, period, onfrac = 0.5, phase_idx = 0., phase_offset = 0., nphases =5, sizex =2048, sizey =1536, SIM_2D = True, eightbit=False):
    '''
    Generates a binary SLM pattern for SIM
    
    Generates a sine wave and then makes it binary with a given duty cycle
    
    Parameters
    ----------
    angle : array_like
        Defines the pattern orientation
    period : float
        Defines the period of the pattern
    onfrac : float
        The fraction of on pixels in the pattern
    phase_idx : int
        The phase of the pattern (see `nphases` for more info)
    phase_offset : float
        The offset in phase, mostly used for aligning patterns of different colors
    nphases : int
        the number of phases
    sizex : int
        size of the pattern
    sizey : int
        size of the pattern

    Returns
    -------
    pattern : ndarray
        A binary array representing a single bitplane to send to the SLM
    '''
    
    if not 0 < onfrac < 1:
        raise ValueError('onfrac must have a value between 0 and 1. onfrac = {}'.format(onfrac))
    
    
    if SIM_2D:
        #Then we only want to take steps of 2pi/n in illumination which means pi/n at the SLM
        phase_step = 2
    else:
        phase_step = 1
    
    yy, xx = meshgrid(arange(sizey),arange(sizex),indexing='ij')
    
    #here's the pattern frequency
    freq = 2*pi/period
    
    phi = (phase_idx/nphases/phase_step)*(2*pi)+phase_offset
    
    #our sine goes from -1 to 1 while on frac goes from 0,1 so move onfrac into the right range
    onfrac = onfrac*2-1
    
    #given angle is the angle of the fourier spots from the 'i' vector, so pattern should be 90 degrees from this
    a=angle
    
    #do the evaluation
    toreturn = np.sin(freq*(np.cos(a)*xx+np.sin(a)*yy)+phi)
    if eightbit:
        toreturn = toreturn - toreturn.min()
        toreturn /= toreturn.max()
        toreturn = (toreturn * 255).astype("uint8")
    else:
        toreturn = toreturn > onfrac
    
    return toreturn

In [3]:
my_angle = deg2rad(4.76)
size=128
pat =  pattern_gen(my_angle, 8, 0.5, 0, 0, 5,size,size)
pat_fft = fftshift(fftn(ifftshift(pat)))
l = size/8
matshow(pat,cmap='Greys')
arrow(size/2,size/2,cos(my_angle)*l, sin(my_angle)*l,width=2/2048*size,ec='c')
grid('off')
matshow(abs(pat_fft),cmap='Greys',norm=LogNorm())
grid('off')
matshow(gaussian_filter(abs(pat_fft),10/2048*size),norm=LogNorm(),cmap='Greys')
arrow(size/2,size/2,cos(my_angle)*l, sin(my_angle)*l,width=2/2048*size,ec='c',length_includes_head=True)
grid('off')



In [6]:
my_angle = deg2rad(60)
size=128
pat =  pattern_gen(my_angle, 8, 0.5, 0, 0, 5,size,size, eightbit=True)
pat_fft = fftshift(fftn(ifftshift(pat)))
l = size/8
matshow(pat,cmap='Greys')
arrow(size/2,size/2,cos(my_angle)*l, sin(my_angle)*l,width=2/2048*size,ec='c')
grid('off')
matshow(abs(pat_fft),cmap='Greys',norm=LogNorm())
grid('off')
matshow(gaussian_filter(abs(pat_fft),10/2048*size),norm=LogNorm(),cmap='Greys')
arrow(size/2,size/2,cos(my_angle)*l, sin(my_angle)*l,width=2/2048*size,ec='c',length_includes_head=True)
grid('off')



In [4]:
#%%px --local
def ideal_period(wavelength, NA = 0.85):
    '''
    All units are in mm
    '''
    pixel_size = 8.2/1000 #pixel size in mm for QXGA display (4DD)
    fl = 250 #focal length of lens in mm
    fl2 = 300 #focal length of the second lens
    ftube = 200 #focal length of the tube lens, for Nikon this is 200 mm
    wl = wavelength/10**6 #wavelength of light
    mag = 1/100
    sigma = sqrt(2) * 12/pixel_size/4 #std dev of gaussian beam in units of pixels at the SLM
    pupil_diameter = 2*NA*mag*fl2    #Size of pupil image at first fourier plane
    hole_radius = 2*wl*fl/(2* pi * sigma *sqrt(2) * pixel_size) #this is the limit of hole size
    hole_radius = 0.1/2# this is more reasonable (50 um)
    period = wl * fl * (1/(pupil_diameter/2 - hole_radius))/ pixel_size #in mm
    
    return period

In [5]:
fig, axs = subplots(3,2, figsize=(12,18))
for i, (ax1, ax2) in enumerate(axs):
    my_angle= deg2rad(4.76-i*60)
    pat =  pattern_gen(my_angle, 7, 0.5, 0, 0, 5,2048,2048)
    pat_fft = fftshift(fftn(ifftshift(pat)))
    ax1.matshow(pat)
    ax2.matshow(abs(pat_fft),norm=LogNorm())



In [6]:
def peak_amp(angle, period, onfrac = 0.5, phaseInd = 0., phaseOffset = 0., nphases = 5, sizex =2048, sizey =1536):
    '''
    Find the precise pattern period
    
    Using 2nd order polynomial fit along either axis
    '''
    
    #this is a standin for nextpow2, I don't know if there's a more efficient one, but its not going to
    #be a bottle neck
    n = 1<<(min(sizex,sizey)-1).bit_length()
    
    my_pat = pattern_gen(angle, period, onfrac, phaseInd, phaseOffset, nphases, n, n)
    
    my_pat_fft = fftshift(fftn(ifftshift(my_pat)))
    
    my_pat_fft_abs = abs(my_pat_fft)
    
    dc_loc = ((np.array(my_pat_fft.shape)+1)//2)
    
    #make mask
    mask = ones_like(my_pat_fft_abs)
    mask[circle(*dc_loc,radius=10)]=0
    mask[:dc_loc[0]]=0
    
    #mask data and find next biggest peak
    masked_fft_data = mask*my_pat_fft_abs
    max_loc = unravel_index(masked_fft_data.argmax(),masked_fft_data.shape)
    #pull the 3x3 region around the peak
    region_size = 3
    start = -(region_size-1)//2
    end = region_size+start

    my_pat_fft_suby = masked_fft_data[max_loc[0]+start:max_loc[0]+end,max_loc[1]]
    my_pat_fft_subx = masked_fft_data[max_loc[0],max_loc[1]+start:max_loc[1]+end]
    
    x = arange(start,end)
    
    try:
        xfit = polyfit(x, my_pat_fft_subx,2)
        yfit = polyfit(x, my_pat_fft_suby,2)
    except TypeError as e:
        print('start = {}, end = {}'.format(start,end))
        print(max_loc)
        print(my_pat_fft_subx)
        amp=nan
        x0 = nan
        y0 = nan
    else:

        xpoly = poly1d(xfit)

        x0 = -xfit[1]/(2*xfit[0])
        y0 = -yfit[1]/(2*yfit[0])
        
        amp = xpoly(0)/xpoly(x0)*yfit[2]

    precisepeak = array(max_loc)+array([x0,y0])-array(dc_loc)
    
    preciseangle = arctan2(precisepeak[0],precisepeak[1])
    
    precise_period = n/norm(precisepeak)
    
    return dict(precise_period=precise_period,pat = my_pat, pat_fft = my_pat_fft,amp = amp)

In [7]:
fig, axs = subplots(3,2, figsize=(12,18))
for i, (ax1, ax2) in enumerate(axs):
    angle =deg2rad(4.76-i*60)
    pat_dict =  peak_amp(array(angle), 7, 0.5, 0, 0, 5,2048,2048)
    ax1.matshow(pat_dict['pat'])
    ax2.matshow(abs(pat_dict['pat_fft']),norm=LogNorm())
    print(pat_dict['amp'])


810951.300313
1001514.03236
964339.981437