Perform a simulation on a tilted detector using Odin.

fit_ellipse.py

"""
Some simple functions TJ was using to fit ellipses. May be
useful, may not be.
"""


import numpy as np
from scipy import optimize

def polar_to_cart(polar_points):
    cart_points = np.zeros_like(polar_points)
    cart_points[:,0] = polar_points[:,0] * np.cos(polar_points[:,1])
    cart_points[:,1] = polar_points[:,0] * np.sin(polar_points[:,1])
    return cart_points

    
def cart_to_polar(cart_points):
    polar_points = np.zeros_like(cart_points)
    polar_points[:,0] = np.sqrt( np.sum( np.power( cart_points, 2 ), axis=1) )
    polar_points[:,1] = np.arctan2( cart_points[:,1], cart_points[:,0] )
    return polar_points


def fit_ellipse(points, a0=1.0, b0=1.0, x00=0.0, y00=0.0, beta0=0.0):
    """
    Fit an set of points in 2D to the equation for an ellipse by least-squares 
    regression.
    
    Parameters
    ----------
    points : np.ndarray
        An N x 2 array describing the (x,y) coordinates of the points to fit.
        
    Optional Parameters
    -------------------
    Correspond to the outputs -- they are initial guesses.
        
    Returns
    -------
    a, b : float
        The major and minor axis parameters of the ellipse
        
    x_0, y_0 : float
        The center (x_0, y_0) of the ellipse.
        
    beta : float
        The angle of rotation from the x-axis to the ellipse major axis.
    """
    
    if not ( (points.shape[1] == 2) and (len(points.shape) == 2) ):
        raise ValueError('mishapen `points`')
    
    def objective(ellipse_params):
        a, b, x_0, y_0, beta = ellipse_params
        polar_points = cart_to_polar( points - np.array([x_0, y_0]) )
        r_ellipse = evaluate_ellipse(polar_points[:,1], a, b, beta)
        differences = r_ellipse - polar_points[:,0]
        return differences
    
    x0 = np.array([a0, b0, x00, y00, beta0])
    opt, err = optimize.leastsq(objective, x0)
    a, b, x_0, y_0, beta = opt
    
    return a, b, x_0, y_0, beta


def evaluate_ellipse(theta, a, b, beta):
    """
    Evaluate r(theta) for an ellipse parameterized by (a, b, beta).
    
    Parameters
    ----------
    theta : float (or np.ndarray of floats)
        The angle with respect to the x-axis, in radians.
        
    a, b : float
        The major and minor axis parameters of the ellipse
        
    x_0, y_0 : float
        The center (x_0, y_0) of the ellipse.
        
    beta : float
        The angle of rotation from the x-axis to the ellipse major axis.
        
    Returns
    -------
    r : float (or np.ndarray of floats)
        The radius for each `theta`.
    """
    angle = theta + beta
    r = a*b / np.sqrt( np.power(b*np.cos(angle), 2) + np.power(a*np.sin(angle), 2) )
    return r

tilted_simulation.py

#!/usr/bin/env/python

"""
A simple script to perform a simulation on a tilted detector 
using Odin. The tilt is applied in the slow-scan (y-axis)
direction. Note that small tilts can make a big difference
in the simulated pattern!

The output is a rings file, containing the simulated stuff.

-- TJL May 14 2013
"""

import numpy as np
from odin import xray
from odin import structure
from odin.math2 import ER_rotation_matrix

# --------------------------------------------------------- #
# Manually set parameters -- change these                   #

traj = structure.load_coor('gold1k.coor') # or whatever

shp           = (201,201) # choose number of pixels
distance      = 20.0      # detector distance, pixel units
energy        = 10.0      # keV

detector_tilt = 0.0       # in degrees

num_molecules = 10
num_shots     = 1

q_values      = [2.66, 3.07] # inv. Ang.
num_phi       = 3600

output_fn     = 'tilted.ring'

# --------------------------------------------------------- #

s = np.array([0.0, 1.0, 0.0])
f = np.array([1.0, 0.0, 0.0])
c = np.array([0.0, 0.0, distance])

# rotate s backwards `detector_tilt` degrees
R = ER_rotation_matrix(f, np.radians(detector_tilt))
s = np.dot(R, s)
print "Rotated s-vector to:", s

bg = xray.BasisGrid()
bg.add_grid_using_center(c, s, f, shp)
b = xray.Beam(None, energy=energy)
d = xray.Detector(bg, b)

print "Simulating onto tilted detector..."
ss = xray.Shotset.simulate(traj, d, num_molecules, num_shots)
r = ss.to_rings(q_values, num_phi)

r.save(output_fn)
print 'Wrote: %s' % output_fn