In [1]:
%matplotlib inline
import numpy as np
import scipy as sp
from scipy import ndimage, spatial
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from matplotlib.widgets import Cursor, MultiCursor, RectangleSelector, Slider
from matplotlib.patches import Rectangle
import gwyddion_import
import tesselect as ts
In [94]:
%matplotlib inline
In [88]:
file_path = 'data/'
file_name = '2013-01-16_№42_l.clossiana.mdt-h15.txt'
raw_data, data_width, data_height, data_step, data_unit = gwyddion_import.ReadData(file_path + file_name)
data_MinusND = ts.MinusND(raw_data, 7, 'hv')
data_bin = (data_MinusND > 1.1*np.mean(data_MinusND)) # coefficient = 1.1 is up to you
#=================================================================================================================================
print('data_width =', data_width, 'm')
print('data_height =', data_height, 'm')
print('data_step =', data_step, 'm/pixel')
plt.figure(figsize=(18, 9), dpi=80)
plt.subplot(121), plt.imshow(data_MinusND, cmap='gray')
plt.title('data_MinusND', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(122), plt.imshow(data_bin, cmap='gray')
plt.title('data_bin', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.show()
In [95]:
data_bin_op = ndimage.binary_opening(data_bin, iterations=1)
#=================================================================================================================================
plt.figure(figsize=(18, 12), dpi=80)
plt.subplot(231), plt.imshow(ndimage.binary_opening(data_bin, iterations = 1), cmap='gray')
plt.title('binary opening (1 iteration)', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(232), plt.imshow(ndimage.binary_opening(data_bin, iterations = 2), cmap='gray')
plt.title('binary opening (2 iterations)', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(233), plt.imshow(ndimage.binary_opening(data_bin, iterations = 3), cmap='gray')
plt.title('binary opening (3 iterations)', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(234), plt.imshow(data_MinusND, cmap='gray')
plt.title('data_MinusND', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(235), plt.imshow(data_bin, cmap='gray')
plt.title('data_bin', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(236), plt.imshow(data_bin_op, cmap='gray')
plt.title('data_bin_op (1 iteration)', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.show()
In [96]:
%matplotlib qt
In [97]:
data_bin_good = data_bin_op.copy()
#data_bin_op = data_bin_good.copy()
fig, ax = plt.subplots(1, 2, figsize=(30,15))
cursor = MultiCursor(fig.canvas, ax, color='y', linewidth=1, horizOn=True, vertOn=True)
ax[0].imshow(data_MinusND, cmap='gray')
im1 = ax[1].imshow(data_bin_op, cmap='gray', interpolation='nearest')
numrows, numcols = raw_data.shape
brush_size = 1
# Coordinates formatting =========================================
def format_coord(x, y):
col = int(x + 0.5)
row = int(y + 0.5)
if col >= 0 and col < numcols and row >= 0 and row < numrows:
z = data_bin_good[row,col]
return 'x = %d, y = %d, z = %d'%(x + 0.5, y + 0.5, z)
ax[0].format_coord = format_coord
ax[1].format_coord = format_coord
# Mouse events ====================================================
def mouse_draw(event, col, row):
global brush_size
num = int((brush_size - 1)/2)
if event.button == 1:
data_bin_good[row-num:row+num+1,col-num:col+num+1] = False
elif event.button == 3:
data_bin_good[row-num:row+num+1,col-num:col+num+1] = True
im1.set_data(data_bin_good)
plt.draw()
def on_mouse_press(event):
if brush_size == 0: return
col = int(event.xdata + 0.5)
row = int(event.ydata + 0.5)
if col >= 0 and col < numcols and row >= 0 and row < numrows:
mouse_draw(event, col, row)
cid_mouse_move = fig.canvas.mpl_connect('motion_notify_event', on_mouse_move)
def on_mouse_move(event):
if brush_size == 0: return
col = int(event.xdata + 0.5)
row = int(event.ydata + 0.5)
if col >= 0 and col < numcols and row >= 0 and row < numrows:
mouse_draw(event, col, row)
# Slider events =================================================
def brush_change(val):
global brush_size
if int(val) == 0:
brush_size = int(val)
elif int(val) % 2 == 0:
brush_size = int(val)
else:
brush_size = int(val)
slider_ax = plt.axes([0.12, 0.05, 0.78, 0.02])
slider = Slider(slider_ax, "Brush Width", 0, 25, valinit=1, color='b', valfmt='%0.2f')
slider.on_changed(brush_change)
# The Beginning =================================================
cid_mouse_press = fig.canvas.mpl_connect('button_press_event', on_mouse_press)
fig.canvas.mpl_disconnect(fig.canvas.manager.key_press_handler_id)
plt.show()
In [98]:
%matplotlib inline
In [99]:
plt.imsave(file_path + file_name + '.png', data_bin_good)
plt.figure(figsize=(10, 10), dpi=80)
plt.imshow(data_bin_good, cmap='gray')
plt.title('data_bin_good', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
Out[99]:
In [147]:
#'''
data_bin_good = plt.imread(file_path + file_name + '.png')[:,:,0]
plt.figure(figsize=(10, 10), dpi=40)
plt.imshow(data_bin_good, cmap='gray')
plt.title('data_bin_good', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
#'''
pass
In [148]:
label_data, num_of_labels = ndimage.label(data_bin_good)
data_centers = ndimage.measurements.center_of_mass(data_bin_good, label_data, range(1, num_of_labels+1)) # list of tuples
data_centers = np.rint(data_centers).astype(int)
x = data_centers[:, 1]
y = data_centers[:, 0]
data_centers_pic = np.zeros(data_bin_good.shape)
data_centers_pic[y, x] = 1 # because plt.imshow() reverts axis
#=================================================================================================================================
plt.figure(figsize=(20, 20), dpi=80)
plt.subplot(221), plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray')
plt.title('data_bin_good and data_centers', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(222), plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray')
plt.title('data_bin_good (numerated)', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
for i in range(len(x)):
plt.annotate(str(i), (x[i]-2, y[i]-1), color='c')
plt.subplot(223), plt.imshow(data_MinusND, cmap='gray')
plt.title('data_MinusND', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.subplot(224), plt.imshow(data_MinusND, cmap='gray'), plt.plot(x, y, 'o')
plt.title('data_MinusND and data_centers', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
plt.show()
In [149]:
tri = ts.FindTriangulation(x, y)
edges, edge_lengths, coordination_numbers, vertex_neighbours = ts.TriEdges(tri)
mask_center, mask_boundary = ts.BoundaryVertices(coordination_numbers, vertex_neighbours)
In [150]:
plt.figure(figsize=(20,20), dpi=80)
# -----------------------------------------------------------------------------------------------
plt.subplot(221)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(mtri.Triangulation(x, y), 'g-')
plt.plot(x, y, 'ro')
plt.title('Not Filtered triangulation', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
# -----------------------------------------------------------------------------------------------
plt.subplot(222)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(tri, 'g-')
plt.plot(x, y, 'ro')
plt.title('Filtered triangulation', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
#for i in range(len(data_centers)):
# plt.annotate(str(i), (x[i]+2, y[i]), color='c')
# -----------------------------------------------------------------------------------------------
plt.subplot(223)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(tri, 'g-')
plt.plot(x[mask_center], y[mask_center], 'ro')
plt.plot(x[mask_boundary], y[mask_boundary], 'bo')
plt.title('Filtered triangulation with boundary', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
#for i in range(len(data_centers)):
# plt.annotate(str(i), (x[i]+2, y[i]), color='c')
# -----------------------------------------------------------------------------------------------
plt.subplot(224)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(tri, 'g-')
mask = np.logical_and((coordination_numbers != 6), mask_center)
plt.plot(x, y, 'bo')
plt.plot(x[mask], y[mask], 'ro')
plt.title('Vertices with coordination number != 6', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
#for i in range(len(data_centers)):
# plt.annotate(str(i), (x[i]+2, y[i]), color='c')
Out[150]:
In [60]:
'''k = 404
mask_boundary[k] = False
mask_center[k] = True'''
#376,381,397,408,404
In [151]:
edge_angles, grain_angles = ts.FindAngles(x, y, edges, vertex_neighbours)
In [152]:
#=================================================================================================================================
plt.figure(figsize=(20,10), dpi=80)
plt.subplot(121)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(tri, 'g-')
for i in range(len(x)):
plt.plot(x[i], y[i], 'o', markersize=8, color=plt.cm.hsv(grain_angles[i]/60+0.5))
#plt.scatter(x, y, s=200, c=grain_angles, cmap='hsv')
#plt.colorbar(shrink=0.8)
#plt.clim(-30,30)
plt.axis((0,255,255,0))
plt.title('Grain angles distribution', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
# ---------------------------------------------------------------------------------------------------------------
plt.subplot(122)
plt.imshow(data_bin_good + 2*data_centers_pic, cmap='gray', interpolation='nearest')
plt.triplot(tri, 'g-')
mask = np.logical_and((coordination_numbers != 6), mask_center)
plt.plot(x, y, 'bo', markersize=8)
plt.plot(x[mask], y[mask], 'ro', markersize=8)
plt.title('Vertices with coordination number != 6', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
Out[152]:
In [153]:
plt.hist(grain_angles)
plt.title('grain_angles', fontsize=20)
plt.xlabel('angle')
plt.show()
In [154]:
label_max = ndimage.maximum(data_MinusND, label_data, range(1, num_of_labels + 1))
height_mins_tri, argmin_array, argmin_xy = ts.FindMinsTri(data_MinusND, label_data, x, y, edges)
In [155]:
plt.figure(figsize=(10, 10), dpi=80)
plt.imshow(data_MinusND, cmap='gray')
plt.plot(argmin_xy[0], argmin_xy[1], 'yo')
plt.title('data_MinusND and coordinates of mins', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
Out[155]:
In [156]:
out_edges = ts.FindOutEdges(edges)
label_data_prec = ts.PreciseLabels(raw_data.shape, argmin_xy, out_edges, mask_center)
In [157]:
# ===========================================================================================
plt.figure(figsize=(20, 20), dpi=80)
plt.subplot(221)
plt.imshow(label_data_prec > 0, cmap='gray')
plt.title('label_data_prec', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
plt.subplot(222)
plt.imshow(label_data_prec > 0, cmap='gray')
plt.plot(argmin_xy[0], argmin_xy[1], 'yo')
plt.title('label_data_prec and coordinates of mins', fontsize=20), plt.xlabel('(in pixels)'), plt.ylabel('(in pixels)')
plt.axis((0,255,255,0))
plt.subplot(223)
plt.imshow(data_bin_good, cmap='gray')
plt.title('data_bin_good', fontsize=20)
plt.xlabel('(in pixels)')
plt.ylabel('(in pixels)')
Out[157]:
In [158]:
label_prec_max = ndimage.maximum(data_MinusND, label_data_prec, range(1, num_of_labels + 1))
height_mins = np.asarray([np.mean(i) for i in height_mins_tri])
height_range = label_prec_max - height_mins
# ==============================================================================
plt.figure(figsize=(20,6))
plt.subplot(121)
plt.hist((label_max - height_mins)[mask_center])
plt.title('label_max - height_mins', fontsize=20)
plt.xlabel('(in meters)')
plt.subplot(122)
plt.hist((label_prec_max - height_mins)[mask_center])
plt.title('label_prec_max - height_mins', fontsize=20)
plt.xlabel('(in meters)')
plt.show()
In [159]:
height_half = height_mins + height_range/2
grain_area, grain_area_half = ts.FindDiameters(data_MinusND, label_data_prec, height_half, data_step)
In [160]:
plt.figure(figsize=(20,6))
plt.subplot(121)
plt.hist(np.sqrt(grain_area[mask_center]))
plt.title('sqrt(grain_area)', fontsize=20)
plt.xlabel('(in meters)')
plt.subplot(122)
plt.hist(np.sqrt(grain_area_half[mask_center]), bins=10)
plt.title('sqrt(grain_area_half)', fontsize=20)
plt.xlabel('(in meters)')
plt.show()
In [161]:
%matplotlib qt
grain_density_list = []
In [168]:
fig, ax = plt.subplots(figsize=(15,15))
cursor = Cursor(ax, useblit=True, color='y', linewidth=1)
plt.plot(x, y, 'ro')
ax.set_axis_bgcolor('k')
plt.axis((0,255,255,0))
def selector_callback(eclick, erelease):
global rect, grain_density, area_, x, y, grain_count
try:
rect.remove()
except:
pass
x1, x2 = min(eclick.xdata, erelease.xdata), max(eclick.xdata, erelease.xdata)
y1, y2 = min(eclick.ydata, erelease.ydata), max(eclick.ydata, erelease.ydata)
rect = Rectangle((x1,y1), width=(x2-x1), height=(y2-y1), color='yellow', alpha=0.4)
ax.add_patch(rect)
plt.draw()
area_ = (x2-x1)*(y2-y1)*data_step**2
grain_count = 0
for i in range(len(x)):
if (x1 < x[i] < x2) and (y1 < y[i] < y2):
grain_count += 1
grain_density = grain_count/area_
selector = RectangleSelector(ax, selector_callback, drawtype='box', useblit=True,
button=[1,3], minspanx=5, minspany=5, spancoords='pixels')
In [169]:
grain_density_list.append(grain_density)
print('grain_density = {:.2f} um-2'.format(grain_density * data_unit**2))
In [170]:
grain_density_list
Out[170]:
In [ ]:
In [171]:
fig, ax = plt.subplots(figsize=(30,30))
cursor = Cursor(ax, useblit=True, color='y', linewidth=1)
im1 = ax.imshow(data_MinusND, cmap='gray')
ax.plot(x, y, 'bo')
plt.axis((0,255,255,0))
numrows, numcols = raw_data.shape
brush_size = 1
data_findarea = data_MinusND.copy()
# Mouse events ====================================================
def mouse_draw(event, col, row):
global brush_size
num = int((brush_size - 1)/2)
if event.button == 1:
data_findarea[row-num:row+num+1,col-num:col+num+1] = False
elif event.button == 3:
data_findarea[row-num:row+num+1,col-num:col+num+1] = True
im1.set_data(data_findarea)
plt.draw()
def on_mouse_press(event):
if brush_size == 0: return
col = int(event.xdata + 0.5)
row = int(event.ydata + 0.5)
if col >= 0 and col < numcols and row >= 0 and row < numrows:
mouse_draw(event, col, row)
cid_mouse_move = fig.canvas.mpl_connect('motion_notify_event', on_mouse_move)
def on_mouse_move(event):
if brush_size == 0: return
col = int(event.xdata + 0.5)
row = int(event.ydata + 0.5)
if col >= 0 and col < numcols and row >= 0 and row < numrows:
mouse_draw(event, col, row)
# Slider events =================================================
def brush_change(val):
global brush_size
if int(val) == 0:
brush_size = int(val)
elif int(val) % 2 == 0:
brush_size = int(val)
else:
brush_size = int(val)
slider_ax = plt.axes([0.12, 0.05, 0.78, 0.02])
slider = Slider(slider_ax, "Brush Width", 0, 25, valinit=1, color='b', valfmt='%0.2f')
slider.on_changed(brush_change)
# The Beginning =================================================
cid_mouse_press = fig.canvas.mpl_connect('button_press_event', on_mouse_press)
fig.canvas.mpl_disconnect(fig.canvas.manager.key_press_handler_id)
plt.show()
In [172]:
area_total = np.sum(data_findarea == 1) * data_step**2
print('area_total = ', area_total)
In [173]:
%matplotlib inline
In [174]:
rdf_x, rdf_y = ts.RDF(x, y, data_step, area_total)
plt.figure(figsize=(8,5))
plt.plot(rdf_x*1e9, rdf_y, 'b-')
plt.plot(np.array([0,8000]), np.array([1,1]), 'r--')
plt.xlim([0,1500])
plt.xlabel('(in nm)')
plt.show()
In [175]:
data_fft1 = np.log(np.abs(np.fft.fftshift(np.fft.fft2(raw_data))))
data_fft2 = np.log(np.abs(np.fft.fftshift(np.fft.fft2(data_MinusND))))
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(data_fft1, cmap='gray')
plt.subplot(122)
plt.imshow(data_fft2, cmap='gray')
plt.show()
In [49]: