Tillmann,
I've added a jiffy script to synthesize pseudo-precession photos from a
rotation dataset, to the latest PHENIX package. We build this package
nightly, so any PHENIX bundle with a version number greater than
"dev-402" will work....however I see that dev-402 is not yet on the
public Web site (http://www.phenix-online.org) so I'll attach the python
script to this email for standalone use.
Here's what you'll need:
1) Make sure any recent PHENIX is installed and on path.
2) Make sure Mosflm is on path as "ipmosflm"
3) Index using two frames in the dataset:
labelit.index /data/project/dataset1_1_E1_###.img 1 90
4) Generate the precession photo:
labelit.python precession_photo.py bravais_choice=1 \ #id number for
the bravais setting as listed by labelit.index
pixel_width=600 resolution_outer=3.0 intensity_full_scale=1000 \
plot_section="H,K,0" image_range=1,90 pdf_output.file="HKL.pdf"
5) These command line options should be self explanatory, but "--help"
gives a fuller description.
Nick Sauter
Tillmann Heinisch wrote:
Hi,
I have problems solving the structure of a protein crystal which seems to be
disordered. In order to investigate the disorder it would be useful to have a
precision photograph that shows reflections only in the [0kl] plane. Does
anyone know software that can transform raw data to give intensity distribution
in distinct zones of hkl?
Many Thanks for your help,
Tillmann
def usage():
return """Python script to generate a synthetic photograph.
Requirements: Phenix must be installed and on path.
Mosflm must be installed and "ipmosflm" on path.
Usage:
1. First index the dataset using labelit.index. Typical syntax to index two
images is
labelit.index /home/data/project/dataset1_1_E1_###.img 1 90
2. Then
labelit.precession_photo <command line arguments>
Example:
labelit.precession_photo bravais_choice=1 \\
pixel_width=600 \\
resolution_outer=3.0 \\
intensity_full_scale=1000 \\
plot_section="H,K,0" \\
image_range=1,90 \\
pdf_output.file="HKL.pdf"
3. Command line arguments explained:
labelit.precession_photo --help
Known limitations:
a) Since the source images are rotation photographs, there is an inherent
approximation in assuming that the data was obtained at the central
rotation angle (0.5 degree in a 0.0 to 1.0 rotation).
b) Real precession photographs are obtained with a finite-width annulus;
however the approximation here is one of an infinitesmal-width opening.
c) This is a naive implementation that reads all of the data into memory
first before calculating the image. More efficient algorithms will be
implemented later.
Author: Nick Sauter, LBNL; May 8, 2010.
Released under the CCTBX license; see cctbx.sf.net.
"""
import math,pickle,os
import libtbx.phil
precession_master_phil = libtbx.phil.parse("""
labelit_precession {
pixel_width = 1500
.type = int
.help = "Width of pixel grid on which to calculate the synthesized image;
must be sufficiently large to sample the Bragg spots"
resolution_outer = 4.0
.type = float
.help = "The high-resolution limit (Angstroms) for the synthesized image"
bravais_choice = None
.type = int
.help = "Crystal setting (integer ID#), as enumerated in labelit.index
output, chosen for the synthesized section; print a list with
labelit.stats_index"
plot_section = None
.type = str
.help = "Miller indices of the section to be synthesized, comma separated,
no spaces, such as H,0,L"
image_range = None
.type = ints(size=2)
.help = "Range of data frames (inclusive) to be used for synthesized image,
i.e., min frame no.,max frame no., no spaces"
intensity_full_scale = 255
.type = int
.help = "Image intensity to be used as full scale (black)"
pdf_output {
file = None
.type = str
.help = "File name for pdf output"
}
}
""")
from scitbx import matrix
from scitbx.array_family import flex
from cctbx.crystal_orientation import crystal_orientation
from cctbx.uctbx import unit_cell
from spotfinder.diffraction.imagefiles import
Spotspickle_argument_module,image_files
from labelit.dptbx import AutoIndexEngine, Parameters
from labelit.command_line.imagefiles import QuickImage
from annlib_ext import AnnAdaptor
def get_precession_parameters(sources=[]):
parameters = precession_master_phil.fetch(sources=sources)
#validation not done here
return parameters
def ai_factory(inputpd):
subpd = inputpd["best_integration"]
ai =
AutoIndexEngine(inputpd["endstation"],inputpd["recommended_grid_sampling"])
P = Parameters(xbeam=float(inputpd["xbeam"]),ybeam=float(inputpd["ybeam"]),
distance=float(inputpd["distance"]),twotheta=float(inputpd["twotheta"]))
ai.setBase(P)
ai.setWavelength(float(inputpd['wavelength']))
ai.setMaxcell(float(inputpd['ref_maxcel']))
ai.setDeltaphi(float(inputpd['deltaphi'])*math.pi/180.)
ai.setMosaicity(subpd["mosaicity"])
ai.setOrientation(subpd["orient"])
return ai
def get_col(n,sqr):
return matrix.col((sqr[n],sqr[n+3],sqr[n+6]))
class calculate_miller:
def __init__(self,plot_section,ori):
labels = "HKL"
hkl_list = plot_section.split(",")
self.section_number = None
for trial_plot_normal_index in [0,1,2]:
try:
self.section_number = int(hkl_list[trial_plot_normal_index])
self.idx1 = (trial_plot_normal_index+1)%3
assert hkl_list[self.idx1]==labels[self.idx1]
self.idx2 = (trial_plot_normal_index+2)%3
assert hkl_list[self.idx2]==labels[self.idx2]
self.idx0 = (trial_plot_normal_index+3)%3
break
except:
pass # one of three settings must succeed
assert self.section_number!=None
assert self.idx0!=None and self.idx1!=None and self.idx2!=None
print self.idx0,self.idx1,self.idx2
rmat = matrix.sqr(ori.reciprocal_matrix())
self.unit_recip_basis1 = get_col(self.idx1,rmat).normalize() #like b*-hat
self.unit_recip_basis2 = get_col(self.idx2,rmat).normalize() #like c*-hat
self.a_hat_perp =
self.unit_recip_basis1.cross(self.unit_recip_basis2).normalize()
# normal to the plotted plane, like b*-hat x c*-hat
self.c_hat_perp = self.a_hat_perp.cross(self.unit_recip_basis1)
# second direction of plot, like a-hat-perp x b*-hat
self.recip_basis1 = get_col(self.idx1,rmat)
self.recip_basis2 = get_col(self.idx2,rmat)
self.recip_basis3 = get_col(self.idx0,rmat)
#reciprocal angle between recip_basis1 & recip_basis2
angle = math.acos(self.unit_recip_basis1.dot(
self.unit_recip_basis2))
#re-expression matrix; (f1.e1 f2.e1)
# (f1.e2 f2.e2)
self.reexpress =
matrix.sqr((self.recip_basis1.length(),self.recip_basis2.length()*math.cos(angle),0.,
0.,
self.recip_basis2.length()*math.sin(angle),0.,
0.,0.,1.)
).inverse()
def set_dmax(self,dmax):
self.dmax=dmax
class precession_plot:
def HKL_from_xy(self,x,y):
# assumes x,y are page coordinates, origin at top left, x down (slow), y
left to right (fast)
# reciprocal space starts in the middle of page, axis1 to the right, axis2
up
width = float(self.P.labelit_precession.pixel_width)/2.0
coord1 = ((float(y) - width)/width)*(1./self.miller_calc.dmax)
coord2 = ((-float(x) + width)/width)*(1./self.miller_calc.dmax)
#plotting vector, represented in the orthonormal basis e = (b*-hat,
c-hat-perp)
r = matrix.col((coord1,coord2,0.0))
#re-express in the crystal basis, f = (b*, c*)
miller = self.miller_calc.reexpress * r
miller3 = self.miller_calc.section_number
hkl_list = [0.,0.,0.]
hkl_list[self.miller_calc.idx0]=miller3
#print hkl_list
#hkl_zero = self.projection_onto_section_zero(hkl_list)
hkl_list[self.miller_calc.idx1]=miller[0]
hkl_list[self.miller_calc.idx2]=miller[1]
#correct for the zero of this section
hkl_list = matrix.col(hkl_list) - self.hkl_zero
#print "hkl_list",hkl_list
if x%50==0 and y%50==0:
print "%6d%6d"%(x,y),"%5.1f %5.1f %5.1f"%tuple(hkl_list),
return tuple(hkl_list)
def projection_onto_section_zero(self,hkl):
#compute this point, P, in orthonormal reciprocal coordinates
P = matrix.sqr(self.rstbx_ori.reciprocal_matrix()) * matrix.col(hkl)
# the plane normal
N = self.miller_calc.a_hat_perp
#projection of P onto the plane, given that a point on the plane O is
0,0,0 (otherwise Q=P-((P-O).dot(N))*N)
Q = P - (P.dot(N))*N
# these are still orthonormal reciprocal coordinates. Want to convert them
back to Miller indices
M = matrix.sqr(self.rstbx_ori.direct_matrix()) * Q
return M
def xy_from_00(self):
width = float(self.P.labelit_precession.pixel_width)/2.0
miller =[self.hkl_zero[self.miller_calc.idx1],
self.hkl_zero[self.miller_calc.idx2],
self.hkl_zero[self.miller_calc.idx0]]
r = self.miller_calc.reexpress.inverse()*miller
coord1 = r[0]
coord2 = r[1]
x = -width*(coord2*self.miller_calc.dmax - 1.)
y = width * (coord1*self.miller_calc.dmax + 1.)
return
x/float(self.P.labelit_precession.pixel_width),y/float(self.P.labelit_precession.pixel_width)
def get_labelit_index_parameters(self,P):
self.P = P # program parameters
H = pickle.load(open("LABELIT_possible"))
counter = [i["counter"] for i in H]
orient = [i["orient"] for i in H]
setting_index = counter.index(self.P.labelit_precession.bravais_choice)
self.rstbx_ori = orient[setting_index]
#analyze the HKL section request
self.miller_calc =
calculate_miller(self.P.labelit_precession.plot_section,self.rstbx_ori)
miller3 = self.miller_calc.section_number
hkl_list = [0.,0.,0.]
hkl_list[self.miller_calc.idx0]=miller3
self.hkl_zero = self.projection_onto_section_zero(hkl_list)
G = pickle.load(open("LABELIT_pickle"))
w = float( G['wavelength'] )
self.ai = ai_factory(G)
self.ai.setOrientation(self.rstbx_ori)
self.pixelsz = float(G["pixel_size"])
template = G["file"][G["file"].keys()[0]]
I = image_files(Spotspickle_argument_module(template))
a_file = I.filenames.FN[0]
self.template = os.path.join(a_file.cwd,a_file.template)
self.wildcard_count = self.template.count("#")
return self
def construct_image_lookup(self):
self.reverse_lookup = {}
self.all_image_centers = flex.double()
self.deltaphi = None
for x in xrange(self.P.labelit_precession.image_range[0],
self.P.labelit_precession.image_range[1]+1):
a_file = self.template.replace("#"*self.wildcard_count,
("%%0%sd"%self.wildcard_count)%x)
Q = QuickImage(a_file)
if self.deltaphi == None: self.deltaphi = Q.deltaphi
assert self.deltaphi == Q.deltaphi
phivalue = (Q.osc_start + (Q.deltaphi/2.0))%360.0
#print "image",x,"phi",phivalue
self.reverse_lookup[phivalue]=x
self.all_image_centers.append(phivalue)
self.distance_tree = AnnAdaptor(self.all_image_centers,1)
def imagefunc(self,deg):
modular_deg = deg%360.
self.distance_tree.query(flex.double([modular_deg]))
if math.sqrt(self.distance_tree.distances[0]) > self.deltaphi/2.0: return
None
image_number =
self.reverse_lookup.get(self.all_image_centers[self.distance_tree.nn[0]])
return image_number
def get_plot_values(self):
param = self.P
all_images={} # a very large memory sink. Might need a more clever
implementation later.
# relies on having all the images in memory at the same time.
init_value = 235
self.all_values =
flex.int(param.labelit_precession.pixel_width**2,init_value)
self.all_values.resize(flex.grid(param.labelit_precession.pixel_width,
param.labelit_precession.pixel_width))
scale_factor = 255./param.labelit_precession.intensity_full_scale
uc = self.rstbx_ori.unit_cell()
self.miller_calc.set_dmax(param.labelit_precession.resolution_outer)
for x in xrange(param.labelit_precession.pixel_width):
for y in xrange(param.labelit_precession.pixel_width):
#print x,y
H,K,L = self.HKL_from_xy(x,y)
resolution = uc.d((int(H),int(K),int(L)))
if x%50==0 and y%50==0:
print "%6.1f"%resolution
if resolution < param.labelit_precession.resolution_outer :
continue
for item in self.ai.predict_hkl((H,K,L),self.miller_calc.dmax):
#this is a kludge! if the HKL is recorded through the Ewald sphere
twice,
# only pick up the first time.
deg = item[2]*180./math.pi
image_number = self.imagefunc(deg)
if image_number == None: continue
pixelx = int(item[0]/self.pixelsz)
pixely = int(item[1]/self.pixelsz)
if not all_images.has_key(image_number):
filenm = self.template.replace("#"*self.wildcard_count,
("%%0%sd"%self.wildcard_count)%image_number)
Q = QuickImage(filenm)
Q.read()
all_images[image_number]=Q
this_data = all_images[image_number].linearintdata
pxl_value = this_data[pixelx,pixely]
#print "%4d %4d %7.3f%7.3f%7.3f,%7.3f,%7.2f->%3d x%4d y%4d val %6d"%(
# x,y,H,K,L,resolution,deg,image_number,pixelx,pixely,pxl_value)
plotted_value = int(pxl_value * scale_factor)
if plotted_value > 255:
plotted_value = 255
if plotted_value < 0:
plotted_value = 0
self.all_values[x,y] = 255 - plotted_value
break
from reportlab.pdfgen.canvas import Canvas
from reportlab.lib.pagesizes import letter
from reportlab.lib.units import inch,cm,mm
page_origin = (20.,190.)
boxedge = 500.
class Graph:
def __init__(self,fileout):
self.c = Canvas(fileout,pagesize=letter)
def title(self,text):
lines = text.split('\n')
self.c.setFont('Helvetica',12)
for x in xrange(len(lines)):
self.c.drawString(2*cm,(26-0.5*x)*cm,lines[x])
def image(self,pi):
self.c.drawInlineImage(pi,page_origin[0],page_origin[1],
width=boxedge,height=boxedge)
self.c.setFillColorRGB(1,0,0)
def dot(self,dc):
self.c.circle(page_origin[0]+(dc[1]*boxedge),page_origin[1]+(boxedge-dc[0]*boxedge),r=1.,stroke=0,fill=1)
def __del__(self):
self.c.save()
class GeneratePDF:
def __init__(self,project):
self.image_size = project.P.labelit_precession.pixel_width
self.R = Graph(project.P.labelit_precession.pdf_output.file)
self.data = project.all_values
self.section = project.P.labelit_precession.plot_section
self.project = project
def make_image_plots_detail(self):
import Image
mode = "L" # black and white
HWsizes = (self.image_size, self.image_size)
pil_image = Image.new(mode,HWsizes)
data = list(self.data)
pil_image.putdata(data)
uc = self.project.rstbx_ori.unit_cell()
axis = ["a*","b*","c*"][self.project.miller_calc.idx1]
axis2 = ["a*","b*","c*"][self.project.miller_calc.idx2]
angle = math.acos(self.project.miller_calc.unit_recip_basis1.dot(
self.project.miller_calc.unit_recip_basis2))*180./math.pi
self.R.title("One reciprocal lattice plane, %s. Red spot is the
%s"%(self.section,
self.section.replace("H","0").replace("K","0").replace("L","0"))+"""
Unit cell: %s"""%(uc)+"""
%s points in the --> direction starting at the red spot"""%axis+"""
%s is in the plane at a %5.1f degree angle"""%(axis2,angle) )
self.R.image(pil_image)
dotcoords = self.project.xy_from_00()
self.R.dot(dotcoords)
self.R.c.showPage()
return self
def run(args):
if args==[]:
print usage()
return
argument_interpreter = libtbx.phil.command_line.argument_interpreter(
master_phil=precession_master_phil,home_scope="hidden_symmetry")
parsed_arguments = []
for arg in args:
if arg.find("help")>=0:
param_master = get_precession_parameters(sources=[])
param_master.show(attributes_level=1)
return
else:
try:
a = argument_interpreter.process(arg=arg)
except:
print "Unknown file or keyword: %s"%arg
print
return
else:
parsed_arguments.append(a)
param_master = get_precession_parameters(sources=parsed_arguments)
param_master.show()
param = param_master.extract()
PP = precession_plot().get_labelit_index_parameters(param)
PP.construct_image_lookup()
PP.get_plot_values()
PDF = GeneratePDF(PP)
PDF.make_image_plots_detail()
if __name__=="__main__":
import sys
run(sys.argv[1:])
