lemon-sandbox/260721-atomic-column-thickness/260727-graphing.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 [21]:

plt.rcParams['font.size'] = 14

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 [3]:

results = {
    '2V': loadmat('./Si2V1-Niter1000.mat'),
    '5V': loadmat('./Si5V1-Niter1000.mat'),
    '8V': loadmat('./Si8V2-Niter1000.mat'),
    '30V': loadmat('./Si30V3-Niter1000.mat'),
}

In [4]:

DZ = {}
recons = {}

for key in results.keys():
    result = results[key]
    roi = result['p']['object_ROI'][0][0][0]
    roix = roi[0].T
    roiy = roi[1]

    DZ.update({key: result['outputs']['z_distance'].item().item()})
    recons.update({key: np.unwrap(np.angle(result['object']))[roix, roiy]})

In [5]:

column_positions = {}

for key in results.keys():
    recon = recons[key]
    image = recon[:,:,recon.shape[-1] // 2]
    blurred_image = gaussian_filter(image, 1)
    distance = 2
    threshold = 0.5
    c0 = kemstem.find_columns(blurred_image, distance=distance, threshold=threshold)

    column_positions.update({key: c0})

In [6]:

WINDOW = 1

column_lists = {}

for key in results.keys():

    c0 = column_positions[key]
    recon = recons[key]

    column_list = []
    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 = recon[wx_start:wx_end, wy_start:wy_end, :].mean(axis=(0, 1))
        column_list.append(column)

    column_lists.update({key: column_list})

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

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

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

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

In [7]:

def double_erf(x, p, q, x0, x1):
    return p * (erf(q*(x - x0)) - erf(q*(x - x1)))

def double_erf_plus_const(x, p0, p1, q, x0, x1):
    return p0 + double_erf(x, p1, q, x0, x1)

In [8]:

crystalline_heights = {}
crystalline_height_errs = {}

for key in results.keys():

    columns = column_lists[key]
    dz = DZ[key]

    crystalline = []
    crystalline_err = []

    for column in tqdm(columns):
        popt, perr, _ = fit(
            double_erf, np.arange(len(column)), column, plot=False,
            p0=[0.1,  0.3,  3, 12],
            # bounds=bounds,
            maxfev=10000
            )
        p1, q, x0, x1 = popt
        dp1, dq, dx0, dx1 = perr

        crystalline.append((x1 - x0) * dz)
        crystalline_err.append(dz * (dx0 + dx1) / 2)

    crystalline_heights.update({key: np.array(crystalline)})
    crystalline_height_errs.update({key: np.array(crystalline_err)})

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

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

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

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

In [9]:

safe_height_indices = {}

for key in results.keys():
    crystalline = crystalline_heights[key]
    crystalline_err = crystalline_height_errs[key]

    # safe_height_index = np.where(crystalline_err < crystalline / 10)
    safe_height_index = np.where(crystalline)
    safe_height_indices.update({key: safe_height_index})

In [10]:

for key, recon in recons.items():
    print(key)
    print(recon.shape[:-1])
    print([(recon.shape[0] * i) // 8 for i in [1, 3, 5, 7]])

2V
(573, 573)
[71, 214, 358, 501]
5V
(580, 580)
[72, 217, 362, 507]
8V
(581, 581)
[72, 217, 363, 508]
30V
(541, 541)
[67, 202, 338, 473]

In [11]:

thickness_rois = {
    '2V' : [[214, 358], [214, 358]],
    '5V' : [[217, 362], [ 72, 217]],
    '8V' : [[217, 363], [ 72, 217]],
    '30V': [[338, 473], [ 202, 338]]
}

In [12]:

roi_indices = {}

for key in results.keys():
    c0 = column_positions[key]
    indices = []
    (x0, x1), (y0, y1) = thickness_rois[key]
    idx = np.where(
        (c0[:,0] >= y0) & (c0[:,0] <= y1) & (c0[:,1] >= x0) & (c0[:,1] <= x1) 
    )
    indices.append(idx)
    roi_indices.update({key: indices})

In [13]:

roi_crystalline_height_means = {}
roi_crystalline_height_stds = {}
roi_amorphous_height_means = {}
roi_amorphous_height_stds = {}

for key in results.keys():

    safe_idx = safe_height_indices[key]
    roi_idx = roi_indices[key]
    total_thickness = recons[key].shape[-1] * DZ[key]

    idx = np.intersect1d(safe_idx, roi_idx)

    height = crystalline_heights[key]

    crystalline_mean = height[idx].mean()
    crystalline_std = height[idx].std()
    amorphous_mean = (total_thickness - crystalline_mean) / 2
    amorphous_std = crystalline_std / 2

    roi_crystalline_height_means.update({key: crystalline_mean})
    roi_crystalline_height_stds.update({key: crystalline_std})
    roi_amorphous_height_means.update({key: amorphous_mean})
    roi_amorphous_height_stds.update({key: amorphous_std})

In [ ]:

key = '8V'
roi = thickness_rois[key]
(xi, xf), (yi, yf) = roi
recon = recons[key][yi:yf, xi:xf]


img_roi = recon[:,:,recon.shape[-1]//2]

cols = column_positions[key][roi_indices[key][0][0]]
cx, cy = cols.T
cx -= yi
cy -= xi


plt.imshow(img_roi)
for i in range(len(cx)):
    plt.text(cy[i], cx[i], str(i))

# plt.plot(*cols[43].T[::-1], 'r.')

In [16]:

#  block for analysing divergent regressions


idx = np.where(crystalline_height_errs[key][roi_indices[key]] > crystalline_heights[key][roi_indices[key]] / 10)[0]

for i in idx[::]:
    c = column_lists[key][i]
    z = np.arange(len(c)) * DZ[key]

    popt, perr, _ = fit(
        double_erf, np.arange(len(c)), c, plot=True,
        p0=[0.1,  0.3,  3, 12],
        # bounds=bounds,
        maxfev=10000
        )
    
    p1, q, x0, x1 = popt
    dp1, dq, dx0, dx1 = perr

    print(f"Column index: {i}")
    print(f"Thickness: {(x1 - x0) * dz:.2f}, Error: {dz * (dx0 + dx1) / 2:.2f}, Ratio: {(dz * (dx0 + dx1) / 2) /  (x1 - x0) * dz:.2f}")

    plt.show()

Column index: 0
Thickness: 244.47, Error: 8.91, Ratio: 17.55

In [18]:

for key in results.keys():

    recon = recons[key]
    heights = crystalline_heights[key]
    c0 = column_positions[key]

    crystalline_mean = roi_crystalline_height_means[key]
    crystalline_std = roi_crystalline_height_stds[key]
    amorphous_mean = roi_amorphous_height_means[key]
    amorphous_std = roi_amorphous_height_stds[key]


    safe_idx = safe_height_indices[key]
    roi_idx = roi_indices[key]
    idx = np.intersect1d(safe_idx, roi_idx)

    print(key)
    print(f"Crystalline thickness: {crystalline_mean:.3f} ± {crystalline_std:.3f}")
    print(f"Amorphous thickness:    {amorphous_mean:.3f} ±  {amorphous_std:.3f}")
    print(f"(# of counted columns = {len(idx)}/{len(roi_idx[0][0])})")

    (x0, x1), (y0, y1) = thickness_rois[key]
    plt.imshow((recon).sum(axis=-1), cmap='gray')
    plt.plot([x0, x0, x1, x1, x0], [y0, y1, y1, y0, y0], lw=1, color='white', linestyle='--')

    plt.axis('off')
    plt.scatter(
        c0[:, 1][idx], c0[:, 0][idx], c=np.array(heights)[idx], s=8, cmap='magma_r',
        vmin=heights[safe_idx].min(),
        vmax=heights[safe_idx].max()
        )
    cbar = plt.colorbar()
    cbar.ax.set_ylabel('Atomic column thickness [A]')
    plt.show()

2V
Crystalline thickness: 151.257 ± 64.435
Amorphous thickness:    46.872 ±  32.217
(# of counted columns = 102/102)

5V
Crystalline thickness: 175.627 ± 8.677
Amorphous thickness:    52.687 ±  4.339
(# of counted columns = 102/102)

8V
Crystalline thickness: 183.862 ± 19.102
Amorphous thickness:    76.069 ±  9.551
(# of counted columns = 118/118)

30V
Crystalline thickness: 75.191 ± 77.787
Amorphous thickness:    302.404 ±  38.894
(# of counted columns = 101/101)

In [19]:

voltages = []
damage = []
damage_err = []

for key in results.keys():
    voltages.append(float(key[:-1]))
    damage.append(roi_amorphous_height_means[key])
    damage_err.append(roi_amorphous_height_stds[key])

plt.figure(figsize=(4, 4), dpi=300)

plt.errorbar(voltages, damage, damage_err,
    fmt="o", capsize=3, capthick=1, elinewidth=1, label="Current study",
    zorder=100, color='k'
    )
plt.plot([30, 8, 5, 2], [220, 70, 40, 10], '*', label='Uzuhashi et al. (2024)\n(Amorphous)')
plt.plot([30, 8, 5, 2], [290, 100, 70, 70], 'p', label='Uzuhashi et al. (2024)\n(Total damage)')
plt.plot([30, 5, 2], [220, 66, 31], 's', label='Burnett et al. (2015)')
plt.plot([30, 5, 2], [220, 25, 10], 'D', label='Mayer et al. (2007)')
plt.plot([30, 5, 2], [210, 20, 10], '<', label='Giannuzzi et al. (2005)')
plt.plot([30, 10], [200, 100], '>', label='Kato et al. (1999)')


plt.xlabel("FIB voltage [keV]")
plt.ylabel("Si damage [Å]")
plt.legend(fontsize=10)
plt.show()

In [29]:

for key in recons.keys():
    roi = thickness_rois[key]
    (xi, xf), (yi, yf) = roi
    recon = recons[key][yi:yf, xi:xf]


    img_roi = recon[:,:,recon.shape[-1]//2]

    cols = column_positions[key][roi_indices[key][0][0]]
    cx, cy = cols.T
    cx -= yi
    cy -= xi


    plt.imshow(img_roi)
    for i in range(len(cx)):
        plt.text(cy[i], cx[i], str(i))
    plt.show()

    # plt.plot(*cols[43].T[::-1], 'r.')

In [23]:

from scipy.stats import norm
from scipy.special import erf

gaussian = lambda x, mu, sig : norm.pdf(x, loc=mu, scale=sig)
double_erf = lambda x, p, q, x0, x1 : p * (erf(q*(x - x0)) - erf(q*(x - x1)))

In [41]:

x = np.linspace(0, 1, 1000)
ps = np.linspace(50/300, 250/300, 15)
ys = []

for p in ps:
    y = gaussian(x, p, 10/300)
    ys.append(y)
ys = np.array(ys)

Y = ys.sum(axis=0)


fig, axs = plt.subplots(2, 2, figsize=(10, 7), dpi=300)


axs[0,0].plot(ps, np.zeros_like(ps), 'ko')
axs[0,0].set_xlim(0, 1)
axs[0,0].set_ylim(-1, 30)
axs[0,0].plot(x, Y, 'r-')
for i in [0, 1, 2, -1]:
    axs[0,0].plot(x, ys[i], 'k-')
    axs[0,0].plot([ps[i]]*2, [0, ys[i].max()], 'k--', lw=1)
axs[0,0].set_xticks([0, 1])
axs[0,0].set_yticks([])
axs[0,0].set_xlabel("Normalized distance")

fig.tight_layout()


c = column_lists['2V'][65]
z = np.arange(len(c)) * DZ['2V']
axs[0,1].plot(z, c, 'k.')

c = column_lists['8V'][25]
z = np.arange(len(c)) * DZ['8V']
axs[1,0].plot(z, c, 'k.')

c = column_lists['30V'][62]
z = np.arange(len(c)) * DZ['30V']
axs[1,1].plot(z, c, 'k.')

plt.show()

In [28]:

# from scipy.ndimage import map_coordinates

# arr = np.angle(recon)

# x1, x2 = 79, 79
# y1, y2 = 31, 55

# # Number of samples along the line
# n = int(np.hypot(x2 - x1, y2 - y1)) + 1

# y = np.linspace(x1, x2, n)
# x = np.linspace(y1, y2, n)
# z = np.arange(arr.shape[2])

# # Sample arr[x, y, layer] along the line
# xx, zz = np.meshgrid(x, z, indexing="xy")
# yy, _ = np.meshgrid(y, z, indexing="xy")

# cross_section = map_coordinates(
#     arr,
#     [xx, yy, zz],
#     order=1,
#     mode="nearest"
# )

# plt.figure(figsize=(2, 8))
# plt.imshow(
#     cross_section,
#     # origin="upper",
#     aspect="auto",
#     # extent=[0, np.hypot(x2 - x1, y2 - y1), arr.shape[2] - 1, 0]
# )
# plt.xlabel("Distance along line")
# plt.ylabel("Layer")
# # plt.colorbar()
# plt.show()

In [226]:

vmin = recon.min()
vmax = recon.max()

fig = plt.figure(figsize=(10, 6), dpi=300)

outer = fig.add_gridspec(
    2, 1,
    height_ratios=[1, 1],
    # hspace=0.3
)

# Top row: 4:1
gs_top = outer[0].subgridspec(1, 3, width_ratios=[6, 1, 1], wspace=0.25)

# Bottom row: 1:1:1 + colorbar
gs_bot = outer[1].subgridspec(
    1, 5,
    width_ratios=[1, 1, 1, 0.05, 0.1],
)

ax_blank = fig.add_subplot(gs_top[0])
ax_blank.axis("off")

ax_cs = fig.add_subplot(gs_top[1])
ax_cs.imshow(cross_section, cmap="magma",
             vmin=vmin, vmax=vmax, aspect="auto")
# ax_cs.axis("off")
# ax_cs.set_xlabel('y')
ax_cs.set_ylabel('Depth [Å]')
ax_cs.set_yticks(np.arange(0, 201, 50) / DZ[key])
ax_cs.set_yticklabels(np.arange(0, 201, 50))
ax_cs.set_xticks([])

xlim = ax_cs.get_xlim()
ax_cs.set_xlim(xlim)
ylim = ax_cs.get_ylim()

ax_cs.plot(xlim, [0, 0], 'r-')
ax_cs.plot(xlim, [2, 2], 'c-')
ax_cs.plot(xlim, [10, 10], 'g-')
# ax_cs.plot(xlim, [ylim[0]]*2, 'y-', lw=5)

ax_p = fig.add_subplot(gs_top[2])
py = np.array(column_lists[key])[roi_idx].squeeze()[66]
px = np.arange(len(py))
ax_p.sharey(ax_cs)
ax_p.plot(py, px, 'k.')
ax_p.set_xticks([0, 0.5])

xlim = ax_p.get_xlim()
ax_p.set_xlim(xlim)
ylim = ax_p.get_ylim()

t = np.linspace(*ylim)
ax_p.plot(double_erf(t, *popt), t, color=plt.colormaps['magma'].colors[2 * 256//3])


ax_p.plot(xlim, [0, 0], 'r-')
ax_p.plot(xlim, [2, 2], 'c-')
ax_p.plot(xlim, [10, 10], 'g-')

ax_p.set_xlabel('Phase [rad]')

axs = [fig.add_subplot(gs_bot[i]) for i in range(3)]
cax = fig.add_subplot(gs_bot[3])

im = axs[0].imshow(recon[:, :, 0], cmap="magma", vmin=vmin, vmax=vmax)
axs[1].imshow(recon[:, :, 2], cmap="magma", vmin=vmin, vmax=vmax)
axs[2].imshow(recon[:, :, 10], cmap="magma", vmin=vmin, vmax=vmax)

x1, x2 = axs[0].get_xlim()
y1, y2 = axs[0].get_ylim()

axs[0].plot([x1, x1, x2, x2, x1], [y1, y2, y2, y1, y1], 'r-', lw=5)
axs[1].plot([x1, x1, x2, x2, x1], [y1, y2, y2, y1, y1], 'c-', lw=5)
axs[2].plot([x1, x1, x2, x2, x1], [y1, y2, y2, y1, y1], 'g-', lw=5)

for ax in axs:
    ax.plot([73+79-0.5, 73+79-0.5], [73+31, 73+55], 'y-', lw=3)
    ax.axis("off")

fig.colorbar(im, cax=cax)
cax.set_ylabel('Phase [rad]')

plt.tight_layout()
plt.show()