Reads in an HRTEM image of pyroxene and allows one to rotate and draw an intensity plot along the a-axis to quantify ortho/clino/half plane stackings.
# coding: utf-8
# Created 2014, Zack Gainsforth
# Definitions:
# CLEN: CLinoENstatite, unit cell a = 9 Å
# OREN: ORthoENstatite, unit cell a = 18 Å
# HALF: Half unit cell, a = 4.5 Å
# PEN: ProtoENstatite, unit cell a = 18 Å
# An OREN -> CLEN transition always has to transform an entire 18 Å cell. Hence always even numbers of CLEN cells.
# For PEN origin, you can get half-planes. Don't know why.
import matplotlib
matplotlib.use('Qt4Agg')
import matplotlib.pyplot as plt
import numpy as np
import pyximport; pyximport.install()
from scipy import ndimage
from PIL import Image
### INPUTS SECTION
# Open the tif of the HRTEM image
FileName = '/Volumes/Stardust/Desktop/TEM Cecil unprocessed/20131219 - Libra - Cecil A6 S2/Cecil Grid A6 S2 - 0087-89 SUM.tif'
# If the scale on the image is off a little bit, then correct it.
ScaleCorr = 0.9651 / ((4.41 - 0.339)/5) # (22.399/22.3181)
ScaleCorr = 1.18
print 'ScaleCorr = ' + str(ScaleCorr)
# Stacking sequence. To start with, everything is a CLEN.
# OREN has a double unit cell, and some half cells can occur for a PEN origin.
CLEN = 0.9651
HALF = CLEN/2 # 0.51 #0.4405
OREN = 1.9 #1.8233
Sequence = np.ones(30)*CLEN
Sequence[[6, 12, 15, 17, 20, 24]] = HALF
Sequence[[25]] = OREN
# We'll draw lines at regular intervals starting at x0.
x0 = 1.907 #0.339
### CODE SECTION
def QuickPlot(x, y, xlabel=None, ylabel=None, title=None, legendstrs=None, savefile=None, boldlevel=1, figax=None):
"""
:param x: Can be a 1D numpy array, or a 2D numpy array. If 1D, then there is one plot. If 2D, every row is the x-axis for another line on the plot.
:param y: Identical dimension to x.
:param xlabel: String (can include TeX) for the label of the x-axis.
:param ylabel: String (can include TeX) for the label of the y-axis.
:param title: String (can include TeX) for the label of the title.
:param legendstrs: List of strings with the same length as x and y, axis 1.
:param savefile: A filename to save to if desired.
:param boldlevel: Sets the linewidths fatter and texts bigger. 1 is normal, 4 is usually good for publication.
:param figax: None makes a new figure and axis. Or plots on an existing one if a tuple is given: (fig, ax)
:return: (fig, ax) The figure and axis so the user can do stuff of his own with it if he wants to.
"""
# Set the bold level.
FontSizeBasis = (boldlevel+2)*4 # Fonts get bigger as boldlevel increases
TickMajorBasis = boldlevel*4 # As fonts get bigger, they need a larger padding from the axis.
# Increase the size of the tick label fonts.
matplotlib.rc('xtick', labelsize=FontSizeBasis)
matplotlib.rc('ytick', labelsize=FontSizeBasis)
# Increase their padding.
matplotlib.rc('xtick.major', pad=TickMajorBasis)
matplotlib.rc('ytick.major', pad=TickMajorBasis)
# Create/reuse a figure to plot on.
if figax is None:
fig, ax = plt.subplots()
else:
fig, ax = figax
# Figure out how many lines are going on this plot.
dim = np.shape(x)
if len(dim) == 1:
# Just one line on this plot.
ax.plot(x,y, linewidth=boldlevel)
elif len(dim) == 2:
# Multiple lines on this plot.
for i in range(dim[0]):
ax.plot(x[i,:], y[i,:], linewidth=boldlevel)
else:
raise TypeError('x is not a 1D or 2D numpy array.')
# Write the x and y labels and the title.
if xlabel is not None:
plt.xlabel(xlabel, fontsize=FontSizeBasis)
if ylabel is not None:
plt.ylabel(ylabel, fontsize=FontSizeBasis)
if title is not None:
plt.title(title, fontsize=FontSizeBasis)
# Show a legend if input.
if legendstrs is not None:
ax.legend(legendstrs)
# Resize the figure appropriately
plt.tight_layout()
# Save it to a file if a name is given.
if savefile is not None:
plt.savefig(savefile)
return fig, ax
def GetLinePlotx(Image, Startx, y, width, height):
""" Given a numpy 2D array representing an image, this creates an integrated (along the y-axis)
linescan from Startx to Endx. It returns the linescan and the section of the image that was integrated.
"""
z = np.zeros(width)
Endx = Startx + width
# Height has to be an odd number (middle plus two sides)
if not height % 2:
height += 1
for i in range(len(z)):
z[i] = sum(Image[y-height/2:y+height/2, Startx+i])
SubImage = Image[y-height/2:y+height/2, Startx:Endx]
return z, SubImage
# Open the trace of the HRTEM image.
# FileName = '/Volumes/Stardust/Desktop/TEM Cecil unprocessed/20131219 - Libra - Cecil A6 S2/Cecil Grid A6 S2 - 0087-89 SUM Trace.txt'
# d = np.genfromtxt(FileName)
# FileName = '/Volumes/Stardust/Desktop/TEM Cecil unprocessed/20131219 - Libra - Cecil A6 S2/testimage.tif'
img = plt.imread(FileName)
imgtags = Image.open(FileName)
# TIFF tag 282 is the number of pixels/unit
dx = imgtags.tag.tags[282]
dx = float(dx[0][1])/float(dx[0][0])
# Now dx is units/pixel. We will assume units are nm.
#print dx
# Rotate it, plot it and draw some vertical lines so we can adjust the rotation iteratively.
imgrot = ndimage.rotate(img, -127.15)
plt.gray()
plt.imshow(imgrot)
plt.plot([800, 800], [0,np.shape(imgrot)[0]], 'r')
plt.plot([1500, 1500], [0,np.shape(imgrot)[0]], 'r')
plt.gcf().set_size_inches(15, 15, forward=True)
plt.gcf().canvas.toolbar.pan()
# Get the linescan.
z, SubImage = GetLinePlotx(imgrot, 435, 1150, 1000, 100)
# Make an x-axis to go with the extracted line plot.
x = np.array(range(len(z)))*dx
# Apply the scale correction if our image scale is a little off.
x *= ScaleCorr
# Now make a figure with two subplots. The top has our linescan. The bottom
# has the subimage from which the linescan was derived.
fig = plt.figure()
ax1 = plt.subplot(2,1,1)
QuickPlot(x,z, xlabel='nm', boldlevel=2, figax=(fig,ax1))
ax2 = plt.subplot(2,1,2)
plt.imshow(SubImage)
# Lines will be drawn vertically from 0 up to the maximum intensity in the plot.
linetop = max(z)
# And loop through the sequence.
xnow = x0
for i in range(len(Sequence)):
# Choose the color for the line.
if Sequence[i] == CLEN:
c = 'r'
elif Sequence[i] == OREN:
c = 'g'
else:
c = 'k'
# And draw it.
ax1.plot([xnow, xnow], [0, linetop], c)
# Add an index so we can type in the new sequence
ax1.text(xnow+0.05, 0, i)
# Increment our x position depending on what width this cell was.
xnow += Sequence[i]
fig.set_size_inches(24, 6, forward=True)
ax1.set_xlim(min(x),max(x))
plt.draw()
plt.show()