%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
%matplotlib inline
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()
data_width = 5.04e-06 m data_height = 5.04e-06 m data_step = 1.96875e-08 m/pixel
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()
%matplotlib qt
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()
%matplotlib inline
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)')
<matplotlib.text.Text at 0x7f54a68660b8>
#'''
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
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()
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)
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')
(0, 255, 255, 0)
'''k = 404
mask_boundary[k] = False
mask_center[k] = True'''
#376,381,397,408,404
edge_angles, grain_angles = ts.FindAngles(x, y, edges, vertex_neighbours)
#=================================================================================================================================
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))
(0, 255, 255, 0)
plt.hist(grain_angles)
plt.title('grain_angles', fontsize=20)
plt.xlabel('angle')
plt.show()
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)
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))
(0, 255, 255, 0)
out_edges = ts.FindOutEdges(edges)
label_data_prec = ts.PreciseLabels(raw_data.shape, argmin_xy, out_edges, mask_center)
Calculated: 0%... 5%... 10%... 15%... 20%... 25%... 30%... 35%... 40%... 45%... 50%... 55%... 60%... 65%... 70%... 75%... 80%... 85%... 90%... 95%... 100%
# ===========================================================================================
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)')
<matplotlib.text.Text at 0x7f54bc057e10>
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()
height_half = height_mins + height_range/2
grain_area, grain_area_half = ts.FindDiameters(data_MinusND, label_data_prec, height_half, data_step)
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()
%matplotlib qt
grain_density_list = []
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')
grain_density_list.append(grain_density)
print('grain_density = {:.2f} um-2'.format(grain_density * data_unit**2))
grain_density = 27.45 um-2
grain_density_list
[30457166247780.469, 32251714801570.062, 28789176969749.867, 27447384255553.895]
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()
area_total = np.sum(data_findarea == 1) * data_step**2
print('area_total = ', area_total)
area_total = 2.16260112305e-11
%matplotlib inline
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()
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()