lemon-sandbox/260721-atomic-column-thickness/260721.ipynb

In [1]:

import os

from tqdm.notebook import tqdm

import numpy as np
from matplotlib import pyplot as plt
from scipy.io import loadmat
from scipy.ndimage import gaussian_filter
from scipy.optimize import curve_fit
from scipy.special import erf

import kemstem

In [2]:

def calculate_R(raw_data, regressed):
    sse = np.sum(np.square(raw_data - regressed))
    sst = np.sum(np.square(raw_data - raw_data.mean()))
    return 1 - sse/sst

def fit(f, xdata, ydata, plot=False, **kwargs):
    xdata = np.array(xdata).flatten()
    ydata = np.array(ydata).flatten()
    if plot:
        plt.plot(xdata, ydata, 'k.')
        if 'sigma' in kwargs:
            plt.errorbar(xdata, ydata, kwargs['sigma'], fmt='none', c='k')
    popt, pcov = curve_fit(f, xdata, ydata, **kwargs)
    perr = np.sqrt(np.diag(pcov))
    R = calculate_R(ydata, f(xdata, *popt))
    if plot:
        xlim, ylim = plt.xlim(), plt.ylim()
        plt.xlim(*xlim)
        plt.ylim(*ylim)
        x = np.linspace(*xlim, 100)
        plt.plot(x, f(x, *popt), 'r--')
        plt.show()
    return popt, perr, R

In [22]:

results = loadmat('./Si2V1-Niter1000.mat')
# results = loadmat('./Si5V2-Niter1000.mat')

In [23]:

roi = results['p']['object_ROI'][0][0][0]
roix = roi[0].T
roiy = roi[1]

DZ = results['outputs']['z_distance'].item().item()

recon_object = results['object'][roix, roiy]

In [24]:

LAYER = 10

plt.imshow(np.angle(recon_object[:, :, LAYER]), cmap='gray')
plt.axis('off')

Out [24]:

(np.float64(-0.5), np.float64(572.5), np.float64(572.5), np.float64(-0.5))

In [25]:

image = np.angle(recon_object[:,:,LAYER])

blurred_image = gaussian_filter(image, 1)
distance = 2
threshold = 0.5

c0 = kemstem.find_columns(blurred_image, distance=distance, threshold=threshold)

plt.imshow(image, cmap='gray')
plt.axis('off')
plt.plot(c0[:, 1], c0[:, 0], 'r.', markersize=2)
plt.show()
print(f'{c0.shape[0]} columns found')

1464 columns found

In [26]:

# test_it = 100
# window_dim = 5
# cf,errs,opts,data_fits = kemstem.refine_columns(blurred_image,c0[test_it,:],window_dim)

# visualization_window = 5*window_dim
# fig,ax = plt.subplots(2,2,constrained_layout=True)
# ax[0,0].matshow(data_fits[:,:,0,0],cmap='gray')
# ax[0,1].matshow(data_fits[:,:,0,1],cmap='gray')
# ax[1,0].matshow(blurred_image,cmap='gray')
# ax[1,0].plot(c0[test_it,1],c0[test_it,0],'b.') # original unfit position shown in blue
# ax[1,0].plot(cf[0,1],cf[0,0],'r.') # fit position shown in red
# ax[1,0].set_xlim([c0[test_it,1]-visualization_window,c0[test_it,1]+visualization_window])
# ax[1,0].set_ylim([c0[test_it,0]+visualization_window,c0[test_it,0]-visualization_window])
# ax[1,1].matshow(data_fits[:,:,0,0]-data_fits[:,:,0,1],cmap='gray',vmin=-.1,vmax=.1) # fit residual
# _ = [tax.axis('off') for tax in ax.ravel()]

In [27]:

# cf, errs, opts, data_fits = kemstem.refine_columns(image, c0, window_dim)

In [28]:

# plt.imshow(image, cmap='gray')
# plt.axis('off')
# plt.plot(cf[:, 1], cf[:, 0], 'r.', markersize=2)
# plt.show()
# print(f'Refined columns')

In [29]:

WINDOW = 1
cx, cy = c0[100]
wx_start, wx_end = cx-WINDOW, cx+WINDOW+1
wy_start, wy_end = cy-WINDOW, cy+WINDOW+1
plt.imshow(np.angle(recon_object)[wx_start:wx_end, wy_start:wy_end, LAYER], cmap='gray')
plt.axis('off')

Out [29]:

(np.float64(-0.5), np.float64(2.5), np.float64(2.5), np.float64(-0.5))

In [30]:

columns = []

for i, (cx, cy) in enumerate(tqdm(c0)):
    wx_start, wx_end = cx-WINDOW, cx+WINDOW+1
    wy_start, wy_end = cy-WINDOW, cy+WINDOW+1
    column = np.angle(recon_object)[wx_start:wx_end, wy_start:wy_end, 1:].mean(axis=(0, 1))
    columns.append(column)

  0%|          | 0/1464 [00:00<?, ?it/s]

In [92]:

def double_erf(x, p0, p1, q, x0, h):
    return p0 + p1 * (erf(q*(x - x0)) - erf(q*(x - x0 - h)))

In [114]:

bounds = [[-0.2, 0, 0, -5, 0], [0.3, 1, 10, 15, 25]]

for i in [700]:
    popt, _, _ = fit(
        double_erf, np.arange(len(columns[0])), columns[i], plot=True,
        p0=[0.1, 0.05, 1, 1, 15],
        bounds=bounds)
    print(popt)

[ 0.09843992  0.17781052  0.36135616  3.09665732 12.8743372 ]

In [115]:

heights = []
heights_error = []

for column in tqdm(columns):
    popt, perr, _ = fit(
        double_erf, np.arange(len(column)), column,
        plot=False, p0=popt, bounds=bounds, maxfev=10000)
    h, dh = popt[-1], perr[-1]
    heights.append(h * DZ)
    heights_error.append(dh * DZ)

heights = np.array(heights)
heights_error = np.array(heights_error)

  0%|          | 0/1464 [00:00<?, ?it/s]

In [116]:

plt.plot(np.sort(heights), 'k')

Out [116]:

[<matplotlib.lines.Line2D at 0x7f2226feea20>]

In [117]:

plt.imshow(np.angle(recon_object).sum(axis=-1), cmap='gray')
plt.axis('off')
plt.scatter(c0[:, 1], c0[:, 0], c=np.array(heights), s=8, cmap='magma_r', vmax=200, vmin=100)

cbar = plt.colorbar()

plt.show()

In [ ]: