Sunday, May 13, 2012



Kikuchi bands from Silicon.

Saturday, May 12, 2012

A rotating cube...

I have been learning Python. I want to work towards writing some scripts that together take EBSD-based orientation data which are relative to a sample surface reference frame and rotate the data for visualization relative to a kinematic or geographic reference frame.
I decided to start simple – learn something about quaternions and rotate 2D plots of a cube.

Below is some code to create an animation (.avi file)  by rotating cube axes.
Practicing. Learning. Far from the goal.

NOTE: IF you want to use this code, you'll need mencoder for the script to successfully create the animation file.

A glimpse of a video that this script produced: cube movie


__________________________________________________________________________




#!/usr/bin/python
##########################################################################
## Cube animation program using python, numpy, os, and matplotlib.
## This script demonstrates how rotate plots using quaternion conversion.
## 
## Written by Zach Michels. Last Edited 5/12/2012.
## Cheers!
##########################################################################

from pylab import *
import numpy as np
import subprocess                 # For issuing commands to the OS.
import os, time, stat, Image
import sys                        # For determining the Python version.

##########################################################################
## Produces a rotation matrix from the 3 last components of a quaternion.
def quaternion_to_matrix(myx):
    # '''converts from a quaternion representation (last 3 values) to rotation matrix.'''
    xb,xc,xd = myx

    xnormsq = xb*xb+xc*xc+xd*xd

    if xnormsq < 1:
        ## Inside the unit sphere, these are the components,
        ## and we only have to calculate 'a' to normalize
        ## the quaternion.
        b,c,d = xb,xc,xd
        a = np.sqrt(1-xnormsq)
    else:
        ## If we have invalid inputs, we reflect the vectors
        ## outside the unit sphere to the other side, and
        ## use the inverse norm and negative values of 'a'.
        b,c,d = -xb/xnormsq,-xc/xnormsq,-xd/xnormsq
        a = -np.sqrt(11.0/xnormsq  )
        ## This generally defeats the whole concept that
        ## small (b,c,d) vector norms have small rotation angles. It's
        ## really just to work near the borders of the sphere. 
        ## Optimization algorithms should work with initalizations 
        ## just inside the spere, and avoid wandering outside of it.

    assert a >= -1
    assert a <= 1

    ## Notice we return a transpose matrix, because we work with line-vectors

    return np.array([ [(a*a+b*b-c*c-d*d), (2*b*c-2*a*d),     (2*b*d+2*a*c)      ],
                         [(2*b*c+2*a*d),     (a*a-b*b+c*c-d*d), (2*c*d-2*a*b)      ],
                         [(2*b*d-2*a*c),     (2*c*d+2*a*b),     (a*a-b*b-c*c+d*d)] ]  ).T  \
                         / (a*a+b*b+c*c+d*d)
                         

###########################################################################


## Calculate quaternion values using sinusoids over the frame index
## number, and get rotation matrix. Notice that the second parameter
## is the number of frames. That allows us to easily create an
## animation that can be played as a loop.
def crazy_rotation(ind,Nind):
    return quaternion_to_matrix(0.5*sin(pi*2*ind*Nind**-1.*array([1,2,3])))

## This function calculates the projected image coordinates form the 3D
## coordinates of the object vertices.
def project(D, vecs):
    return vvs[:,:2]/(vvs[:,[2,2]]-D)

## The cube vertices
vs = reshape(mgrid[-1:2:2,-1:2:2,-1:2:2].T, (8,3))

## Generate the list of connected vertices
ed=[(j,k)
    for j in range(8)
    for k in range(j,8)
    if sum(abs(vs[j]-vs[k]))==2 ]


## 'Camera' position
D=-5

## Create the figure. figsize needed to set the image sizes at the end...
figure(1, figsize=(6.4,4.8))

## create/get the axes object.
ax=subplot(1,1,1)

## Set image title.
suptitle("Zach's rotating cube")

## Number of frames. Lets user set number.

Nind = input("How many frames would you like to make in your animation? Each frame corresponds to a percent of the overall intended rotations, such that a greater number of frames ultimately produces a slower animation. This script will write a single .PNG image for each frame into the same directory as the .py file. These .PNG files are used to make the .AVI file, but are not automatically deleted. You will be prompted later for the directory PATH if you would like to delete them. NOTE: a greater number of frames will require longer time for script processing. NUMBER OF FRAMES (200+ recommended): ")


##
## Now start a loop over the animation frames. The index is used both
## to name the images and to calculate the rotation matrix. You could
## use this index to make some physics stuff...
for ind in range(Nind):
    print ind
    ## This is crucial, clears the figure for the new plot. 
    ax.clear() 

    ## Get the rotation matrix for the current frame.
    rotM = crazy_rotation(ind, Nind)

    ## Calculate the 3D coordinates of the vertices of the rotated
    ## cube. Just multiply the vectors by the rotation matrix...
    vvs=dot(vs,rotM)

    ## Now calculate the image coordinates of the points.
    pt = project(D,vvs)

    ## Plot the edges.
    for j,k in ed:
        ax.plot(pt[[j,k],0], pt[[j,k],1], 'm-', lw=3)

    ## Plot the vertices.
    ax.plot(pt[:,0], pt[:,1], 'bo')

    ## Set axes limits
    ax.axis('equal')
    ax.axis([-0.5,0.5,-0.5,0.5])

    ## Save the current frame. We need the dpi (along with the figure
    ## size up there) to set the image size.
    savefig('anim%03d.png'%ind, dpi=100)

# Now stitch graphed images together using Mencoder to create a movie.
# Each image will become a single frame in the movie.
#
# Use Python to make what would normally be a command line
# call to Mencoder.  Specifically, the command line call we want to
# emulate is (without the initial '#'):
# mencoder mf://*.png -mf type=png:w=800:h=600:fps=25 -ovc lavc -lavcopts 
# vcodec=mpeg4 -oac copy -o output.avi
# See the MPlayer and Mencoder documentation for details.
#

command = ('mencoder',
           'mf://*.png',
           '-mf',
           'type=png:w=800:h=600:fps=25',
           '-ovc',
           'lavc',
           '-lavcopts',
           'vcodec=mpeg4',
           '-oac',
           'copy',
           '-o',
           'cube_zm.avi')

# os.spawnvp(os.P_WAIT, 'mencoder', command)

print "\n\nabout to execute:\n%s\n\n" % ' '.join(command)
subprocess.check_call(command)

# explain compiled video output file and format
print "\n\n The movie was written to 'cube_zm.avi'"


print "\n\n You may want to delete *.png now.\n\n"

# opportunity to delete .PNG files used to create .AVI animation. 
# Note: If the directory PATH supplied is incorrect, the script 
# will return an error and the files will not be dleted. 
# However, the file cube.avi should still be created.

if raw_input("Delete all .PNG files in the directory? (y/n): ") == "y":
path = raw_input("Path to directory? (ex: '/Users/Usernname/folder'):  ")
if raw_input("are you sure you want to delete these files? Ths cannot be undone. (y/n): ") == "y":

dircount=0
pngcount=0

# deletes .PNG files by stepping through working directory 
# (the directory that contains the .py file).
for file in os.listdir(path):
    dircount=dircount+1
    basename=os.path.basename(file)
    NoExtName=os.path.splitext(file)[0]
   
    if basename.endswith('.png'):
pngcount=pngcount+1
print 'Deleting '+basename
os.remove(file)

Monday, May 7, 2012

Go Micro.