From 87d0c7f38381c93d331c063210b17d24c2dbdad4 Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Tue, 22 Mar 2005 14:06:07 +0000 Subject: [PATCH] Initial revision --- ImageD11/__init__.py | 0 ImageD11/bisplev.py | 185 ++++ ImageD11/blobcorrector.py | 193 ++++ ImageD11/data.py | 70 ++ ImageD11/guiindexer.py | 160 +++ ImageD11/guimaker.py | 36 + ImageD11/guipeaksearch.py | 479 +++++++++ ImageD11/guitransformer.py | 346 +++++++ ImageD11/indexing.py | 558 +++++++++++ ImageD11/license.py | 327 +++++++ ImageD11/listdialog.py | 95 ++ ImageD11/opendata.py | 299 ++++++ ImageD11/peakmerge.py | 161 +++ ImageD11/plot3d.py | 143 +++ ImageD11/simplex.py | 283 ++++++ ImageD11/transform.py | 170 ++++ ImageD11/twodplot.py | 308 ++++++ ImageD11/unitcell.py | 333 +++++++ LICENSE | 325 +++++++ MANIFEST | 30 + PKG-INFO | 10 + README.TXT | 17 + scripts/ImageD11_gui.py | 136 +++ scripts/peaksearch.py | 115 +++ scripts/recoveromega.py | 57 ++ setup.py | 56 ++ src/bispev.f | 1880 ++++++++++++++++++++++++++++++++++++ src/closest.c | 219 +++++ src/connectedpixels.c | 455 +++++++++ src/dset.h | 98 ++ src/splines.c | 218 +++++ 31 files changed, 7762 insertions(+) create mode 100644 ImageD11/__init__.py create mode 100644 ImageD11/bisplev.py create mode 100644 ImageD11/blobcorrector.py create mode 100644 ImageD11/data.py create mode 100644 ImageD11/guiindexer.py create mode 100644 ImageD11/guimaker.py create mode 100644 ImageD11/guipeaksearch.py create mode 100644 ImageD11/guitransformer.py create mode 100644 ImageD11/indexing.py create mode 100644 ImageD11/license.py create mode 100644 ImageD11/listdialog.py create mode 100644 ImageD11/opendata.py create mode 100644 ImageD11/peakmerge.py create mode 100644 ImageD11/plot3d.py create mode 100644 ImageD11/simplex.py create mode 100644 ImageD11/transform.py create mode 100644 ImageD11/twodplot.py create mode 100644 ImageD11/unitcell.py create mode 100644 LICENSE create mode 100644 MANIFEST create mode 100644 PKG-INFO create mode 100644 README.TXT create mode 100644 scripts/ImageD11_gui.py create mode 100644 scripts/peaksearch.py create mode 100644 scripts/recoveromega.py create mode 100644 setup.py create mode 100644 src/bispev.f create mode 100644 src/closest.c create mode 100644 src/connectedpixels.c create mode 100644 src/dset.h create mode 100644 src/splines.c diff --git a/ImageD11/__init__.py b/ImageD11/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/ImageD11/bisplev.py b/ImageD11/bisplev.py new file mode 100644 index 00000000..53698b0f --- /dev/null +++ b/ImageD11/bisplev.py @@ -0,0 +1,185 @@ + + +# THIS CODE COMES FROM THE SCIPY PACKAGE AT www.scipy.org +# +# IT HAS BEEN COPIED TO SAVE YOU HAVING TO INSTALL ALL THE +# OTHER STUFF FROM THERE. +# +import _splines +import Numeric as n + +def myasarray(a): + if type(a) in [type(1.0),type(1L),type(1),type(1j)]: + return n.asarray([a]) + elif type(a) is n.ArrayType and len(a)==1: + # Takes care of mapping array(number) to array([number]) + return n.asarray([a[0]]) + else: + return n.asarray(a) + +def bisplev(x,y,tck,dx=0,dy=0): + """Evaluate a bivariate B-spline and its derivatives. + Description: + Return a rank-2 array of spline function values (or spline derivative + values) at points given by the cross-product of the rank-1 arrays x and y. + In special cases, return an array or just a float if either x or y or + both are floats. + Inputs: + x, y -- Rank-1 arrays specifying the domain over which to evaluate the + spline or its derivative. + tck -- A sequence of length 5 returned by bisplrep containing the knot + locations, the coefficients, and the degree of the spline: + [tx, ty, c, kx, ky]. + dx, dy -- The orders of the partial derivatives in x and y respectively. + Outputs: (vals, ) + vals -- The B-pline or its derivative evaluated over the set formed by + the cross-product of x and y. + """ + tx,ty,c,kx,ky=tck + if not (0<=dx1: return z + if len(z[0])>1: return z[0] + return z[0][0] + + + +_surfit_cache = {'tx': n.array([],'d'),'ty': n.array([],'d'), + 'wrk': n.array([],'d'), 'iwrk':n.array([],'i')} + +def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,s=None, + eps=1e-16,tx=None,ty=None,full_output=0,nxest=None,nyest=None,quiet=1): + """Find a bivariate B-spline representation of a surface. + + Description: + + Given a set of data points (x[i], y[i], z[i]) representing a surface + z=f(x,y), compute a B-spline representation of the surface. + + Inputs: + + x, y, z -- Rank-1 arrays of data points. + w -- Rank-1 array of weights. By default w=ones(len(x)). + xb, xe -- End points of approximation interval in x. + yb, ye -- End points of approximation interval in y. + By default xb, xe, yb, ye = x[0], x[-1], y[0], y[-1] + kx, ky -- The degrees of the spline (1 <= kx, ky <= 5). Third order + (kx=ky=3) is recommended. + task -- If task=0, find knots in x and y and coefficients for a given + smoothing factor, s. + If task=1, find knots and coefficients for another value of the + smoothing factor, s. bisplrep must have been previously called + with task=0 or task=1. + If task=-1, find coefficients for a given set of knots tx, ty. + s -- A non-negative smoothing factor. If weights correspond to the inverse + of the standard-deviation of the errors in z, then a good s-value + should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m=len(x) + eps -- A threshold for determining the effective rank of an over-determined + linear system of equations (0 < eps < 1) --- not likely to need + changing. + tx, ty -- Rank-1 arrays of the knots of the spline for task=-1 + full_output -- Non-zero to return optional outputs. + nxest, nyest -- Over-estimates of the total number of knots. If None then + nxest = max(kx+sqrt(m/2),2*kx+3), + nyest = max(ky+sqrt(m/2),2*ky+3) + quiet -- Non-zero to suppress printing of messages. + + Outputs: (tck, {fp, ier, msg}) + + tck -- A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and + coefficients (c) of the bivariate B-spline representation of the + surface along with the degree of the spline. + + fp -- The weighted sum of squared residuals of the spline approximation. + ier -- An integer flag about splrep success. Success is indicated if + ier<=0. If ier in [1,2,3] an error occurred but was not raised. + Otherwise an error is raised. + msg -- A message corresponding to the integer flag, ier. + + Remarks: + + SEE bisplev to evaluate the value of the B-spline given its tck + representation. + """ + x,y,z=map(myasarray,[x,y,z]) + x,y,z=map(n.ravel,[x,y,z]) # ensure 1-d arrays. + m=len(x) + if not (m==len(y)==len(z)): raise TypeError, 'len(x)==len(y)==len(z) must hold.' + if w is None: w=n.ones(m,'d') + else: w=myasarray(w) + if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m) + if xb is None: xb=x[0] + if xe is None: xe=x[-1] + if yb is None: yb=y[0] + if ye is None: ye=y[-1] + if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1' + if s is None: s=m-sqrt(2*m) + if tx is None and task==-1: raise TypeError, 'Knots_x must be given for task=-1' + if tx is not None: _curfit_cache['tx']=myasarray(tx) + nx=len(_surfit_cache['tx']) + if ty is None and task==-1: raise TypeError, 'Knots_y must be given for task=-1' + if ty is not None: _curfit_cache['ty']=myasarray(ty) + ny=len(_surfit_cache['ty']) + if task==-1 and nx<2*kx+2: + raise TypeError, 'There must be at least 2*kx+2 knots_x for task=-1' + if task==-1 and ny<2*ky+2: + raise TypeError, 'There must be at least 2*ky+2 knots_x for task=-1' + if not ((1<=kx<=5) and (1<=ky<=5)): + raise TypeError, \ + 'Given degree of the spline (kx,ky=%d,%d) is not supported. (1<=k<=5)'%(kx,ky) + if m<(kx+1)*(ky+1): raise TypeError, 'm>=(kx+1)(ky+1) must hold' + if nxest is None: nxest=kx+sqrt(m/2) + if nyest is None: nyest=ky+sqrt(m/2) + nxest,nyest=max(nxest,2*kx+3),max(nyest,2*ky+3) + if task>=0 and s==0: + nxest=int(kx+sqrt(3*m)) + nyest=int(ky+sqrt(3*m)) + if task==-1: + _surfit_cache['tx']=myasarray(tx) + _surfit_cache['ty']=myasarray(ty) + tx,ty=_surfit_cache['tx'],_surfit_cache['ty'] + wrk=_surfit_cache['wrk'] + iwrk=_surfit_cache['iwrk'] + u,v,km,ne=nxest-kx-1,nyest-ky-1,max(kx,ky)+1,max(nxest,nyest) + bx,by=kx*v+ky+1,ky*u+kx+1 + b1,b2=bx,bx+v-ky + if bx>by: b1,b2=by,by+u-kx + lwrk1=u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1 + lwrk2=u*v*(b2+1)+b2 + tx,ty,c,o = _fitpack._surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,task,s,eps, + tx,ty,nxest,nyest,wrk,lwrk1,lwrk2) + _curfit_cache['tx']=tx + _curfit_cache['ty']=ty + _curfit_cache['wrk']=o['wrk'] + ier,fp=o['ier'],o['fp'] + tck=[tx,ty,c,kx,ky] + if ier<=0 and not quiet: + print _iermess2[ier][0] + print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx), + len(ty),m,fp,s) + ierm=min(11,max(-3,ier)) + if ierm>0 and not full_output: + if ier in [1,2,3,4,5]: + print "Warning: "+_iermess2[ierm][0] + print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx), + len(ty),m,fp,s) + else: + try: + raise _iermess2[ierm][1],_iermess2[ierm][0] + except KeyError: + raise _iermess2['unknown'][1],_iermess2['unknown'][0] + if full_output: + try: + return tck,fp,ier,_iermess2[ierm][0] + except KeyError: + return tck,fp,ier,_iermess2['unknown'][0] + else: + return tck + diff --git a/ImageD11/blobcorrector.py b/ImageD11/blobcorrector.py new file mode 100644 index 00000000..d6acb055 --- /dev/null +++ b/ImageD11/blobcorrector.py @@ -0,0 +1,193 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +import Numeric as n +import bisplev, math + +class correctorclass: + + """ + Applies a spatial distortion to a peak position using a fit2d splinefile + """ + + def __init__(self,splinefile): + """ + Argument is the name of a fit2d spline file + """ + self.splinefile=splinefile + self.readfit2dspline(splinefile) + self.tolerance=1e-5 + + + def correct(self,x,y): + """ + Transform x,y in raw image coordinates into x,y of an + idealised image. Returns a tuple (x,y), expects a + pair of floats as arguments + """ + ynew=bisplev.bisplev(y,x,self.tck1)+y + xnew=bisplev.bisplev(y,x,self.tck2)+x + return xnew, ynew + + def distort(self,xnew,ynew): + """ + Distort a pair of points xnew, ynew to find where they + would be in a raw image + """ + yold=ynew-bisplev.bisplev(ynew,xnew,self.tck1) + xold=xnew-bisplev.bisplev(ynew,xnew,self.tck2) + # First guess, assumes distortion is constant + y=ynew-bisplev.bisplev(yold,xold,self.tck1) + x=xnew-bisplev.bisplev(yold,xold,self.tck2) + # Second guess should be better + error=math.sqrt((x-xold)*(x-xold)+(y-yold)*(y-yold)) + ntries=0 + while error > self.tolerance: + ntries=ntries+1 + xold=x + yold=y + y=ynew-bisplev.bisplev(yold,xold,self.tck1) + x=xnew-bisplev.bisplev(yold,xold,self.tck2) + error=math.sqrt((x-xold)*(x-xold)+(y-yold)*(y-yold)) + # print error,xold,x,yold,y + if ntries == 10: + raise "Error getting the inverse spline to converge" + return x,y + + def test(self,x,y): + """ + Checks that the correct and distort functions are indeed + inversely related to each other + """ + xnew,ynew = self.correct(x,y) + xold,yold = self.distort(xnew,ynew) + error = math.sqrt( (x-xold)*(x-xold) + (y-yold)*(y-yold)) + if error > self.tolerance: + print "Failed!" + raise "Problem in correctorclass" + + + + + def readfit2dfloats(self,fp,n,debug=0): + """ + Interprets a 5E14.7 formatted fortran line + """ + vals=[] + j=0 + while j < n: + i=0 + line=fp.readline() + while i < 5*14: + if(debug):print line + if(debug):print i,line[i:i+14] + vals.append(float(line[i:i+14]) ) + j=j+1 + i=i+14 + if j == n: break + i=i+1 + return vals + + + # read the fit2d array into a tck tuple + def readfit2dspline(self,name): + """ + Reads a fit2d spline file into a scipy/fitpack tuple, tck + A fairly long and dull routine... + """ + kx=3 + ky=3 + fp=open(name,"r") + line=fp.readline() # SPATIAL DISTORTION SPLINE INTERPOLATION COEFFICIENTS + if line[:7] != "SPATIAL": + raise SyntaxError, name+": file does not seem to be a fit2d spline file" + line=fp.readline() # BLANK LINE + line=fp.readline() # VALID REGION + line=fp.readline() # the actual valid region, assume xmin,ymin,xmax,ymax + vals=line.split() + self.xmin=float(vals[0]) + self.ymin=float(vals[1]) + self.xmax=float(vals[3]) + self.ymax=float(vals[3]) + line=fp.readline() # BLANK + line=fp.readline() # GRID SPACING, X-PIXEL SIZE, Y-PIXEL SIZE + line=fp.readline() + vals=line.split() + self.gridspacing=float(vals[0]) + self.xsize=float(vals[1]) + self.ysize=float(vals[2]) + line=fp.readline() # BLANK + line=fp.readline() # X-DISTORTION + line=fp.readline() # two integers nx1,ny1 + vals=line.split() + nx1=int(vals[0]) + ny1=int(vals[1]) + # Now follow fit2d formatted line 5E14.7 + tx1=n.array(self.readfit2dfloats(fp,nx1),n.Float32) + ty1=n.array(self.readfit2dfloats(fp,ny1),n.Float32) + c1 =n.array(self.readfit2dfloats(fp,(nx1-4)*(ny1-4)),n.Float32) + line=fp.readline() #BLANK + line=fp.readline() # Y-DISTORTION + line=fp.readline() # two integers nx2, ny2 + vals=line.split() + nx2=int(vals[0]) + ny2=int(vals[1]) + tx2=n.array(self.readfit2dfloats(fp,nx2),n.Float32) + ty2=n.array(self.readfit2dfloats(fp,ny2),n.Float32) + c2 =n.array(self.readfit2dfloats(fp,(nx2-4)*(ny2-4)),n.Float32) + fp.close() + self.tck1=(tx1,ty1,c1,kx,ky) + self.tck2=(tx2,ty2,c2,kx,ky) + + + +if __name__=="__main__": + import sys,time + from Numeric import * + start=time.time() + splinefile=sys.argv[1] + infile=sys.argv[2] + outfile=sys.argv[3] + out=open(outfile,"w") + npks=0 + tck1,tck2,xmin,xmax,ymin,ymax=readfit2dspline(splinefile) + print "Time to read fit2d splinefile %f/s"%(time.time()-start) + start=time.time() + for line in open(infile,"r").readlines(): + try: + vals = [float(v) for v in line.split()] + if vals[0]>3: + x=vals[2] + y=vals[3] + out.write(line[:-1]) + ynew=bisplev.bisplev(y,x,tck1)+y + xnew=bisplev.bisplev(y,x,tck2)+x + out.write(" ---> %f %f\n"%(xnew,ynew)) + npks=npks+1 + except ValueError: + pass + except: + raise + + out.close() + print "That took %f/s for %d peaks"%(time.time()-start,npks) + + diff --git a/ImageD11/data.py b/ImageD11/data.py new file mode 100644 index 00000000..04b0342f --- /dev/null +++ b/ImageD11/data.py @@ -0,0 +1,70 @@ + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +# ImageD11 data object + +class data: + """ Generic datatype for handling 2D images + Just a wrapper for Numeric arrays now, so that information + from file headers can be transported around + """ + def __init__(self,array,header=None): + """ Minimal things - dimensions and a numeric array """ + self.rows=array.shape[0] + self.columns=array.shape[1] + self.data=array + if header==None: + self.header={} + self.header["rows"]=array.shape[0] + self.header["columns"]=array.shape[1] + else: + self.header=header + def __add__(self,other): + if type(other) in [type(1),type(1.0)]: + return data(self.header,self.data+other) + if self.header["rows"]==other.header["rows"] and \ + self.header["columns"]==other.header["columns"]: + return data(self.header,self.data+other.data) + else: + raise Exception("Dimensions not compatible") + def __subtract__(self,other): + if type(other) in [type(1),type(1.0)]: + return data(self.header,self.data-other) + if self.header["rows"]==other.header["rows"] and \ + self.header["columns"]==other.header["columns"]: + return data(self.header,self.data-other.data) + else: + raise Exception("Dimensions not compatible") + def __multiply__(self,other): + if type(other) in [type(1),type(1.0)]: + return data(self.header,self.data*other) + if self.header["rows"]==other.header["rows"] and \ + self.header["columns"]==other.header["columns"]: + return data(self.header,self.data*other.data) + else: + raise Exception("Dimensions not compatible") + def __divide__(self,other): + if type(other) in [type(1),type(1.0)]: + return data(self.header,self.data/other) + if self.header["rows"]==other.header["rows"] and \ + self.header["columns"]==other.header["columns"]: + return data(self.header,self.data/other.data) + else: + raise Exception("Dimensions not compatible") + diff --git a/ImageD11/guiindexer.py b/ImageD11/guiindexer.py new file mode 100644 index 00000000..ff4c0db1 --- /dev/null +++ b/ImageD11/guiindexer.py @@ -0,0 +1,160 @@ + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +from Numeric import * +from Tkinter import * + +import time,math,unitcell,sys + +import indexing + +from listdialog import listdialog + + +class guiindexer: + + def __init__(self,parent): + """ + Parent is a hook to features of the parent gui + """ + self.parent=parent +# peaks are in self.parent.finalpeaks + self.menuitems = ( "Indexing", 0, + [ ( "Load g-vectors", 0, self.loadgv), + ( "Plot x/y/z", 5, self.plotxyz ), + ( "Load parameters", 1, self.loadfileparameters), + ( "Edit parameters", 0, self.editparameters), + ( "Assign peaks to powder rings", 0, self.assignpeaks), + ( "Generate trial orientations",0, self.find), + ( "Score trial orientations",0, self.scorethem), + ( "Save parameters", 0, self.saveparameters), + ( "Write out indexed peaks",0,self.saveindexing ) + ] ) + self.parameters={} + self.indexer=indexing.indexer(unitcell=self.parent.unitcell,gv=self.parent.gv) + self.loadparameters() + self.plot3d=None + + def loadgv(self): + if self.parent.gv!=None: + from tkMessageBox import askyesno + if not askyesno("Overwrite current g-vectors?","Loading new will overwrite the ones in memory now, are you sure?"): + return + filename=self.parent.opener.show(title="File containing g-vectors") + self.indexer.readgvfile(filename) + + def scorethem(self): + self.indexer.scorethem() + def assignpeaks(self): + self.indexer.assigntorings() + + def loadfileparameters(self): + filename=self.parent.opener.show(title="File containing indexing parameters") + self.loadparameters(filename) + + def loadparameters(self,filename=None): + if filename==None: + self.parameters['cosine_tol']=self.indexer.cosine_tol + self.parameters['hkl_tol']=self.indexer.hkl_tol + self.parameters['ring_1']=self.indexer.ring_1 + self.parameters['ring_2']=self.indexer.ring_2 + self.parameters['minpks']=self.indexer.minpks + self.parameters['uniqueness']=self.indexer.uniqueness + self.parameters['ds_tol']=self.indexer.ds_tol + self.parameters['wavelength']=self.indexer.wavelength + else: + try: + lines = open(filename,"r").readlines() + for line in lines: + name,value=line.split() + self.parameters[name]=value + self.indexer.cosine_tol=float(self.parameters['cosine_tol']) + self.indexer.hkl_tol= float(self.parameters['hkl_tol']) + self.indexer.ring_1=int(self.parameters['ring_1']) + self.indexer.ring_2=int(self.parameters['ring_2']) + self.indexer.minpks=int(self.parameters['minpks']) + self.indexer.uniqueness=float(self.parameters['uniqueness']) + self.indexer.ds_tol=float(self.parameters['ds_tol']) + self.indexer.wavelength=float(self.parameters['wavelength']) + except IOError: + from tkMessageBox import showinfo + showinfo("Sorry","Could not open/interpret your file") + + + def saveparameters(self,filename=None): + if filename==None: + filename=self.parent.saver.show(title="File to save indexing parameters") + f=open(filename,"w") + keys=self.parameters.keys() + keys.sort() + for key in keys: + try: + f.write("%s %f\n"%(key,self.parameters[key])) + except: + f.write("%s %s\n"%(key,self.parameters[key])) + f.close() + + + def editparameters(self): + self.parameters['cosine_tol']=self.indexer.cosine_tol + self.parameters['hkl_tol']=self.indexer.hkl_tol + self.parameters['ring_1']=self.indexer.ring_1 + self.parameters['ring_2']=self.indexer.ring_2 + self.parameters['minpks']=self.indexer.minpks + self.parameters['uniqueness']=self.indexer.uniqueness + self.parameters['ds_tol']=self.indexer.ds_tol + self.parameters['wavelength']=self.indexer.wavelength + d=listdialog(self.parent,items=self.parameters,title="Indexing parameters") + self.parameters=d.result + self.indexer.cosine_tol=float(self.parameters['cosine_tol']) + self.indexer.hkl_tol= float(self.parameters['hkl_tol']) + self.indexer.ring_1=int(self.parameters['ring_1']) + self.indexer.ring_2=int(self.parameters['ring_2']) + self.indexer.minpks=int(self.parameters['minpks']) + self.indexer.ds_tol=float(self.parameters['ds_tol']) + self.indexer.uniqueness=float(self.parameters['uniqueness']) + try: + self.indexer.wavelength=float(self.parameters['wavelength']) + except: + self.indexer.wavelength=None + + + def plotxyz(self): + """ + Plots the x,y arrays being used + """ + import plot3d + if self.indexer.gv!=None: + if self.plot3d==None: + self.plot3d = plot3d.plot3d(self.parent,self.indexer.gv) + self.plot3d.go() + print self.plot3d + else: + self.plot3d.changedata(self.indexer.gv) + + + def find(self): + self.indexer.find() + + + def saveindexing(self,filename=None): + if filename==None: + filename=self.parent.saver.show(title="File to save indexing parameters") + self.indexer.saveubis(filename) diff --git a/ImageD11/guimaker.py b/ImageD11/guimaker.py new file mode 100644 index 00000000..18f1abfc --- /dev/null +++ b/ImageD11/guimaker.py @@ -0,0 +1,36 @@ + + +# This is copied from a book somewhere +# either programming python or the cookbook. +# look it up and give credit please..! Probably the book from Mark Lutz?? + +import sys +from Tkinter import * + +class GuiMaker(Frame): # Inherit from Tk frame + menuBar=[] + def __init__(self,parent=None): + Frame.__init__(self,parent) + self.pack(expand=YES, fill=BOTH) + self.statusbar=Label(self,text="Welcome to ImageD11_v0.4") + self.statusbar.pack(side=BOTTOM) + self.start() + self.makeWidgets() + self.makeMenuBar() + + def makeMenuBar(self): + menubar = Menu(self.master) + self.master.config(menu=menubar) + for (name, key, items) in self.menuBar: + pulldown = Menu(menubar) + self.addMenuItems(pulldown,items) + menubar.add_cascade(label=name, underline=key, menu=pulldown) + + def addMenuItems(self, menu, items): + for item in items: + if item=="separator": + menu.add_separator({}) + else: + menu.add_command(label = item[0], underline = item[1], command=item[2] ) + + diff --git a/ImageD11/guipeaksearch.py b/ImageD11/guipeaksearch.py new file mode 100644 index 00000000..14a289ff --- /dev/null +++ b/ImageD11/guipeaksearch.py @@ -0,0 +1,479 @@ + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +from Numeric import * +import time + +TOLERANCE = 1.0 # Pixel separation for combining peaks + +class peak: + def __init__(self,line,omega,threshold,num): + self.omega=omega + self.num=num + self.line=line + self.threshold=threshold + v = [float(x) for x in line.split()] + self.np=v[0] + self.avg=v[1] + self.x=v[2] + self.y=v[3] + self.xc=v[4] + self.yc=v[5] + self.sigx=v[6] + self.sigy=v[7] + self.covxy=v[8] + self.forgotten=False + + def combine(self,other): + # + # If the thresholds are different, we see the same peaks + # twice. We will just pick the one with the higher threshold + # + # Otherwise we try to average them + # + if self.threshold > other.threshold and self.omega==other.omega: + other.forgotten=True + return self + else: + if self.threshold < other.threshold and self.omega==other.omega: + self.forgotten=True + return other + else: + np = self.np + other.np + threshold=self.threshold + s = self.avg*self.np + other.avg*other.np + avg = s/np + omega = (self.omega*self.np*self.avg + other.omega*other.np*other.avg)/s + num = (self.num *self.np*self.avg + other.num *other.np*other.avg)/s + x = (self.x *self.np*self.avg + other.x *other.np*other.avg)/s + y = (self.y *self.np*self.avg + other.y *other.np*other.avg)/s + xc = (self.xc *self.np*self.avg + other.xc *other.np*other.avg)/s + yc = (self.yc *self.np*self.avg + other.yc *other.np*other.avg)/s + sigx = (self.sigx *self.np*self.avg + other.sigx *other.np*other.avg)/s + sigy = (self.sigy *self.np*self.avg + other.sigy *other.np*other.avg)/s + covxy = (self.covxy*self.np*self.avg + other.covxy*other.np*other.avg)/s + self.forgotten=True + other.forgotten=True + # Make a new line + line = "%d %f %f %f %f %f %f %f %f"%(np,avg,x,y,xc,yc,sigx,sigy,covxy) + return peak(line,omega,threshold,num) + def __cmp__(self,other): + o=cmp(self.omega,other.omega) + if o == 0: + x=cmp(self.xc,other.xc) + if x ==0 : + return cmp(self.yc,other.yc) + else: + return x + else: + return o + + def __eq__(self,other): + try: + if abs(self.xc - other.xc) < TOLERANCE and abs(self.yc - other.yc) < TOLERANCE : + return True + else: + return False + except: + print self,other + raise + + +class pkimage: + """ + Contains the header information from an image + """ + def __init__(self,name): + self.name=name + self.header={} + self.header["File"]=name + def addtoheader(self,headerline): + if headerline.find("=")>0: + h,v=headerline[1:].split("=") + self.header[h.lstrip().rstrip()]=v + return + + def otherheaderstuff(self,headerline): + if headerline.find("Processed")>0: + return + if headerline.find("Spatial")>0: + return + if headerline.find("Threshold")>0: + return + if headerline.find("Number_of_pixels")>0: + return + else: # No equals sign means a threshold level or titles + print headerline + raise + return + + + +class guipeaksearcher: + + def __init__(self,parent): + """ + Parent is a hook to features of the parent gui + """ + self.parent=parent + self.lines=None + self.allpeaks=None + self.merged=None + self.final=None + self.menuitems = ( "PeakSearching", 4, + [ ( "Search raw images" , 0, self.searchraw ), + ( "Read pks file", 0, self.readpeaks ), + ( "Harvest peaks", 0, self.harvestpeaks), + ( "Merge peaks", 0, self.mergepeaks ), + ( "Filter peaks", 0, self.filter ), + ( "Save good peaks",2, self.savepeaks) ] ) + + pass + + def searchraw(self): + from tkMessageBox import showinfo + showinfo("Sorry","""Not implemented for gui, please use the jonpeaksearch script on amber/crunch for now""") + + def readpeaks(self,filename=None): + if filename == None: + filename = self.parent.opener.show(title="Peaksearch results file") + try: + self.lines = open(filename,"r").readlines() + except IOError: + showinfo("Sorry","Could not read your file"+filename) + return + # Get a list of filenames, omega angles + self.images=[] + i=0 + for line in self.lines: + i=i+1 +# print line[0:5] + if line[0:6]=="# File": + name = line.split()[-1] + currentimage=pkimage(name) + self.images.append(currentimage) + try: + imagenumber = int(name[-8:-4]) + except: + imagenumber = -1 + currentimage.linestart=i + currentimage.imagenumber=imagenumber + continue + if line[0]=="#": + currentimage.addtoheader(line) + continue + self.imagenumbers=zeros(len(self.images),Int) + self.omegas =zeros(len(self.images),Float) + i=0 + while i < len(self.images): + self.imagenumbers[i]=int(self.images[i].imagenumber) + self.omegas[i] =float(self.images[i].header["Omega"]) + i=i+1 + import twodplot + self.parent.twodplotter.adddata( + ("number versus omega", + twodplot.data( + self.imagenumbers,self.omegas, + {"xlabel" : "imagenumber", + "ylabel" : "Omega", + "title" : self.images[0].name+"..."+self.images[-1].name} + ) # data + )# tuple to plotter + ) # call + + def harvestpeaks(self): + # Check we have read in a file already + if self.lines==None: + # we have not read the file yet + self.readpeaks() + # Now we need to select the range of peaks to use + from tkMessageBox import askyesno + ans = askyesno("Have you selected a sensible range of images?",""" +Use the mouse to select the range of image numbers and omega angles that +you want to use from the graph on the screen. + +If all is OK, then say yes now and we will try to harvest the peaks. +Otherwise, say no, select the right range and come back "harvestpeaks" again +""" ) + print ans + # We now have the ranges in imagenumber and omega from + # the graph + numlim=self.parent.twodplotter.a.get_xlim() + omlim=self.parent.twodplotter.a.get_ylim() + start=time.time() + peaks=[] + for image in self.images: + # Check + om=float(image.header["Omega"]) + if image.imagenumber < numlim[1] and \ + image.imagenumber > numlim[0] and \ + om < omlim[1] and om > omlim[0] : + # Read peaks + i=image.linestart+1 + line=self.lines[i] + maxi=len(self.lines)-1 + npks=0 + while line.find("# File")<0 and i < maxi: + if line.find("Threshold")>0: + threshold=float(line.split()[-1]) + if line[0]!='#' and len(line)>10: + peaks.append(peak(line, om, threshold, image.imagenumber)) + npks=npks+1 + i=i+1 + line=self.lines[i] + print "Time to read into lists",time.time()-start + start=time.time() + # Sort by omega, then x, then y + peaks.sort() + print "Time to sort by omega:",time.time()-start + self.allpeaks=peaks + from tkMessageBox import showinfo + showinfo("Harvested peaks","You have a total of %d peaks,"%(len(self.allpeaks))+ + "no peaks have been merged") + + + + def mergepeaks(self): + if self.allpeaks==None: + self.harvestpeaks() + if self.allpeaks==None: + return + start=time.time() + from tkMessageBox import showinfo + # Now go through blocks in omega order, sorting by x,y + npeaks=len(self.allpeaks) + merge1=[self.allpeaks[0]] + i=1 + while i < npeaks: + # merge peaks with same omega values + if self.allpeaks[i] == merge1[-1]: + merge1[-1]=merge1[-1].combine(self.allpeaks[i]) + else: + merge1.append(self.allpeaks[i]) + i=i+1 + peaks=merge1 + npeaks=len(peaks) + print "After merging equivalent omegas",npeaks,time.time()-start + # Now merge the same peak on adjacent omega angles. + # First make a list of unique omega angles + i=0 + start=time.time() + uniq={} + while i < npeaks: + omega = peaks[i].omega + i=i+1 + if uniq.has_key(omega): + continue + else: + uniq[omega]=i + # + nomega=len(uniq.keys()) + print "Number of different omegas=",nomega,time.time()-start + # Now merge peaks with adjacent omega angles + # Need to find peaks which match each other + # Different threshold levels should already have been merged + start=time.time() + i=0 + merged=[] + keys=uniq.keys() + keys.sort() + prevframe=[] + while i < nomega-2: + first=uniq[keys[i]] + last =uniq[keys[i+1]] + # + # Active peaks are present on the current frame + # These can merge with the next frame + # + active=peaks[first:last] + check=len(active)+len(prevframe) + ncarry=len(prevframe) + active=active+prevframe + prevframe=[] + if len(active) != check: + raise "Problem here - lost something" + nextlast =uniq[keys[i+2]] + nextframe=peaks[last:nextlast] + print "Setting %-5d with %-6d peaks on this and %-5d peaks on next, %-5d in buffer\r"%(i,last-first,nextlast-last,ncarry), + for peak1 in active: + m=0 + if peak1.forgotten: + continue + for peak2 in nextframe: + if peak2.forgotten: + continue + if peak1==peak2: + newpeak=peak1.combine(peak2) + m=1 + break # Hope we would only overlap one peak + if m==0: + # This peak is no longer active, it did not overlap + merged.append(peak1) + else: + # This peak did overlap, we need to add it for scanning the next frame + prevframe.append(newpeak) + i=i+1 + # Now we finished the loop + # the last frame just needs to merge with anything in prevframe + print + print "final frame" + if nomega == 2: + active=peaks[keys[0]:keys[1]] + first=uniq[keys[1]] + last=uniq[keys[2]] + if nomega > 2: + first=uniq[keys[i]] + last =uniq[keys[i+1]] + active=prevframe + if nomega >= 2: + for peak1 in active: + m=0 + if peak1.forgotten: + continue + for peak2 in peaks[first:last]: + if peak2.forgotten: + continue + if peak1 == peak2: + newpeak=peak1.combine(peak2) + m=1 + break + if m==1: + merged.append(newpeak) + else: + merged.append(peak1) + else: + merged=peaks + # Everything ought to be covered here + merged.sort() + self.merged=merged + from tkMessageBox import showinfo + showinfo("Finished merging peaks","You have a total of "+str(len(self.merged))+" after merging") + print "That took",time.time()-start + return + + def filter(self): + if self.merged==None: + self.mergepeaks() + # Eventually we should offer filters based on shape/intensity/x/y etc + # + # For now just throw everything into an array + x=[] # To hold xc, yc, omega + y=[] + om=[] + for peak in self.merged: + x.append(peak.xc) + y.append(peak.yc) + om.append(peak.omega) + self.parent.finalpeaks=array([x,y,om]) + import twodplot + self.parent.twodplotter.hideall() + self.parent.twodplotter.adddata( + ( "Filtered peaks", + twodplot.data( + self.parent.finalpeaks[0,:], + self.parent.finalpeaks[1,:], + {"xlabel" : "x", + "ylabel" : "y", + "title" : "Peak positions on detector"} + ) ) ) + # Need to filter based on x,y + # also based on intensity + # also based on shape + + def savepeaks(self,filename=None): + # Write out minimal information + # list of xcorr,ycorr,omega + if self.parent.finalpeaks!=None: + p=self.parent.finalpeaks + if filename==None: + filename=self.parent.saver.show(title="Filtered peak positions") + f=open(filename,"w") + f.write("# xc yc omega\n") + for i in range(p.shape[1]): + f.write("%f %f %f\n"%(p[0,i],p[1,i],p[2,i])) + f.close() + + + def slow(self,filename=None): + if filename == None: + filename = parent.opener.show(title="Name of raw peaks file") + self.images=[] + self.npix=[] + self.avg=[] + self.xraw=[] + self.yraw=[] + self.xc = [] + self.yc = [] + self.sig_x=[] + self.sig_y=[] + self.cov_xy=[] + for line in open(filename,"r").xreadlines(): +# print line[0:5] + if line[0:6]=="# File": + name = line.split()[-1] + currentimage=pkimage(name) + self.images.append(currentimage) + imagenumber = int(name[-8:-4]) + continue + if line.find("Threshold")>0: + threshold = float(line.split()[-1]) + continue + if line.find(" Omega")>0: + omega = float(line.split()[-1]) + continue + if line[0]=="#": + currentimage.addtoheader(line) + continue + if len(line)>5: +# print line + values = [float(x) for x in line.split()] + self.npix.append(values[0]) + self.avg.append(values[1]) + self.xraw.append(values[2]) + self.yraw.append(values[3]) + self.xc.append(values[4]) + self.yc.append(values[5]) + self.sig_x.append(values[6]) + self.sig_y.append(values[7]) + self.cov_xy.append(values[8]) + start=time.time() + self.npix = array(self.npix) + self.avg = array(self.avg) + self.xraw = array(self.xraw) + self.yraw = array(self.yraw) + self.xc = array(self.xc) + self.yc = array(self.yc) + self.sig_x = array(self.sig_x) + self.sig_y = array(self.sig_y) + self.cov_xy= array(self.cov_xy) + print time.time()-start + + + +if __name__=="__main__": + import sys,time + testfile=sys.argv[1] + object = guipeaksearcher(None) + start=time.time() +# import profile +# profile.run('object.readpeaks(testfile)','prof') + object.readpeaks(testfile) + print "That took",time.time()-start,"/s" +# print object.peaks.shape +# print object.peaks[:10,:] diff --git a/ImageD11/guitransformer.py b/ImageD11/guitransformer.py new file mode 100644 index 00000000..80fcd6e5 --- /dev/null +++ b/ImageD11/guitransformer.py @@ -0,0 +1,346 @@ + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +from Numeric import * +from Tkinter import * + +import time,math,unitcell,sys + +import transform + +from listdialog import listdialog + + + +class guitransformer: + + def __init__(self,parent): + """ + Parent is a hook to features of the parent gui + """ + self.parent=parent +# peaks are in self.parent.finalpeaks + self.menuitems = ( "Transformation", 0, + [ ( "Load filtered peaks", 0, self.loadfiltered), + ( "Plot y/z", 5, self.plotyz ), + ( "Load parameters", 1, self.loadfileparameters), + ( "Edit parameters", 0, self.editparameters), + ( "Plot tth/eta", 0, self.plotreta ), + ( "Add unit cell peaks",0, self.addcellpeaks), + ( "Fit",0, self.fit), + ( "Save parameters", 0, self.saveparameters), + ( "Compute g-vectors", 0, self.computegv), + ( "Save g-vectors", 0, self.savegv) + ] ) + self.twotheta=None + self.parameters={} + self.loadparameters() + + def loadfiltered(self): + if self.parent.finalpeaks!=None: + from tkMessageBox import askyesno + if not askyesno("Overwrite current peaks?","Loading new peaks will overwrite the ones in memory now, are you sure?"): + return + filename=self.parent.opener.show(title="File containing filtered peaks") + f=open(filename,"r") + # 0123456789012 + line=f.readline() + if line[0:12] !="# xc yc omega"[0:12]: + print line + from tkMessageBox import showinfo + showinfo("Sorry","That does not seem to be a filter peaks file, output from the peaksearching menu option") + return + x=[] + y=[] + om=[] + for line in f.readlines(): + v=[float(z) for z in line.split()] + x.append(v[0]) + y.append(v[1]) + om.append(v[2]) + f.close() + self.parent.finalpeaks=array([x,y,om]) + + def loadfileparameters(self): + filename=self.parent.opener.show(title="File containing detector parameters") + self.loadparameters(filename) + + def loadparameters(self,filename=None): + if filename==None: + self.parameters['z-center']=1002.832 + self.parameters['y-center']=1005.768 + self.parameters['distance']=303.7 + self.parameters['z-size']=46.77648 + self.parameters['y-size']=48.08150 + self.parameters['tilt-z']=0.0 + self.parameters['tilt-y']=0.0 + self.parameters['fit_tolerance']=0.05 + self.parameters['wavelength']=0.155 + self.parameters['cell__a']=2.9508 + self.parameters['cell__b']=2.9508 + self.parameters['cell__c']=4.6855 + self.parameters['cell_alpha']=90.0 + self.parameters['cell_beta']=90.0 + self.parameters['cell_gamma']=120.0 + self.parameters['cell_lattice_[P,A,B,C,I,F]']="P" + else: + try: + lines = open(filename,"r").readlines() + for line in lines: + name,value=line.split() + self.parameters[name]=value + except IOError: + from tkMessageBox import showinfo + showinfo("Could not open your file") + + + def saveparameters(self,filename=None): + if filename==None: + filename=self.parent.saver.show(title="File to save detector parameters") + f=open(filename,"w") + keys=self.parameters.keys() + keys.sort() + for key in keys: + try: + f.write("%s %f\n"%(key,self.parameters[key])) + except: + f.write("%s %s\n"%(key,self.parameters[key])) + f.close() + + + def editparameters(self): + d=listdialog(self.parent,items=self.parameters,title="Detector parameters") + self.parameters=d.result + + def plotyz(self): + """ + Plots the x,y arrays being used + """ + import twodplot + self.parent.twodplotter.hideall() + self.parent.twodplotter.adddata( + ( "Filtered peaks", + twodplot.data( + self.parent.finalpeaks[0,:], + self.parent.finalpeaks[1,:], + {"xlabel" : "y", "ylabel" : "z", "title" : "Peak positions on detector"} ) ) ) + + + def gof(self,args): + self.tolerance=float(self.parameters['fit_tolerance']) + self.parameters['y-center']=args[0] + self.parameters['z-center']=args[1] + self.parameters['distance']=args[2] + self.parameters['wavelength']=args[3] + self.parameters['tilt-y']=args[4] + self.parameters['tilt-z']=args[5] + self.twotheta, self.eta = transform.compute_tth_eta( self.peaks_xy, + args[0], # yc is the centre in y + float(self.parameters['y-size']) , # ys is the y pixel size + args[4] , # ty is the tilt around y + args[1], # zc is the centre in z + float(self.parameters['z-size']) , # zs is the z pixel size + args[5] , # tz is the tilt around z + args[2]*1e3) # is the sample - detector distance + w=args[3] + gof=0. + npeaks=0 + for i in range(len(self.tthc)):# (twotheta_rad_cell.shape[0]): + self.tthc[i]=transform.degrees(math.asin(self.fitds[i]*w/2)*2) + diff=take(self.twotheta,self.indices[i]) - self.tthc[i] +# print "peak",i,"diff",maximum.reduce(diff),minimum.reduce(diff) + gof=gof+dot(diff,diff) + npeaks=npeaks+len(diff) + gof=gof/npeaks + return gof*1e3 + + def fit(self): + import simplex + if self.theoryds==None: + from tkMessageBox import showinfo + showinfo("Try again","You need to have added the unitcell peaks already") + return + # Assign observed peaks to rings + w=float(self.parameters['wavelength']) + self.indices=[] + self.tthc=[] + self.fitds=[] + self.tolerance=float(self.parameters['fit_tolerance']) + tthmax = self.parent.twodplotter.a.get_xlim()[1] + print "Tolerance for assigning peaks to rings",self.tolerance,"max tth",tthmax + for i in range(len(self.theoryds)): + dsc=self.theoryds[i] + tthcalc=math.asin(dsc*w/2)*360./math.pi # degrees + if tthcalc>tthmax: + break +# print tthcalc + logicals= logical_and( greater(self.twotheta, tthcalc-self.tolerance), + less(self.twotheta, tthcalc+self.tolerance) ) + + if sum(logicals)>0: +# print maximum.reduce(compress(logicals,self.twotheta)),minimum.reduce(compress(logicals,self.twotheta)) + self.tthc.append(tthcalc) + self.fitds.append(dsc) + ind=compress(logicals,range(self.twotheta.shape[0])) + self.indices.append(ind) +# print "Ring",i,tthcalc,maximum.reduce(take(self.twotheta,ind)),minimum.reduce(take(self.twotheta,ind)) +# if raw_input("OK?")[0] not in ["Y","y"]: +# return + guess=[ float(self.parameters['y-center']) , + float(self.parameters['z-center']) , + float(self.parameters['distance']) , + float(self.parameters['wavelength']) , + float(self.parameters['tilt-y']) , + float(self.parameters['tilt-z']) ] + inc=[ .1 , .1 , .1 , 0.001 , transform.radians(0.1) , transform.radians(0.1) ] + s=simplex.Simplex(self.gof,guess,inc) + newguess,error,iter=s.minimize() + self.parameters['y-center']=newguess[0] + self.parameters['z-center']=newguess[1] + self.parameters['distance']=newguess[2] + self.parameters['wavelength']=newguess[3] + self.parameters['tilt-y']=newguess[4] + self.parameters['tilt-z']=newguess[5] + inc=[ .01 , .01 , .01 , 0.0001 , transform.radians(0.01) , transform.radians(0.01) ] + guess=newguess + s=simplex.Simplex(self.gof,guess,inc) + newguess,error,iter=s.minimize() + self.parameters['y-center']=newguess[0] + self.parameters['z-center']=newguess[1] + self.parameters['distance']=newguess[2] + self.parameters['wavelength']=newguess[3] + self.parameters['tilt-y']=newguess[4] + self.parameters['tilt-z']=newguess[5] + print newguess + + def plotreta(self): + self.peaks_xy = self.parent.finalpeaks[0:2,:] + self.x=self.peaks_xy[0,:] + self.y=self.peaks_xy[1,:] + self.omega=self.parent.finalpeaks[2,:] + print self.peaks_xy.shape + try: + self.twotheta, self.eta = transform.compute_tth_eta( self.peaks_xy, + float(self.parameters['y-center']), # yc is the centre in y + float(self.parameters['y-size']) , # ys is the y pixel size + float(self.parameters['tilt-y']) , # ty is the tilt around y + float(self.parameters['z-center']), # zc is the centre in z + float(self.parameters['z-size']) , # zs is the z pixel size + float(self.parameters['tilt-z']) , # tz is the tilt around z + float(self.parameters['distance'])*1e3) # is the sample - detector distance + self.ds = 2*sin(transform.radians(self.twotheta)/2)/float(self.parameters['wavelength']) + except: + keys=self.parameters.keys() + keys.sort() + for key in keys: + print self.parameters[key],type(self.parameters[key]) + raise + import twodplot +# self.parent.twodplotter.hideall() + self.parent.twodplotter.adddata( + ( "2Theta/Eta", + twodplot.data( + self.twotheta, + self.eta, + {"xlabel":"TwoTheta / degrees", + "ylabel":"Azimuth / degrees", + "title" :"Peak positions"} + ))) + + + def addcellpeaks(self): + # + # Given unit cell, wavelength and distance, compute the radial positions + # in microns of the unit cell peaks + # + a =float(self.parameters['cell__a']) + b =float(self.parameters['cell__b']) + c =float(self.parameters['cell__c']) + alpha =float(self.parameters['cell_alpha']) + beta =float(self.parameters['cell_beta']) + gamma =float(self.parameters['cell_gamma']) + lattice=self.parameters['cell_lattice_[P,A,B,C,I,F]'] + self.unitcell=unitcell.unitcell( [a,b,c,alpha,beta,gamma] , lattice) + self.parent.unitcell=self.unitcell + if self.twotheta==None: + self.twotheta, self.eta = transform.compute_tth_eta( self.peaks_xy, + float(self.parameters['y-center']), # yc is the centre in y + float(self.parameters['y-size']) , # ys is the y pixel size + float(self.parameters['tilt-y']) , # ty is the tilt around y + float(self.parameters['z-center']), # zc is the centre in z + float(self.parameters['z-size']) , # zs is the z pixel size + float(self.parameters['tilt-z']) , # tz is the tilt around z + float(self.parameters['distance']*1e3)) # is the sample - detector distance + # Find last peak in radius + highest = maximum.reduce(self.twotheta) + wavelength=float(self.parameters['wavelength']) + ds = 2*sin(transform.radians(highest)/2.)/wavelength + self.dslimit=ds + #print "highest peak",highest,"corresponding d*",ds + self.theorypeaks=self.unitcell.gethkls(ds) + self.unitcell.makerings(ds) + self.theoryds=self.unitcell.ringds + tths = [arcsin(wavelength*dstar/2)*2 for dstar in self.unitcell.ringds] + self.theorytth=transform.degrees(array(tths)) + import twodplot + self.parent.twodplotter.adddata( + ( "HKL peaks", + twodplot.data( + self.theorytth, + zeros(self.theorytth.shape[0]), + {'color':'r', + 'pointtype':'+'} + ))) + + def computegv(self): + """ + Using self.twotheta, self.eta and omega angles, compute x,y,z of spot + in reciprocal space + """ + self.gv = transform.compute_g_vectors(self.twotheta,self.eta,self.omega,float(self.parameters['wavelength'])) +# tthnew,etanew,omeganew=transform.uncompute_g_vectors(self.gv,float(self.parameters['wavelength'])) + self.parent.gv=self.gv +# print "Testing reverse transformations" +# for i in range(50): +# print self.twotheta[i],self.eta[i],self.omega[i], +# print tthnew[i],etanew[0][i],omeganew[0][i] + return + + def savegv(self,filename=None): + """ + Save g-vectors into a file + Use crappy .ass format from previous for now (testing) + """ + if filename==None: + filename=self.parent.saver.show(title="File to save gvectors") + f=open(filename,"w") + f.write(self.unitcell.tostring()) + f.write("\n") + f.write("# wavelength = %f\n"%( float(self.parameters['wavelength']) ) ) + f.write("# ds h k l\n") + for peak in self.theorypeaks: + f.write("%10.7f %4d %4d %4d\n"%(peak[0],peak[1][0],peak[1][1],peak[1][2])) + order = argsort(self.twotheta) + f.write("# xr yr zr xc yc ds phi omega\n") + print maximum.reduce(self.omega),minimum.reduce(self.omega) + for i in order: + f.write("%f %f %f %f %f %f %f %f \n"%(self.gv[0,i],self.gv[1,i],self.gv[2,i],self.x[i],self.y[i],self.ds[i],self.eta[i],self.omega[i])) + f.close() + + diff --git a/ImageD11/indexing.py b/ImageD11/indexing.py new file mode 100644 index 00000000..cc534a57 --- /dev/null +++ b/ImageD11/indexing.py @@ -0,0 +1,558 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +from Numeric import * +import closest + + +import math,time,sys + +def mod_360(theta, target): + """ + Find multiple of 360 to add to theta to be closest to target + """ + diff=theta-target + while diff < -180: + theta=theta+360 + diff=theta-target + while diff > 180: + theta=theta-360 + diff=theta-target + return theta + + +class indexer: + """ + A class for searching for orientation matrices + """ + def __init__(self,unitcell=None,gv=None, + cosine_tol = 0.002, + minpks = 10 , + hkl_tol=0.01, + ring_1=1, + ring_2=2, + ds_tol=0.005, + wavelength=-1, + uniqueness=0.5, + max_grains=100): + """ + Unitcell would be a unitcell object for generating hkls peaks + gv would be a 3*n array of points in reciprocal space + """ + self.unitcell=unitcell + self.gv=gv + if gv !=None: + self.gvflat=reshape(fromstring(self.gv.tostring(),Float),self.gv.shape) # Makes it contiguous in memory, hkl fast index + + self.cosine_tol=cosine_tol + self.wavelength=wavelength + self.hkl_tol=hkl_tol + self.ring_1=ring_1 + self.ring_2=ring_2 + self.uniqueness=uniqueness + self.minpks=minpks + self.ds_tol=ds_tol + self.max_grains=max_grains + + self.ubis=[] + self.scores=[] + + + def assigntorings(self): + """ + Assign the g-vectors to hkl rings + """ + # rings are in self.unitcell + limit = maximum.reduce(self.ds) + print "Maximum d-spacing considered",limit + self.unitcell.makerings(limit) + dsr=self.unitcell.ringds + self.ra=zeros(self.gv.shape[0],Int)-1 + self.na=zeros(len(dsr),Int) + print "Ring assignment array shape",self.ra.shape + tol=float(self.ds_tol) + for i in range(self.ra.shape[0]): + # Assign the peak to a ring (or no ring) + ds=self.ds[i] + best=999. + for j in range(len(dsr)): + diff=abs(ds-dsr[j]) + if diff < tol: + if diff < best: + self.ra[i]=j + best=diff + # Report on assignments + ds=array(self.ds) + print "Ring ( h, k, l) Mult total indexed to_index " + for j in range(len(dsr)): + ind = compress( equal(self.ra,j), arange(self.ra.shape[0]) ) + self.na[j]=ind.shape[0] + n_indexed = sum(where( take(self.ga,ind) > -1, 1, 0)) + n_to_index = sum(where( take(self.ga,ind) == -1, 1, 0)) + diffs = abs(take(ds,ind) - dsr[j]) + h=self.unitcell.ringhkls[dsr[j]][0] + print "Ring %-3d (%3d,%3d,%3d) %3d %5d %5d %5d"%(j,h[0],h[1],h[2],len(self.unitcell.ringhkls[dsr[j]]), + self.na[j],n_indexed,n_to_index) + # We will only attempt to index g-vectors which have been assigned to hkl rings (this gives a speedup if there + # are a lot of spare peaks + ind = compress( greater(self.ra,-1) , arange(self.ra.shape[0]) ) + self.gvr = take(self.gv , ind) + print "Using only those peaks which are assigned to rings for scoring trial matrices" + print "Shape of scoring matrix",self.gvr.shape + self.gvflat=reshape(fromstring(self.gvr.tostring(),Float),self.gvr.shape) # Makes it contiguous in memory, hkl fast index + + + + def find(self): + """ + Dig out the potential hits + """ + # Which are the rings being used for indexing + hkls1 = self.unitcell.ringhkls[self.unitcell.ringds[int(self.ring_1)]] + hkls2 = self.unitcell.ringhkls[self.unitcell.ringds[int(self.ring_2)]] + print "hkls of rings being used for indexing" + print "Ring 1:",hkls1 + print "Ring 2:",hkls2 + cosangles=[] + for h1 in hkls1: + for h2 in hkls2: + ca=self.unitcell.anglehkls(h1,h2) + cosangles.append(ca[1]) + cosangles.sort() + coses=[] + while len(cosangles)>0: + a=cosangles.pop() + if abs(a-1.)<1e-5 or abs(a+1.)<1e-5: # Throw out 180 degree angles + continue + if len(coses)==0: + coses.append(a) + continue + if abs(coses[-1]-a) > 1e-5: + coses.append(a) + print "Possible angles and cosines between peaks in rings:" + for c in coses: + print math.acos(c)*180/math.pi,c + # + # Need indices of gvectors to test + iall = arange(self.gv.shape[0]) + # + # Optionally only used unindexed peaks here? Make this obligatory + i1 = compress(logical_and(equal(self.ra,self.ring_1), self.ga==-1 ) , iall).tolist() + i2 = compress(logical_and(equal(self.ra,self.ring_2), self.ga==-1 ) , iall).tolist() + print "Number of peaks in ring 1:",len(i1) + print "Number of peaks in ring 2:",len(i2) + print "Minimum number of peaks to identify a grain",self.minpks + # print self.gv.shape + ntry=0 + nhits=0 + self.hits=[] + tol=float(self.cosine_tol) + ng=0 + mp=sqrt(sum(self.gv*self.gv,1)) + # print mp.shape + ps1 = take(self.gv,i1,0) + mp1 = take(mp,i1,0) + n1 = ps1.copy() + ps2 = take(self.gv,i2,0) + mp2 = take(mp,i2,0) + n2 = ps2.copy() + # print "mp1.shape",mp1.shape + # print "n1[:,1].shape",n1[:,1].shape + for i in range(3): + n1[:,i]=n1[:,i]/mp1 + n2[:,i]=n2[:,i]/mp2 + cs = array(coses,Float) + found=0 + hits=[] + start = time.time() + onepercent=len(i1)/100 + if onepercent < 1: onepercent=1 + start=time.time() + for i in range(len(i1)): + if i%onepercent == 0: + print "Percent done %6.3f%% ... potential hits %-6d \r"%(i*100./len(i1),len(hits)), + costheta=matrixmultiply(n2,n1[i]) + best,diff = closest.closest(costheta,cs) + if diff < tol: + hits.append( [ diff, i1[i], i2[best] ]) + print "Percent done %6.3f%% ... potential hits %-6d \r"%(i*100./len(i1),len(hits)), + print + print "Number of trial orientations generated",len(hits) + print "Time taken",time.time()-start + self.hits=hits + + + def scorethem(self): + start=time.time() + ts=0 + tor=0 + ng=0 + tol=float(self.hkl_tol) + gv=self.gvflat + all=len(self.hits) + print "Scoring",all,"potential orientations" + prog=0 + ng=0 + nuniq=0 + while len(self.hits) > 0 and ng -1 or self.ga[j]>-1: # skip things which are already assigned + continue + if i==j: + continue + try: + t0=time.time() +# print "\n\n",diff,i,j,self.gv[i,:],self.gv[j,:] + self.unitcell.orient(self.ring_1, self.gv[i,:], self.ring_2, self.gv[j,:],verbose=0) + UBI=self.unitcell.UBI + t1=time.time() + # n=self.score(UBI) + # Function call overhead actually makes a big difference here + n=closest.score(UBI,gv,tol) + t2=time.time() + tor=tor+t1-t0 + ts=ts+t2-t1 +# print self.ring_1,self.ring_2,n + except: + print i,j,self.ring_1,self.ring_2 + print self.gv[i] + print self.gv[j] + raise + if n > self.minpks: +# See if we already have this grain... + ubio=self.refine(self.unitcell.UBI.copy()) # refine the orientation + ind=self.getind(ubio) # indices of peaks indexed + ga=take(self.ga,ind) # previous grain assignments + uniqueness=sum(where(ga==-1,1,0))*1.0/ga.shape[0] + if uniqueness > self.uniqueness: + put(self.ga,ind,len(self.scores)+1) + self.ubis.append(self.unitcell.UBI.copy()) + self.scores.append(n) + ng=ng+1 + else: + nuniq=nuniq+1 +# put(self.ga,ind,ng) + print + print "Number of orientations with more than",self.minpks,"peaks is",len(self.ubis) + print "Time taken",time.time()-start + if len(self.ubis)>0: + bestfitting=argmax(self.scores) + print "UBI for best fitting\n",self.refine(self.ubis[bestfitting]) + print "Indexes",self.scorelastrefined,"peaks, with =",self.fitlastrefined + print "That was the best thing I found so far" + notaccountedfor = sum(where( logical_and(self.ga==-1, self.ra!=-1),1,0)) + print "Number of peaks assigned to rings but not indexed = ",notaccountedfor + else: + print "Try again, either with larger tolerance or fewer minimum peaks" + + def saveubis(self,filename,tol=None): + """ + Save orientation matrices + """ + f=open(filename,"w") + i=0 + import math + from ImageD11 import transform + from LinearAlgebra import inverse + for ubi in self.ubis: + if tol==None: + tol=self.hkl_tol + h=matrixmultiply(ubi,transpose(self.gv)) + hint=floor(h+0.5).astype(Int) # rounds down + gint=matrixmultiply(inverse(ubi),hint) + diff=h-hint + drlv2=sum(diff*diff,0) + ind = compress( less(drlv2,tol*tol) , arange(self.gv.shape[0]) ) + mdrlv= sum(sqrt(take(drlv2,ind)))/ind.shape[0] + f.write("Grain: %d Npeaks=%d =%f\n"%(i,ind.shape[0],mdrlv)) + i=i+1 + f.write("UBI:\n"+str(ubi)+"\n") + f.write("Peak ( h k l ) drlv ") + if self.wavelength < 0: + f.write("\n") + else: + f.write(" Omega_obs Omega_calc Eta_obs Eta_calc tth_obs tth_calc\n") + tc,ec,oc = transform.uncompute_g_vectors(gint,self.wavelength) + for j in ind: + f.write("%-6d ( % 6.4f % 6.4f % 6.4f ) % 12.8f "%(j,h[0,j],h[1,j],h[2,j],sqrt(drlv2[j])) ) + if self.wavelength < 0: + f.write("\n") + else: + # # # These should be equal to + to=math.asin(self.wavelength*self.ds[j]/2)*360/math.pi # tth observed + eo=self.eta[j] + oo=self.omega[j] + tc1=tc[j] + # Choose which is closest in eta/omega, there are two choices, {eta,omega}, {-eta,omega+180} + w=argmin( [ abs(ec[0][j] - eo) , abs(ec[1][j] - eo) ] ) + ec1=ec[w][j] + oc1=oc[w][j] + # Now find best omega within 360 degree intervals + oc1=mod_360(oc1,oo) + f.write(" % 9.4f % 9.4f % 9.4f % 9.4f % 9.4f % 9.4f"% (oo,oc1, eo,ec1 ,to,tc1)) + if self.ra[j]==-1: + f.write(" *** was not assigned to ring\n") + else: + f.write("\n") + f.write("\n\n") + # peaks assigned to rings + in_rings = compress(greater(self.ra,-1),arange(self.gv.shape[0])) + f.write("\n\nAnd now listing via peaks which were assigned to rings\n") + nleft=0 + nfitted=0 + for peak in in_rings: + # Compute hkl for each grain + h=self.gv[peak,:] + f.write("\nPeak= %-5d Ring= %-5d gv=[ % -6.4f % -6.4f % -6.4f ] omega= % 9.4f eta= % 9.4f tth= % 9.4f\n"%(peak,self.ra[peak],h[0],h[1],h[2], + self.omega[peak],self.eta[peak],self.tth[peak])) + m=0 + n=0 + bestubi=999. + for ubi in self.ubis: + hi = matrixmultiply(ubi,h) + hint = floor(hi+0.5).astype(Int) + gint = matrixmultiply(inverse(ubi),hint) + diff=hi-hint + drlv2 = sum(diff*diff,0) + if drlv2 < bestubi: + bestubi=drlv2 + besthi =hi + bestm=m + if drlv2 < tol*tol: + f.write("Grain %-5d (%3d,%3d,%3d)"%(m,hint[0],hint[1],hint[2])) + f.write(" ( % -6.4f % -6.4f % -6.4f ) "%(hi[0],hi[1],hi[2])) + # hint + tt,e,o=transform.uncompute_one_g_vector(gint,self.wavelength) +# print "obs",self.omega[peak],self.eta[peak],self.tth[peak] +# print "calc",tt,e,o + w=[ abs(e[0] - self.eta[peak]) , abs(e[1] - self.eta[peak]) ] + w=argmin( w ) + et=e[w] + om=o[w] + # Now find best omega within 360 degree intervals + om=mod_360(om,self.omega[peak]) + f.write(" omega= % 9.4f eta= %9.4f tth= %9.4f\n"%(om,et,tt) ) + n=n+1 + m=m+1 + if n==0: + f.write("Peak not assigned, closest=[ % -6.4f % -6.4f % -6.4f ] for grain %d\n"%(besthi[0],besthi[1],besthi[2],bestm)) + nleft=nleft+1 + else: + nfitted=nfitted+1 + + f.write("\n\nTotal number of peaks was %d\n"%(self.gv.shape[0])) + f.write("Peaks assigned to grains %d\n"%(nfitted)) + f.write("Peaks assigned to rings but remaining unindexed %d\n"%(nleft)) + + f.write("Peaks not assigned to rings at all %d\n"%(sum(where(self.ra==-1,1,0)))) + f.close() + + + + + def getind(self,UBI,tol=None): + """ + Returns the indices of peaks in self.gv indexed by matrix UBI + """ + if tol==None: + tol=self.hkl_tol + h=matrixmultiply(UBI,transpose(self.gv)) + hint=floor(h+0.5).astype(Int) # rounds down + diff=h-hint + drlv2=sum(diff*diff,0) + drlv2=where(self.ra==-1,tol+1,drlv2) + ind = compress( less(drlv2,tol) , arange(self.gv.shape[0]) ) + return ind + + + def score(self,UBI,tol=None): + """ + Decide which are the best orientation matrices + """ +# t0=time.time() + if tol==None: + return closest.score(UBI,self.gvflat,self.hkl_tol) + else: + return closest.score(UBI,self.gvflat,tol) + t1=time.time() + h=matrixmultiply(UBI,transpose(self.gv)) + hint=floor(h+0.5).astype(Int) # rounds down + diff=h-hint + drlv2=sum(diff*diff,0) + tol = float(self.hkl_tol) + tol = tol*tol +# print "%e"%(tol) +# for i in range(10): +# print h[:,i],hint[:,i],drlv2[i] + ind = compress( less(drlv2,tol) , arange(self.gv.shape[0]) ) +# ind = compress( less(drlv2[:npks],tol) , arange(npks) ) + t2=time.time() + print n-len(ind),"Time in c",t1-t0,"Time in python",t2-t1 +# print "Grain UBI" +# print UBI +# print "Number of peaks",ind.shape[0] + return ind + + def refine(self,UBI): + """ + Refine an orientation matrix and rescore it. + + From Paciorek et al Acta A55 543 (1999) + UB = R H-1 + where: + R = sum_n r_n h_n^t + H = sum_n h_n h_n^t + r = g-vectors + h = hkl indices + """ +# print "Orientation and unit cell refinement of" +# print "UBI\n",UBI +# print "Scores before",self.score(UBI) + # Need to find hkl indices for all of the peaks which are indexed + h=matrixmultiply(UBI,transpose(self.gv)) + hint=floor(h+0.5).astype(Int) # rounds down + diff=h-hint + drlv2=sum(diff*diff,0) + tol = float(self.hkl_tol) + tol = tol*tol + # Only use peaks which are assigned to rings for refinement + ind = compress( logical_and(less(drlv2,tol),greater(self.ra,-1)) , arange(self.gv.shape[0]) ) + scoreb4=ind.shape[0] + contribs = take(drlv2,ind) + try: + fitb4=sum(contribs)/contribs.shape[0] + except: + print "No contributing reflections for\n",UBI + raise + drlv2_old=drlv2 + R=zeros((3,3),Float) + H=zeros((3,3),Float) + for i in ind: + r = self.gv[i,:] + k = hint[:,i].astype(Float) +# print r,k + R = R + outerproduct(r,k) + H = H + outerproduct(k,k) + from LinearAlgebra import inverse + try: + HI=inverse(H) + UBoptimal=matrixmultiply(R,HI) + UBIo=inverse(UBoptimal) + except: + # A singular matrix - this sucks. + UBIo=UBI + h=matrixmultiply(UBIo,transpose(self.gv)) + hint=floor(h+0.5).astype(Int) # rounds down + diff=h-hint + drlv2=sum(diff*diff,0) + ind = compress( logical_and(less(drlv2,tol),greater(self.ra,-1)), arange(self.gv.shape[0]) ) + self.scorelastrefined=ind.shape[0] + contribs = take(drlv2,ind) + try: + self.fitlastrefined=math.sqrt(sum(contribs)/contribs.shape[0]) + except: + print "No contributing reflections for\n",UBI + print "After refinement, it was OK before ???" + raise +# for i in ind: +# print "( %-6.4f %-6.4f %-6.4f ) %12.8f %12.8f"%(h[0,i],h[1,i],h[2,i],sqrt(drlv2[i]),sqrt(drlv2_old[i])) +# print UBIo +# print "Scores after", self.score(UBIo,self.hkl_tol) +# print "diff\n",UBI-UBIo +# print "Mean drlv now",sum(sqrt(drlv2))/drlv2.shape[0], +# print "Mean drlv old",sum(sqrt(drlv2_old))/drlv2_old.shape[0] + return UBIo + + + def coverage(self): + """ + Compute the expected coverage of reciprocal space + use the min/max obs values of xp/yp/omega to work out what was measured in the scan? + No lambda or + """ + pass + + + def readgvfile(self,filename): + f=open(filename,"r") + import unitcell,math + # Lattice!!! + self.unitcell = unitcell.cellfromstring(f.readline()) + while 1: + line=f.readline() + if line[0]=="#": + if line.find("wavelength")>-1: + self.wavelength = float(line.split()[-1]) + print "Got wavelength from gv file of ",self.wavelength + continue + if line.find("ds h k l")>-1: + continue # reads up to comment line + if line.find("omega")>-1: + break + self.eta=[] # Raw peak information + self.omega=[] + self.ds=[] + self.xr=[] + self.yr=[] + self.zr=[] + self.xp=[] + self.yp=[] + for line in f.xreadlines(): + try: + v=[float(x) for x in line.split()] + self.xr.append(v[0]) + self.yr.append(v[1]) + self.zr.append(v[2]) + self.xp.append(v[3]) + self.yp.append(v[4]) + self.ds.append(v[5]) + self.eta.append(v[6]) + self.omega.append(v[7]) + except: + print line + raise +# raise "Problem interpreting the last thing I printed" + f.close() + if self.wavelength > 0: + self.tth=arcsin(array(self.ds)*self.wavelength/2)*360/math.pi + else: + self.tth=zeros(len(self.ds)) + self.gv=transpose(array( [ self.xr , self.yr, self.zr ] ,Float)) + self.ga=zeros(len(self.ds),Int)-1 # Grain assignments + + self.gvflat=reshape(fromstring(self.gv.tostring(),Float),self.gv.shape) # Makes it contiguous in memory, hkl fast index + print "Read your gv file containing",self.gv.shape +# if self.wavelength>0: +# print "First ten peaks" +# from ImageD11 import transform +# import math +# if self.gv.shape[0]>10: +# all=10 +# else: +# all=self.gv.shape[0] +# to,eo,oo=transform.uncompute_g_vectors(transpose(self.gv)[:all,:],self.wavelength) +# print "Peak tth_calc tth_gv_obs eta_obs eta_gv_obs_1 eta_gv_obs_1 omega_obs omega_gv_obs omega_gv_obs" +# for i in range(all): +# tc=math.asin(self.ds[i]*self.wavelength/2)*360/math.pi +# print i,tc,to[i],self.eta[i],eo[0][i],eo[1][i],self.omega[i],oo[0][i],oo[1][i] +# print "...." diff --git a/ImageD11/license.py b/ImageD11/license.py new file mode 100644 index 00000000..874ea920 --- /dev/null +++ b/ImageD11/license.py @@ -0,0 +1,327 @@ +license = """ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your + freedom to share and change it. By contrast, the GNU General Public + License is intended to guarantee your freedom to share and change free + software--to make sure the software is free for all its users. This + General Public License applies to most of the Free Software + Foundation's software and to any other program whose authors commit to + using it. (Some other Free Software Foundation software is covered by + the GNU Library General Public License instead.) You can apply it to + your programs, too. + + When we speak of free software, we are referring to freedom, not + price. Our General Public Licenses are designed to make sure that you + have the freedom to distribute copies of free software (and charge for + this service if you wish), that you receive source code or can get it + if you want it, that you can change the software or use pieces of it + in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid + anyone to deny you these rights or to ask you to surrender the rights. + These restrictions translate to certain responsibilities for you if you + distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether + gratis or for a fee, you must give the recipients all the rights that + you have. You must make sure that they, too, receive or can get the + source code. And you must show them these terms so they know their + rights. + + We protect your rights with two steps: (1) copyright the software, and + (2) offer you this license which gives you legal permission to copy, + distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain + that everyone understands that there is no warranty for this free + software. If the software is modified by someone else and passed on, we + want its recipients to know that what they have is not the original, so + that any problems introduced by others will not reflect on the original + authors' reputations. + + Finally, any free program is threatened constantly by software + patents. We wish to avoid the danger that redistributors of a free + program will individually obtain patent licenses, in effect making the + program proprietary. To prevent this, we have made it clear that any + patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and + modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains + a notice placed by the copyright holder saying it may be distributed + under the terms of this General Public License. The "Program", below, + refers to any such program or work, and a "work based on the Program" + means either the Program or any derivative work under copyright law: + that is to say, a work containing the Program or a portion of it, + either verbatim or with modifications and/or translated into another + language. (Hereinafter, translation is included without limitation in + the term "modification".) Each licensee is addressed as "you". + + Activities other than copying, distribution and modification are not + covered by this License; they are outside its scope. The act of + running the Program is not restricted, and the output from the Program + is covered only if its contents constitute a work based on the + Program (independent of having been made by running the Program). + Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's + source code as you receive it, in any medium, provided that you + conspicuously and appropriately publish on each copy an appropriate + copyright notice and disclaimer of warranty; keep intact all the + notices that refer to this License and to the absence of any warranty; + and give any other recipients of the Program a copy of this License + along with the Program. + + You may charge a fee for the physical act of transferring a copy, and + you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion + of it, thus forming a work based on the Program, and copy and + distribute such modifications or work under the terms of Section 1 + above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + + These requirements apply to the modified work as a whole. If + identifiable sections of that work are not derived from the Program, + and can be reasonably considered independent and separate works in + themselves, then this License, and its terms, do not apply to those + sections when you distribute them as separate works. But when you + distribute the same sections as part of a whole which is a work based + on the Program, the distribution of the whole must be on the terms of + this License, whose permissions for other licensees extend to the + entire whole, and thus to each and every part regardless of who wrote it. + + Thus, it is not the intent of this section to claim rights or contest + your rights to work written entirely by you; rather, the intent is to + exercise the right to control the distribution of derivative or + collective works based on the Program. + + In addition, mere aggregation of another work not based on the Program + with the Program (or with a work based on the Program) on a volume of + a storage or distribution medium does not bring the other work under + the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, + under Section 2) in object code or executable form under the terms of + Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + + The source code for a work means the preferred form of the work for + making modifications to it. For an executable work, complete source + code means all the source code for all modules it contains, plus any + associated interface definition files, plus the scripts used to + control compilation and installation of the executable. However, as a + special exception, the source code distributed need not include + anything that is normally distributed (in either source or binary + form) with the major components (compiler, kernel, and so on) of the + operating system on which the executable runs, unless that component + itself accompanies the executable. + + If distribution of executable or object code is made by offering + access to copy from a designated place, then offering equivalent + access to copy the source code from the same place counts as + distribution of the source code, even though third parties are not + compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program + except as expressly provided under this License. Any attempt + otherwise to copy, modify, sublicense or distribute the Program is + void, and will automatically terminate your rights under this License. + However, parties who have received copies, or rights, from you under + this License will not have their licenses terminated so long as such + parties remain in full compliance. + + 5. You are not required to accept this License, since you have not + signed it. However, nothing else grants you permission to modify or + distribute the Program or its derivative works. These actions are + prohibited by law if you do not accept this License. Therefore, by + modifying or distributing the Program (or any work based on the + Program), you indicate your acceptance of this License to do so, and + all its terms and conditions for copying, distributing or modifying + the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the + Program), the recipient automatically receives a license from the + original licensor to copy, distribute or modify the Program subject to + these terms and conditions. You may not impose any further + restrictions on the recipients' exercise of the rights granted herein. + You are not responsible for enforcing compliance by third parties to + this License. + + 7. If, as a consequence of a court judgment or allegation of patent + infringement or for any other reason (not limited to patent issues), + conditions are imposed on you (whether by court order, agreement or + otherwise) that contradict the conditions of this License, they do not + excuse you from the conditions of this License. If you cannot + distribute so as to satisfy simultaneously your obligations under this + License and any other pertinent obligations, then as a consequence you + may not distribute the Program at all. For example, if a patent + license would not permit royalty-free redistribution of the Program by + all those who receive copies directly or indirectly through you, then + the only way you could satisfy both it and this License would be to + refrain entirely from distribution of the Program. + + If any portion of this section is held invalid or unenforceable under + any particular circumstance, the balance of the section is intended to + apply and the section as a whole is intended to apply in other + circumstances. + + It is not the purpose of this section to induce you to infringe any + patents or other property right claims or to contest validity of any + such claims; this section has the sole purpose of protecting the + integrity of the free software distribution system, which is + implemented by public license practices. Many people have made + generous contributions to the wide range of software distributed + through that system in reliance on consistent application of that + system; it is up to the author/donor to decide if he or she is willing + to distribute software through any other system and a licensee cannot + impose that choice. + + This section is intended to make thoroughly clear what is believed to + be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in + certain countries either by patents or by copyrighted interfaces, the + original copyright holder who places the Program under this License + may add an explicit geographical distribution limitation excluding + those countries, so that distribution is permitted only in or among + countries not thus excluded. In such case, this License incorporates + the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions + of the General Public License from time to time. Such new versions will + be similar in spirit to the present version, but may differ in detail to + address new problems or concerns. + + Each version is given a distinguishing version number. If the Program + specifies a version number of this License which applies to it and "any + later version", you have the option of following the terms and conditions + either of that version or of any later version published by the Free + Software Foundation. If the Program does not specify a version number of + this License, you may choose any version ever published by the Free Software + Foundation. + + 10. If you wish to incorporate parts of the Program into other free + programs whose distribution conditions are different, write to the author + to ask for permission. For software which is copyrighted by the Free + Software Foundation, write to the Free Software Foundation; we sometimes + make exceptions for this. Our decision will be guided by the two goals + of preserving the free status of all derivatives of our free software and + of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY + FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN + OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES + PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED + OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS + TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE + PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, + REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING + WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR + REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, + INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING + OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED + TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY + YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER + PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE + POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest + possible use to the public, the best way to achieve this is to make it + free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest + to attach them to the start of each source file to most effectively + convey the exclusion of warranty; and each file should have at least + the "copyright" line and a pointer to where the full notice is found. + + ImageD11_v0.4 Software for beamline ID11 + Copyright (C) 2005 Jon Wright + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + Also add information on how to contact you by electronic and paper mail. + + If the program is interactive, make it output a short notice like this + when it starts in an interactive mode: + + ImageD11 version 0.4, Copyright (C) 2005 Jon Wright + ImageD11 comes with ABSOLUTELY NO WARRANTY; for details select help, license. + This is free software, and you are welcome to redistribute it + under certain conditions; + + This General Public License does not permit incorporating your program into + proprietary programs. If your program is a subroutine library, you may + consider it more useful to permit linking proprietary applications with the + library. If this is what you want to do, use the GNU Library General + Public License instead of this License. + """ diff --git a/ImageD11/listdialog.py b/ImageD11/listdialog.py new file mode 100644 index 00000000..cbff96a9 --- /dev/null +++ b/ImageD11/listdialog.py @@ -0,0 +1,95 @@ + + +# +# Look up where this comes from. +# It was copied from somewhere and modified slightly + + + +from Numeric import * +from Tkinter import * + + + +class listdialog(Toplevel): + """ + Dialog box for setting detector parameters + Takes a list of strings and numbers + """ + def __init__(self, parent, title = None, items=None): + Toplevel.__init__(self, parent) + self.transient(parent) + if title: + self.title(title) + self.parent = parent + self.result = items + body = Frame(self) + self.initial_focus = self.body(body,items) + body.pack(padx=5, pady=5) + self.buttonbox() + self.grab_set() + if not self.initial_focus: + self.initial_focus = self + self.protocol("WM_DELETE_WINDOW", self.cancel) + self.geometry("+%d+%d" % (parent.winfo_rootx()+50, + parent.winfo_rooty()+50)) + self.initial_focus.focus_set() + self.wait_window(self) + def body(self, master, items): + # create dialog body. return widget that should have + # initial focus. this method should be overridden + self.e=[] + if items!=None: + i=0 + keys=items.keys() + keys.sort() + self.keys=keys + for key in keys: + Label(master,text=key).grid(row=i) + el=Entry(master) + el.insert(END,items[key]) + el.grid(row=i,column=1) + self.e.append(el) + i=i+1 + return self.e[0] + + def buttonbox(self): + # add standard button box. override if you don't want the + # standard buttons + box = Frame(self) + w = Button(box, text="OK", width=10, command=self.ok, default=ACTIVE) + w.pack(side=LEFT, padx=5, pady=5) + w = Button(box, text="Cancel", width=10, command=self.cancel) + w.pack(side=LEFT, padx=5, pady=5) + self.bind("<Return>", self.ok) + self.bind("<Escape>", self.cancel) + box.pack() + # + # standard button semantics + def ok(self, event=None): + if not self.validate(): + self.initial_focus.focus_set() # put focus back + return + self.withdraw() + self.update_idletasks() + self.apply() + self.cancel() + + def cancel(self, event=None): + # put focus back to the parent window + self.parent.focus_set() + self.destroy() + # + # command hooks + def validate(self): + return 1 # override + + def apply(self): + retdict={} + i=0 + for item in self.e: + retdict[self.keys[i]]=item.get() + i=i+1 + self.result=retdict + print self.result + diff --git a/ImageD11/opendata.py b/ImageD11/opendata.py new file mode 100644 index 00000000..20f1c937 --- /dev/null +++ b/ImageD11/opendata.py @@ -0,0 +1,299 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +# ImageD11 generic file opener +# +# +# Returns a data object +# +from data import data +from Numeric import * + +def makename(stem,num,extn,width=4): + if width>0: + fs="%s%0"+"%d"%(width)+"d%s" + else: + fs="%s%d%s" + return fs%(stem,num,extn) + + +def opendata(filename): + if filename[-3:]=="edf": + return openedf(filename) + if filename[-3:] in ["xye","epf","inp"]: + return openxye(filename) + if filename[-5]=="." and filename[-4:].isdigit(): + return openbruker(filename) + +def openxye(filename): + h={} + h['columns']=3 + for line in open(filename,"r"): + xye.append(map(float,line.split())) + a=array(xye) + h['rows']=a.shape[1] + h['xye']=True + return data(h,a) + +def readbytestream(file,offset,x,y,nbytespp,datatype='int',signed='n', + swap='n',typeout=UInt16): + """ + Reads in a bytestream from a file (which may be a string indicating + a filename, or an already opened file (should be "rb")) + offset is the position (in bytes) where the pixel data start + nbytespp = number of bytes per pixel + type can be int or float (4 bytes pp) or double (8 bytes pp) + signed: normally signed data 'y', but 'n' to try to get back the right numbers + when unsigned data are converted to signed (python has no unsigned numeric types.) + swap, normally do not bother, but 'y' to swap bytes + typeout is the Numeric type to output, normally UInt16, but more if overflows occurred + x and y are the pixel dimensions + """ + tin="dunno" + len=nbytespp*x*y # bytes per pixel times number of pixels + if datatype=='int' and signed=='n': + if nbytespp==1 : tin=UInt8 + if nbytespp==2 : tin=UInt16 + if nbytespp==4 : tin=UInt32 + if datatype=='int' and signed=='y': + if nbytespp==1 : tin=Int8 + if nbytespp==2 : tin=Int16 + if nbytespp==4 : tin=Int32 + if datatype=='float': + tin=Float32 + if datatype=='double' : + tin=Float64 + if tin=="dunno" : + raise SyntaxError, "Did not understand what type to try to read" + opened=0 + if type(tin) == type(file): # Did we get a string or a file pointer? + f=open(file,'rb') + opened=1 + else: + f=file + f.seek(offset) + if swap=='y': + ar=array(reshape(byteswapped(fromstring(f.read(len),tin)),(x,y)),typeout) + else: + ar=array(reshape(fromstring(f.read(len),tin),(x,y)),typeout) + if(opened):f.close() + return(ar) + + +def readbrukerheader(file): + """ + Reads a Bruker file header into a Python dictionary + file=filename or file pointer + """ + s="string" # s is a string + if type(s) == type(file): # if arg is a string, open file, else treat as file object + f=open(file,"rb") + opened=1 + else: + f=file + opened=0 # opened var flags to close again if we open the file + i=80 + hs=f.read(512) # always start with a 512 byte header + block=hs + Header={} # dict to take results + while i < 512 : # wander along the 512 bytes + key,val=block[i-80:i].split(":",1) # uses 80 char lines in key : value format + key=key.strip() # remove the whitespace (why?) + val=val.strip() + if Header.has_key(key): # append lines if key already there + Header[key]=Header[key]+'\n'+val + else: + Header[key]=val + i=i+80 # next 80 characters + nhdrblks=int(Header['HDRBLKS']) # we must have read this in the first 512 bytes. + # print "got first chunk, headerblocks is",nhdrblks + # Now read in the rest of the header blocks, appending to what we have + rest=f.read(512*(nhdrblks-1)) + block = block[i-80:512] + rest + hs=hs+rest + j=512*nhdrblks + while i < j : + # print i,"*",block[i-80:i].strip(),"*" + if block[i-80:i].find(":") > 0: # as for first 512 bytes of header + key,val=block[i-80:i].split(":",1) + key=key.strip() + val=val.strip() + if Header.has_key(key): + Header[key]=Header[key]+'\n'+val + else: + Header[key]=val + i=i+80 + Header['datastart']=f.tell() # make a header item called "datastart" + if(opened):f.close() +# print hs + Header['headerstring']=hs + print Header['datastart'],len(hs) + return Header # s + + + +def readbruker(file): + """ + Reads in a Bruker file, returning the data and header + file may be a string or file object + FIXME we should later modify to take ROI ranges somehow (xmin,xmax,ymin,ymax) + """ + s="string" # s is a string + if type(s) == type(file): # if arg is a string, open file, else treat as file object + f=open(file,"rb") + opened=1 + else: + f=file + opened=0 + Header=readbrukerheader(f) # pass in the file pointer, it stays open + npixelb=int(Header['NPIXELB']) # you had to read the Bruker docs to know this! + rows =int(Header['NROWS']) + cols =int(Header['NCOLS']) + # We are now at the start of the image - assuming readbrukerheader worked + size=rows*cols*npixelb + data=readbytestream(f,f.tell(),rows,cols,npixelb,datatype="int",signed='n',swap='n') + no=int(Header['NOVERFL']) # now process the overflows + if no>0: # Read in the overflows + # need at least Int32 sized data I guess - can reach 2^21 + data=data.astype(UInt32) + # 16 character overflows, 9 characters of intensity, 7 character position + for i in range(no): + ov=f.read(16) + intensity=int(ov[0:9]) + position=int(ov[9:16]) + r=position%rows # relies on python style modulo being always + + c=position/rows # relies on truncation down + #print "Overflow ",r,c,intensity,position,data[r,c],data[c,r] + data[c,r]=intensity + f.close() + Header["rows"]=rows + Header["columns"]=cols + return Header,data + + +def openbruker(filename): + return data(readbruker(filename)) + +def edfheader(filename): + f=open(filename,"rb") + # Header is 1024 byte blocks, terminated by "}" + fh=f.read(1024) + i=1023 + while fh.find("}\n")<0: + fh+=f.read(1024) + # Interpret header + headeritems=fh[1:-1].split(";") + hd={} + for item in headeritems: + if item.find("=")>0: + hd[item.split("=")[0].lstrip().rstrip()]=item.split("=")[1].lstrip().rstrip() + hd["rows"]=int(hd["Dim_1"]) + hd["columns"]=int(hd["Dim_2"]) +# print hd['description'] +# line=hd['description'].split(" ") +# print line +# while 1: +# try: +# value=line.pop(-1) +# name=line.pop(-1) +# print "n,v",name,value +# hd[name]=value +# except: +# break + #print hd + # seek back from the end of the file + f.seek( -int(hd["Size"]) , 2 ) + datastring=f.read( int(hd["Size"]) ) + f.close() + return hd + +def openedf(filename): + f=open(filename,"rb") + # Header is 1024 byte blocks, terminated by "}" + fh=f.read(1024) + i=1023 + while fh.find("}\n")<0: + fh+=f.read(1024) + # Interpret header + headeritems=fh[1:-1].split(";") + hd={} + for item in headeritems: + if item.find("=")>0: + hd[item.split("=")[0].lstrip().rstrip()]=item.split("=")[1].lstrip().rstrip() + hd["rows"]=int(hd["Dim_1"]) + hd["columns"]=int(hd["Dim_2"]) +# print hd['description'] +# line=hd['description'].split(" ") +# print line +# while 1: +# try: +# value=line.pop(-1) +# name=line.pop(-1) +# print "n,v",name,value +# hd[name]=value +# except: +# break + #print hd + # seek back from the end of the file + f.seek( -int(hd["Size"]) , 2 ) + datastring=f.read( int(hd["Size"]) ) + + f.close() + # Convert datastring to Numeric array + if hd["DataType"]=="UnsignedShort": numerictype=UInt16 + else: + raise TypeError("Unimplemented edf filetype") + if hd["ByteOrder"]=="LowByteFirst": + ar=reshape( + fromstring(datastring,numerictype), + (hd["rows"],hd["columns"]) ) + else: + ar=reshape( + fromstring(datastring,numerictype).byteswapped(), + (hd["rows"],hd["columns"]) ) + return data(ar,hd) + + +if __name__=="__main__": + import sys,time + if len(sys.argv)!=2: + print "Usage: %s filename" % (sys.argv[0]) + sys.exit() + t1=time.time() + testdata=opendata(sys.argv[1]) + t2=time.time() + print "Time to read file =",t2-t2,"/s" + print "Rows =",testdata.header['rows'] + print "Columns =",testdata.header['columns'] + print "Maximum =",maximum.reduce(ravel(testdata.data)) + print "Minimum =",minimum.reduce(ravel(testdata.data)) + t3=time.time() + print "Time native ops =",t3-t2,"/s" + s =sum( ravel(testdata.data).astype(Float32) ) # 16 bit overflows + sq=sum( pow( ravel(testdata.data).astype(Float32), 2) ) + n=testdata.header['rows']*testdata.header['columns'] + print "Sum =",s + print "Average =",s/n + print "Variance =",sqrt(sq/n - (s/n)*(s/n)) + t4=time.time() + print "Time float ops =",t4-t3,"/s" + print "Total run time =",t4-t1,"/s" + diff --git a/ImageD11/peakmerge.py b/ImageD11/peakmerge.py new file mode 100644 index 00000000..e223e74c --- /dev/null +++ b/ImageD11/peakmerge.py @@ -0,0 +1,161 @@ + + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +import sys + + +TOLERANCE = 1.0 + +class peak: + def __init__(self,xc,yc,omega,avg,np): + self.xc=xc + self.yc=yc + self.avg=avg + self.np=np + self.omega=omega + + self.sorter=omega + self.forget=False + + + def combine(self,other): + np = self.np+other.np + s = (self.avg*self.np + other.avg*other.np) + avg = s/np + omega = (self.omega*self.np*self.avg + other.omega*other.np*other.avg)/s + xc = (self.xc*self.np*self.avg + other.xc*other.np*other.avg)/s + yc = (self.yc*self.np*self.avg + other.yc*other.np*other.avg)/s + newpeak = (xc,yc,omega,avg,np) + # print "Returning:",newpeak + self.forget=True + other.forget=True + return peak(xc,yc,omega,avg,np) + + def __eq__(self,other): + try: + if abs(self.xc - other.xc) < TOLERANCE and abs(self.yc - other.yc) < TOLERANCE : + return True + else: + return False + except: + print self,other + raise + + def __str__(self): + return "xc= %f yc= %f Omega= %f Avg= %f np= %d"%(self.xc,self.yc,self.omega,self.avg,self.np) + def __repr__(self): + return "xc= %f yc= %f Omega= %f Avg= %f np= %d"%(self.xc,self.yc,self.omega,self.avg,self.np) + def __cmp__(self,other): + o=cmp(self.omega,other.omega) + if o == 0: + x=cmp(self.xc,other.xc) + if x ==0 : + return cmp(self.yc,other.yc) + else: + return x + else: + return o + + +peakdict = {} + +try: + infile = sys.argv[1] + outfile = sys.argv[2] +except: + print "%s infile outfile minimum_threshold"%(sys.argv[0]) +try: + th=float(sys.argv[3]) +except: + th=0. + print "Using zero minimum threshold" + +Omega = -9999. + +lines = open(infile,"r").xreadlines() +i=0 +t=0 +print "# Opened file %s and read it"%(infile) +while i < len(lines): + if lines[i].find(" Omega")>0: + Omega = lines[i].split()[-1] + if lines[i].find("Threshold level")>0: +# print lines[i],lines[i].split() + t = float(lines[i].split()[-1]) + if th > t: + i=i+1 + continue + if lines[i][0] != "#" and len(lines[i])>10: + vals = [ float(x) for x in lines[i].split() ] + p=peak(vals[4],vals[5],float(Omega),vals[1],vals[0]) + try: + found=0 + for index in range(len(peakdict[Omega])): + if p == peakdict[Omega][index]: + peakdict[Omega][index]=p.combine(peakdict[Omega][index]) + found=1 + if found==0: + peakdict[Omega].append(p) + except KeyError: + peakdict[Omega]=[p] + i=i+1 # while loop + +print "# Made a dictionary of peaks" +keys = peakdict.keys() +keys.sort(lambda x,y : cmp(float(x),float(y))) +print "# Sorted peaks " +peaklist=[] + + +#print peakdict["81"] + +i=0 +while i < len(keys)-1: + # Check if this peak is on the next frame + for k in range(len(peakdict[keys[i]])): # peaks on this frame + for j in range(len(peakdict[keys[i+1]])): + if peakdict[keys[i]][k]==peakdict[keys[i+1]][j]: # peaks on next frame + peakdict[keys[i+1]][j]=peakdict[keys[i]][k].combine(peakdict[keys[i+1]][j]) + peakdict[keys[i]][k].forget=True + i=i+1 + +print "# Merged peak on adjacent frames" + +peaklist=[] +i=0 +while i < len(keys)-1: + # Check if this peak is on the next frame + for mypeak in peakdict[keys[i]]: # peaks on this frame + if not mypeak.forget: + peaklist.append(mypeak) +# print mypeak + i=i+1 + +print "# Made final peak list" +print "# Number of merged peaks",len(peaklist) +peaklist.sort() +outfile=open(sys.argv[2],"w") +for p in peaklist: + if p.avg>1000: +# print p + outfile.write("%s\n"%(p)) + +print "# All done" diff --git a/ImageD11/plot3d.py b/ImageD11/plot3d.py new file mode 100644 index 00000000..34fc6b5b --- /dev/null +++ b/ImageD11/plot3d.py @@ -0,0 +1,143 @@ +#! + +# This is statement is required by the build system to query build info +if __name__ == '__build__': + raise Exception + + +import string +__version__ = string.split('$Revision$')[1] +__date__ = string.join(string.split('$Date$')[1:3], ' ') +__author__ = 'Jon Wright from example by Tarn Weisner Burton ' + +try: + from Numeric import * +except: + import sys + print "This demo requires the Numeric extension, sorry." + sys.exit() + + +from OpenGL.GL import * +from OpenGL.Tk import * + + +import sys + + +class plot3d(Toplevel): + def __init__(self,parent,data=None): + Toplevel.__init__(self,parent) + if data!=None: + xyz=data.copy() + else: + xyz=array([0,0,0]) + self.ps=StringVar() + self.ps.set('1.') + self.pointsize=1. + self.npeaks=xyz.shape[0] + + self.o = Opengl(self, width = 400, height = 400, double = 1) + self.o.redraw = self.redraw + self.o.autospin_allowed = 1 + self.o.fovy=5 + self.o.near=1e6 + self.o.far=1e-6 + import math + self.o.distance=maximum.reduce(ravel(xyz))*4/math.tan(self.o.fovy*math.pi/180) + glVertexPointerd(xyz) + glEnableClientState(GL_VERTEX_ARRAY) + f=Frame(self) + Button(f,text="Help",command=self.o.help).pack(side=LEFT) + Button(f,text="Reset",command=self.o.reset).pack(side=LEFT) + Button(f,text="Pointsize",command=self.setps).pack(side=LEFT) + Entry(f,textvariable=self.ps).pack(side=LEFT) + Button(f,text="Quit",command=self.goaway).pack(side=RIGHT) + self.dataoff=0 + self.o.pack(side = 'top', expand = 1, fill = 'both') + f.pack(side=BOTTOM,expand=NO,fill=X) + Label(self,text="Red=[1,0,0] Green=[0,1,0] Blue=[0,0,1]").pack(side=BOTTOM,expand=NO,fill=X) + + + def go(self): + """ + Allow the toplevel to return a handle for changing data + """ + self.o.mainloop() + + def goaway(self): + print "Called goaway" + self.o.destroy() + self.destroy() + print "Ought to be gone now..." + + def changedata(self,xyz): + self.xyz=xyz.copy() + self.npeaks=xyz.shape[0] + glDisableClientState(GL_VERTEX_ARRAY) + glVertexPointerd(self.xyz) + glEnableClientState(GL_VERTEX_ARRAY) + self.o.tkRedraw() + + def setps(self): + self.pointsize=float(self.ps.get()) + self.o.tkRedraw() + + + def redraw(self,o): + glClearColor(0., 0., 0., 0) + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + glOrtho(-1,1,-1,1,-1,1) + glDisable(GL_LIGHTING) + glColor3f(1.0, 1.0, 1.0) # white + glPointSize(self.pointsize) +# glEnable(GL_POINT_SMOOTH) + glDrawArrays(GL_POINTS, 0, self.npeaks ) + glBegin(GL_LINE_LOOP) + glColor3f(1.0, 0.0, 0.0) # red + glVertex3f(0.,0.,0.) + glVertex3f(1.,0.,0.) + glEnd() + glBegin(GL_LINE_LOOP) + glColor3f(0.0, 1.0, 0.0) # green + glVertex3f(0.,0.,0.) + glVertex3f(0.,1.,0.) + glEnd() + glBegin(GL_LINE_LOOP) + glColor3f(0.0, 0.0, 1.0) # blue + glVertex3f(0.,0.,0.) + glVertex3f(0.,0.,1.) + glEnd() + glEnable(GL_LIGHTING) + + + + + + + +if __name__=="__main__": + + try: + lines=open(sys.argv[1],"r").readlines() + except: + print "Usage %s gvector_file"%(sys.argv[0]) + sys.exit() + + on=0 + xyz=[] + for line in lines: + if on==1: + try: + vals=[float(x) for x in line.split()] + xyz.append( [ vals[0],vals[1],vals[2] ]) + except: + pass + if line.find("xr yr zr")>0: + on=1 + + npeaks = len(xyz) + + xyz=array(xyz) + main(xyz) + diff --git a/ImageD11/simplex.py b/ImageD11/simplex.py new file mode 100644 index 00000000..3f160d75 --- /dev/null +++ b/ImageD11/simplex.py @@ -0,0 +1,283 @@ +#!/usr/bin/env python +# +# Copyright (c) 2001 Vivake Gupta (v@omniscia.org). All rights reserved. +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License as +# published by the Free Software Foundation; either version 2 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +# USA +# +# This software is maintained by Vivake (v@omniscia.org) and is available at: +# http://www.omniscia.org/~vivake/python/Simplex.py + +# Modified (debugged?) 7/16/2004 Michele Vallisneri (vallis@vallis.org) + +# Copied into ImageD11 by Jon Wright - This works really really well!!! + +""" Simplex - a regression method for arbitrary nonlinear function minimization + +Simplex minimizes an arbitrary nonlinear function of N variables by the +Nelder-Mead Simplex method as described in: + +Nelder, J.A. and Mead, R., "A Simplex Method for Function Minimization", + Computer Journal, Vol. 7, 1965, pp. 308-313. + +It makes no assumptions about the smoothness of the function being minimized. +It converges to a local minimum which may or may not be the global minimum +depending on the initial guess used as a starting point. +""" + +import math +import copy +import sys + +class Simplex: + def __init__(self, testfunc, guess, increments, kR = -1, kE = 2, kC = 0.5): + """Initializes the simplex. + INPUTS + ------ + testfunc the function to minimize + guess[] an list containing initial guesses + increments[] an list containing increments, perturbation size + kR reflection constant (alpha =-1.0) + kE expansion constant (gamma = 2.0) + kC contraction constant (beta = 0.5) + """ + self.testfunc = testfunc + self.guess = guess + self.increments = increments + self.kR = kR + self.kE = kE + self.kC = kC + self.numvars = len(self.guess) + self.simplex = [] + + self.lowest = -1 + self.highest = -1 + self.secondhighest = -1 + + self.errors = [] + self.currenterror = 0 + + # Initialize vertices + # MV: the first vertex is just the initial guess + # the other N vertices are the initial guess plus the individual increments + # the last two vertices will store the centroid and the reflected point + # the compute errors at the ... vertices + + for vertex in range(0, self.numvars + 3): + self.simplex.append(copy.copy(self.guess)) + + for vertex in range(0, self.numvars + 1): + for x in range(0, self.numvars): + if x == (vertex - 1): + self.simplex[vertex][x] = self.guess[x] + self.increments[x] + self.errors.append(0) + self.calculate_errors_at_vertices() + + def minimize(self, epsilon = 0.0001, maxiters = 250, monitor = 1): + """Walks to the simplex down to a local minima. + INPUTS + ------ + epsilon convergence requirement + maxiters maximum number of iterations + monitor if non-zero, progress info is output to stdout + + OUTPUTS + ------- + an array containing the final values + lowest value of the error function + number of iterations taken to get here + """ + + iter = 0 + + for iter in range(0, maxiters): + # Identify highest and lowest vertices + + self.highest = 0 + self.lowest = 0 + for vertex in range(1, self.numvars + 1): + if self.errors[vertex] > self.errors[self.highest]: + self.highest = vertex + if self.errors[vertex] < self.errors[self.lowest]: + self.lowest = vertex + + # Identify second-highest vertex + + self.secondhighest = self.lowest + for vertex in range(0, self.numvars + 1): + if vertex == self.highest: + continue + elif vertex == self.secondhighest: + continue + elif self.errors[vertex] > self.errors[self.secondhighest]: + self.secondhighest = vertex + + # Test for convergence: + # compute the average merit figure (ugly) + + S = 0.0 + for vertex in range(0, self.numvars + 1): + S = S + self.errors[vertex] + F2 = S / (self.numvars + 1) + + # compute the std deviation of the merit figures (ugly) + + S1 = 0.0 + for vertex in range(0, self.numvars + 1): + S1 = S1 + (self.errors[vertex] - F2)**2 + T = math.sqrt(S1 / self.numvars) + + # Optionally, print progress information + + if monitor: + print '\r' + 72 * ' ', + print '\rIteration = %d Best = %f Worst = %f' % \ + (iter,self.errors[self.lowest],self.errors[self.highest]), + sys.stdout.flush() + + if T <= epsilon: + # We converged! Break out of loop! + + break; + else: + # Didn't converge. Keep crunching. + + # Calculate centroid of simplex, excluding highest vertex + # store centroid in element N+1 + + # loop over coordinates + for x in range(0, self.numvars): + S = 0.0 + for vertex in range(0, self.numvars + 1): + if vertex == self.highest: + continue + S = S + self.simplex[vertex][x] + self.simplex[self.numvars + 1][x] = S / self.numvars + + # reflect the simplex across the centroid + # store reflected point in elem. N + 2 (and self.guess) + + self.reflect_simplex() + self.currenterror = self.testfunc(self.guess) + + if self.currenterror < self.errors[self.highest]: + self.accept_reflected_point() + + if self.currenterror <= self.errors[self.lowest]: + self.expand_simplex() + self.currenterror = self.testfunc(self.guess) + + # at this point we can assume that the highest + # value has already been replaced once + if self.currenterror < self.errors[self.highest]: + self.accept_expanded_point() + elif self.currenterror >= self.errors[self.secondhighest]: + # worse than the second-highest, so look for + # intermediate lower point + + self.contract_simplex() + self.currenterror = self.testfunc(self.guess) + + if self.currenterror < self.errors[self.highest]: + self.accept_contracted_point() + else: + self.multiple_contract_simplex() + + # Either converged or reached the maximum number of iterations. + # Return the lowest vertex and the currenterror. + + for x in range(0, self.numvars): + self.guess[x] = self.simplex[self.lowest][x] + self.currenterror = self.errors[self.lowest] + return self.guess, self.currenterror, iter + + # same as expand, but with alpha < 1; kC = 0.5 fine with NR + + def contract_simplex(self): + for x in range(0, self.numvars): + self.guess[x] = self.kC * self.simplex[self.highest][x] + (1 - self.kC) * self.simplex[self.numvars + 1][x] + return + + # expand: if P is vertex and Q is centroid, alpha-expansion is Q + alpha*(P-Q), + # or (1 - alpha)*Q + alpha*P; default alpha is 2.0; agrees with NR + def expand_simplex(self): + for x in range(0, self.numvars): + self.guess[x] = self.kE * self.guess[x] + (1 - self.kE) * self.simplex[self.numvars + 1][x] + return + + # reflect: if P is vertex and Q is centroid, reflection is Q + (Q-P) = 2Q - P, + # which is achieved for kR = -1 (default value); agrees with NR + def reflect_simplex(self): + # loop over variables + for x in range(0, self.numvars): + self.guess[x] = self.kR * self.simplex[self.highest][x] + (1 - self.kR) * self.simplex[self.numvars + 1][x] + # store reflected point in elem. N + 2 + self.simplex[self.numvars + 2][x] = self.guess[x] + return + + # multiple contraction: around the lowest point; agrees with NR + + def multiple_contract_simplex(self): + for vertex in range(0, self.numvars + 1): + if vertex == self.lowest: + continue + for x in range(0, self.numvars): + self.simplex[vertex][x] = 0.5 * (self.simplex[vertex][x] + self.simplex[self.lowest][x]) + self.calculate_errors_at_vertices() + return + + def accept_contracted_point(self): + self.errors[self.highest] = self.currenterror + for x in range(0, self.numvars): + self.simplex[self.highest][x] = self.guess[x] + return + + def accept_expanded_point(self): + self.errors[self.highest] = self.currenterror + for x in range(0, self.numvars): + self.simplex[self.highest][x] = self.guess[x] + return + + def accept_reflected_point(self): + self.errors[self.highest] = self.currenterror + for x in range(0, self.numvars): + self.simplex[self.highest][x] = self.simplex[self.numvars + 2][x] + return + + def calculate_errors_at_vertices(self): + for vertex in range(0, self.numvars + 1): + if vertex == self.lowest: + # compute the error unless we're the lowest vertex + continue + for x in range(0, self.numvars): + self.guess[x] = self.simplex[vertex][x] + self.currenterror = self.testfunc(self.guess) + self.errors[vertex] = self.currenterror + return + +def myfunc(args): + return abs(args[0] * args[0] * args[0] * 5 - args[1] * args[1] * 7 + math.sqrt(abs(args[0])) - 118) + +def main(): + s = Simplex(myfunc, [1, 1, 1], [2, 4, 6]) + values, err, iter = s.minimize() + print 'args = ', values + print 'error = ', err + print 'iterations = ', iter + +if __name__ == '__main__': + main() + + diff --git a/ImageD11/transform.py b/ImageD11/transform.py new file mode 100644 index 00000000..041438f3 --- /dev/null +++ b/ImageD11/transform.py @@ -0,0 +1,170 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +from Numeric import * +import math as m + +def degrees(x): + """Convenience function""" + return x*180.0/math.pi + +def radians(x): + """Convenience function""" + return x*math.pi/180.0 + + +def compute_tth_eta(peaks,yc,ys,ty,zc,zs,tz,dist): + """ + Peaks is a 2 d array of x,y + yc is the centre in y + ys is the y pixel size + ty is the tilt around y + zc is the centre in z + zs is the z pixel size + tz is the tilt around z + dist is the sample - detector distance + """ + r1 = array( [ [ m.cos(tz) , m.sin(tz) , 0 ], + [ -m.sin(tz), m.cos(tz) , 0 ], + [ 0 , 0 , 1 ]],Float) + r2 = array( [ [ m.cos(ty) , 0 , m.sin(ty) ], + [ 0 , 1 , 0 ], + [-m.sin(ty) , 0 , m.cos(ty) ]],Float) + r2r1=matrixmultiply(r1,r2) + vec =array( [ zeros(peaks.shape[1]) , # place detector at zero, sample at -dist + (peaks[0,:]-yc)*ys , # x in search + (peaks[1,:]-zc)*zs ] , Float) # y in search + #print vec.shape + rotvec=matrixmultiply(r2r1,vec) + magrotvec=sqrt(sum(rotvec*rotvec,0)) + #print magrotvec.shape + eta=degrees(arctan2(rotvec[2,:],rotvec[1,:])) # bugger this one for tilts + # cosine rule a2 = b2+c2-2bccos(A) + # a is distance from (0,0,0) to rotvec => magrotvec + # b is distance from sample to detector + # c is distance from sample to pixel + a=magrotvec # 1D + b=dist # 0D + c=rotvec # 3D + #print c.shape + c[0,:]=c[0,:]-dist# 3D + c=sqrt(sum(c*c,0))# 1D + #print a.shape,c.shape + # print a.shape,c.shape + costwotheta = (b*b + c*c - a*a)/2/b/c + # print self.costwotheta.shape + twothetarad=arccos(costwotheta) + twotheta=degrees(twothetarad) + return twotheta, eta + + + +def compute_g_vectors(tth,eta,omega,wavelength): + """ + Generates spot positions in reciprocal space from twotheta, wavelength, omega and eta + Assumes single axis vertical + """ + tth=radians(tth) + eta=radians(eta) + om =radians(omega) + # Compute k vector - the scattering in the laboratory + c=cos(tth/2) # cos theta + s=sin(tth/2) # sin theta + ds=2*s/wavelength + k=zeros((3,tth.shape[0]),Float) + # x - along incident beam + k[0,:] = -ds*s + # y - towards door + k[1,:] = -ds*c*sin(eta) # plus or minus sin or cos + # z - towards roof + k[2,:] = ds*c*cos(eta) + # G-vectors - rotate k onto the crystal axes + g=zeros((3,ds.shape[0]),Float) + # + # g = R . k where: + # R = ( cos(omega) , sin(omega), 0 ) + # (-sin(omega) , cos(omega), 0 ) + # ( 0 , 0 , 1 ) + g[0,:] = cos(om)*k[0,:]+sin(om)*k[1,:] + g[1,:] =-sin(om)*k[0,:]+cos(om)*k[1,:] + g[2,:] = k[2,:] + return g + + +def uncompute_g_vectors(gv,wavelength): + """ + Given g-vectors compute tth,eta,omega + assert uncompute_g_vectors(compute_g_vector(tth,eta,omega))==tth,eta,omega + """ + # k=array(g.shape,Float) + # k[2,:]=g[2,:] # Z component in laboratory + # Length gives two-theta + gv=gv.copy().astype(Float) + ds = sqrt(sum(gv*gv,0)) + s = ds*wavelength/2 # sin theta + # Two theta from Bragg's law + try: + th = arcsin(s) + except: + print "Problem getting theta from sin(theta)" + print "Maximum and minimum in sin(theta)=",maximum.reduce(s),minimum.reduce(s) + c = cos(th) # cos theta + # Two theta in degrees + k=zeros(gv.shape,Float) + k[2,:]=gv[2,:] # Last line above == ds*c*cos(eta) + k[0,:] = -ds*s # From above + # Should check if this is true or not? + bottom=ds*c + try: + coseta = gv[2,:]/bottom + except: # Either a (0,0,0) peak or cos(theta) = 0 -> 90 degrees + ind=compress(bottom<1e-6,range(ds.shape[0])) + put(bottom,ind,1) + coseta = gv[2,:]/bottom + put(coseta,ind,1.) + #print "max and min in coseta",maximum.reduce(coseta),minimum.reduce(coseta) + coseta=clip(coseta,-1,1) + tth = degrees(2*th) + eta1 = degrees(arccos(coseta)) + eta2 = -eta1 + k[1,:] = -ds*c*sin(radians(eta1)) # Known also + # Know k and g for each peak - need to solve for omega + omega_crystal = arctan2(gv[0,:],gv[1,:]) + omega_laboratory = arctan2(k[0,:],k[1,:]) + omega1 = degrees( omega_crystal - omega_laboratory ) + k[1,:] = -ds*c*sin(radians(eta2)) # Known also + # Know k and g for each peak - need to solve for omega + omega_crystal = arctan2(gv[0,:],gv[1,:]) + omega_laboratory = arctan2(k[0,:],k[1,:]) + omega2 = degrees( omega_crystal - omega_laboratory ) + return tth,[eta1,eta2],[omega1,omega2] + +def uncompute_one_g_vector(gv,wavelength): + """ + Given g-vectors compute tth,eta,omega + assert uncompute_g_vectors(compute_g_vector(tth,eta,omega))==tth,eta,omega + """ + t,e,o=uncompute_g_vectors(transpose(array([gv,gv])),wavelength) + return t[0],[e[0][0],e[1][0]],[o[0][0],o[1][0]] + + + + diff --git a/ImageD11/twodplot.py b/ImageD11/twodplot.py new file mode 100644 index 00000000..3f81c0d2 --- /dev/null +++ b/ImageD11/twodplot.py @@ -0,0 +1,308 @@ + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +# This was broadly copied from an example in matplotlib, but then +# modified extensively. +# Should go back to something cleaner later... - JPW 2005 + + +""" +From the matplotlib examples - modified for mouse +""" +import matplotlib +matplotlib.use('TkAgg') +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg +from matplotlib.figure import Figure + +import Tkinter as Tk +import tkFileDialog, tkMessageBox +import sys,os,time + +class data: + def __init__(self,x,y,d={}): + self.x=x + self.y=y + self.d=d + +class twodplot(Tk.Frame): + def __init__(self,parent=None,data=None): + Tk.Frame.__init__(self,parent) + self.f = Figure(figsize=(8,5), dpi=100) + self.a = self.f.add_subplot(111) + self.plotitems={} + self.maxpoints=1000 # work around slow plotting + # print data + if data != None: + self.plotitems[data[0]]=data[1] + self.title=None + self.xr=self.yr=None + # a tk.DrawingArea + self.canvas = FigureCanvasTkAgg(self.f, master=self) + self.canvas.show() + self.tkc=self.canvas.get_tk_widget() + self.tkc.pack(side=Tk.TOP, fill=Tk.BOTH, expand=1) + self.tkc.bind("",self.on_down) + self.tkc.bind("",self.on_up) + self.tkc.bind("",self.on_move) + self.tkc.bind("",self.on_2) + self.tkc.bind("",self.on_3) + self.bindkeys() + self.rubberbandbox=None + self.bf1=Tk.Frame(self) + Tk.Button(master=self.bf1, text='Save Plot', command=self.printplot).pack(side=Tk.LEFT) + Tk.Button(master=self.bf1, text='LogY', command=self.logy).pack(side=Tk.LEFT) + Tk.Button(master=self.bf1, text='LogX', command=self.logx).pack(side=Tk.LEFT) +# FIXME - buttons for panx/y and zoomx/y + Tk.Button(master=self.bf1,text='>' , + command=lambda : self.keypress(self.a.panx,1 ) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf1,text='<', + command=lambda : self.keypress(self.a.panx,-1) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf1,text='^' , + command=lambda : self.keypress(self.a.pany,-1 ) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf1,text='v' , + command=lambda : self.keypress(self.a.pany,1) ).pack(side=Tk.LEFT) + self.bf1.pack(side=Tk.BOTTOM) + self.bf2=Tk.Frame(self) + Tk.Button(master=self.bf2,text='UnZoomX' , + command=lambda : self.keypress(self.a.zoomx,-1 ) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf2,text='ZoomX', + command=lambda : self.keypress(self.a.zoomx,1) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf2,text='ZoomY' , + command=lambda : self.keypress(self.a.zoomy,1 ) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf2,text='UnZoomY' , + command=lambda : self.keypress(self.a.zoomy,-1) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf2,text='Autoscale' , + command=lambda : self.keypress(self.autoscale) ).pack(side=Tk.LEFT) + Tk.Button(master=self.bf2,text='Autoscale Y', + command=lambda : self.keypress(self.autoscaley, e ) ).pack(side=Tk.LEFT) + self.bf2.pack(side=Tk.BOTTOM) + self.label=Tk.Label(self,text="Plot window messages") + self.label.pack(side=Tk.BOTTOM,fill=Tk.X, expand=0) + self.pack(side=Tk.TOP,fill=Tk.BOTH,expand=Tk.YES) + self.hidden=[] + self.replot() + self.xd=None + self.yd=None + + def printplot(self): + fn = tkFileDialog.asksaveasfilename(title="File name to print to", + defaultextension="png") + print fn + f,e = os.path.splitext(fn) + extns=['png','ps','eps','bmp','raw','rgb'] + print e + if e.lower() in ['.png','.ps','.eps','.bmp','.raw','.rgb']: + self.update_idletasks() # Try to get screen redrawn + self.canvas.print_figure(fn, dpi=300, orientation='landscape') + else: + tkMessageBox.showinfo("Sorry","I can only make output in these formats"+str(extns)) + + def keypress(self,*arg): + if len(arg)>1: + arg[0](*arg[1:]) + else: + arg[0]() + self.canvas.show() + + + def bindkeys(self): + return + self.bind_all('' ,lambda e: self.keypress(self.a.panx,1 ) ) + self.bind_all('',lambda e: self.keypress(self.a.panx,-1) ) + self.bind_all('' ,lambda e: self.keypress(self.a.pany,-1 ) ) + self.bind_all('' ,lambda e: self.keypress(self.a.pany,1) ) + self.bind_all('' ,lambda e: self.keypress(self.a.zoomx,-1 ) ) + self.bind_all('',lambda e: self.keypress(self.a.zoomx,1) ) + self.bind_all('' ,lambda e: self.keypress(self.a.zoomy,1 ) ) + self.bind_all('' ,lambda e: self.keypress(self.a.zoomy,-1) ) + self.bind_all('' , lambda e : self.keypress(self.autoscale) ) + self.bind_all('', lambda e : self.keypress(self.autoscaley, e ) ) + + def autoscaley(self,e): + print dir(self.a.dataLim) + print self.a.dataLim + yr=self.a.get_ylim() + self.a.set_ylim(yr) + + def adddata(self,data): + """ + Takes a tuple of name, data object + """ + self.plotitems[data[0]]=data[1] + if data[0] in self.hidden: + self.hidden.remove(data[0]) + self.replot() + + def hideall(self): + self.xr = self.yr = None + for item in self.plotitems.keys(): + self.hidden.append(item) + + + def removedata(self,name): + try: + self.plotitems.pop(name) + except KeyError: + pass + + + def replot(self): + self.a.clear() + if self.title!= None: + self.a.set_title(self.title) +# b : blue +# g : green +# r : red +# c : cyan +# m : magenta +# y : yello + c = ['g','r','c','m','y','b'] + for name in self.plotitems.keys(): + if name in self.hidden: + continue + item=self.plotitems[name] + print self.plotitems[name].d + if item.d.has_key('color'): + pc=item.d['color'] + else: + c.append(c[0]) + pc=c.pop(0) + if item.d.has_key('pointtype'): + pc=item.d['pointtype']+pc + else: + pc="."+pc + if item.d.has_key("xlabel"): + self.a.set_xlabel(item.d["xlabel"]) + if item.d.has_key("ylabel"): + self.a.set_ylabel(item.d["ylabel"]) + if item.d.has_key("title"): + self.a.set_title(item.d["title"]) + if item.x.shape[0]>self.maxpoints: + if tkMessageBox.askyesno("Slow plotting workaround","Shall I plot only the first %d points for increased speed?"%(self.maxpoints)): + ret = self.a.plot(item.x[:self.maxpoints],item.y[:self.maxpoints],pc) + else: + ret = self.a.plot(item.x,item.y,pc) + else: + ret = self.a.plot(item.x,item.y,pc) + if self.xr!=None: + self.a.set_xlim(self.xr) + if self.yr!=None: + self.a.set_ylim(self.yr) + self.canvas.show() + + def logy(self): +# FIXME - clip negative values before making logscaled? + if self.a.yaxis.is_log(): + self.a.set_yscale('linear') + else: + self.a.set_yscale('log') + self.canvas.show() + + def logx(self): +# FIXME - clip negative values before making logscaled? + if self.a.xaxis.is_log(): + self.a.set_xscale('linear') + else: + self.a.set_xscale('log') + self.canvas.show() + + def on_3(self,event): + self.autoscale() + + def autoscale(self): + self.a.cla() + self.xr = self.yr = None + self.replot() + + def on_2(self,event): + height = self.f.bbox.height() + x, y = event.x, height-event.y + (xd,yd)= self.a.transData.inverse_xy_tup( (x,y) ) + self.label.config(text="Clicked at x=%f y=%f"%(xd,yd)) + + # Callback functions for mouse + def on_down(self,event): + # get the x and y coords, flip y from top to bottom + height = self.f.bbox.height() + x, y = event.x, height-event.y + # transData transforms data coords to display coords. Use the + # inverse method to transform back + (self.xd,self.yd)= self.a.transData.inverse_xy_tup( (x,y) ) + # print "print mouse down at", t, val + # rubber banding: + if self.rubberbandbox!=None: self.tkc.delete(self.rubberbandbox) + self.startx=self.tkc.canvasx(event.x) + self.starty=self.tkc.canvasx(event.y) + + def on_move(self,event): + x = self.tkc.canvasx(event.x) + y = self.tkc.canvasy(event.y) + if (self.startx != event.x) and (self.starty != event.y) : + if self.rubberbandbox!=None: + self.tkc.delete(self.rubberbandbox) + self.rubberbandbox = self.tkc.create_rectangle(self.startx, self.starty, x, y, outline='green') + # this flushes the output, making sure that + # the rectangle makes it to the screen + # before the next event is handled + + def on_up(self,event): + # get the x and y coords, flip y from top to bottom + if self.xd==None: + return + self.tkc.delete(self.rubberbandbox) + height = self.f.bbox.height() + x, y = event.x, height-event.y + # transData transforms data coords to display coords. Use the + # inverse method to transform back + (self.xu,self.yu) = self.a.transData.inverse_xy_tup( (x,y) ) + if self.xu != self.xd and self.yu != self.yd: + # rescale + xr=[self.xd,self.xu];xr.sort() + yr=[self.yd,self.yu];yr.sort() + self.xr=xr + self.yr=yr + self.a.set_xlim(xr) + self.a.set_ylim(yr) + self.canvas.show() + + +if __name__=="__main__": + import epffile, powbase, mcadata, ciidata + if len(sys.argv)<3: + print "Usage: %s filename format"%(sys.argv[0]) + x=arange(0.0,3.0,0.01) + dat=epffile.powderdata(x,sin(2*pi*x)+5,sqrt(sin(2*pi*x)+5),{ "title":"sin x" } ) + else: + try: + if sys.argv[2]=="powbase": + dat=powbase.powbase(sys.argv[1]) + if sys.argv[2]=="epf": + dat=epffile.epffile(sys.argv[1]) + if sys.argv[2]=="mca": + dat=mcadata.mcadata(sys.argv[1]) + except: + print "Could not read your file %s" % (sys.argv[1]) + raise + + root = Tk.Tk() + root.wm_title("Two dimensional plotting") + p=twodplot(root,data=dat) + + Tk.mainloop() + diff --git a/ImageD11/unitcell.py b/ImageD11/unitcell.py new file mode 100644 index 00000000..98225bf6 --- /dev/null +++ b/ImageD11/unitcell.py @@ -0,0 +1,333 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +# +# +# + + +from Numeric import * +from LinearAlgebra import inverse +import math + +def radians(x): + return x*math.pi/180. + +def degrees(x): + return x*180./math.pi + + + +def cross(a,b): + """ + a x b has length |a||b|sin(theta) + """ + return array([ a[1]*b[2]-a[2]*b[1] ,a[2]*b[0]-b[2]*a[0], a[0]*b[1]-b[0]*a[1] ],Float) + +def unit(a): + """ + Normalise vector a to unit length + """ + return a/sqrt(dot(a,a)) + + + + +# Systematic absences + +def P(h,k,l): + return False + +def A(h,k,l): + return (k+l)%2 != 0 + +def B(h,k,l): + return (h+l)%2 != 0 + +def C(h,k,l): + return (h+k)%2 != 0 + +def I(h,k,l): + return (h+k+l)%2 != 0 + +def F(h,k,l): + return (h+k)%2!=0 or (h+l)%2!=0 or (k+l)%2!=0 + +outif = { + "P" : P , + "A" : I , + "B" : B , + "C" : C , + "I" : I , + "F" : F } + +def cellfromstring(s): + items=s.split() + # print items + latt = [float(x) for x in items[0:6]] + try: + symm = items[6] + except IndexError: + symm = 'P' + return unitcell(latt,symm) + +class unitcell: + # Unit cell stuff + # Generate a list of peaks from a unit cell + def __init__(self,lattice_parameters,symmetry="P",verbose=0): + """ + Unit cell class + supply a list (tuple etc) of a,b,c,alpha,beta,gamma + optionally a symmetry, one of "P","A","B","C","I","F" + """ + self.lattice_parameters=array(lattice_parameters) + if self.lattice_parameters.shape[0]!=6: + raise "You must supply 6 lattice parameters, a,b,c,alpha,beta,gamma" + self.symmetry=symmetry + if self.symmetry not in ["P","A","B","C","I","F"]: + raise "Your symmetry "+self.symmetry+" was not recognised" + # assigning a function here! + self.absent=outif[self.symmetry] + a = self.lattice_parameters[0] + b = self.lattice_parameters[1] + c = self.lattice_parameters[2] + self.alpha=radians(self.lattice_parameters[3]) + ca= math.cos(radians(self.lattice_parameters[3])) + cb= math.cos(radians(self.lattice_parameters[4])) + cg= math.cos(radians(self.lattice_parameters[5])) + if verbose==1: print "Unit cell",self.lattice_parameters + self.g = array( [[ a*a , a*b*cg, a*c*cb ], + [ a*b*cg , b*b , b*c*ca ], + [ a*c*cb , b*c*ca, c*c ]],Float) + if verbose==1: print "Metric tensor\n",self.g + try: + self.gi = inverse(self.g) + except: + raise "Unit cell was degenerate, could not determine reciprocal metric tensor" + if verbose==1: print "Reciprocal Metric tensor\n",self.gi + self.as=sqrt(self.gi[0,0]) + self.bs=sqrt(self.gi[1,1]) + self.cs=sqrt(self.gi[2,2]) + + self.alphas=degrees(math.acos(self.gi[1,2]/self.bs/self.cs)) + self.betas =degrees(math.acos(self.gi[0,2]/self.as/self.cs)) + self.gammas=degrees(math.acos(self.gi[0,1]/self.as/self.bs)) + if verbose==1: print "Reciprocal cell" + if verbose==1: print self.as,self.bs,self.cs,self.alphas,self.betas,self.gammas + # Equation 3 from Busing and Levy + self.B = array ( [ [ self.as , self.bs*cos(radians(self.gammas)) , self.cs*cos(radians(self.betas)) ] , + [ 0 , self.bs*sin(radians(self.gammas)) , -self.cs*sin(radians(self.betas))*ca ], + [ 0 , 0 , 1./c ] ] , Float) + if verbose==1: print self.B + if verbose==1: print matrixmultiply(transpose(self.B),self.B)-self.gi # this should be zero + self.hkls=None + self.peaks=None + self.limit=0 + + def tostring(self): + """ + Write out a line containing unit cell information + """ + return "%f %f %f %f %f %f %s"%(self.lattice_parameters[0], + self.lattice_parameters[1], + self.lattice_parameters[2], + self.lattice_parameters[3], + self.lattice_parameters[4], + self.lattice_parameters[5], + self.symmetry) + + def anglehkls(self,h1,h2): + """ + Compute the angle between reciprocal lattice vectors h1, h2 + """ + g1 = dot(h1,matrixmultiply(self.gi,h1)) + g2 = dot(h2,matrixmultiply(self.gi,h2)) + g12= dot(h1,matrixmultiply(self.gi,h2)) + costheta = g12/sqrt(g1*g2) + try: + return degrees(math.acos(costheta)),costheta + except: + if abs(costheta-1) < 1e-6: + return 0.,1.0 + if abs(costheta+1) < 1e-6: + return 180.,-1.0 + print "Error in unit cell class determining angle" + print "h1",h1,"h2",h2,"Costheta=",costheta + raise + + def gethkls(self,dsmax): + """ + Generate hkl list + Argument dsmax is the d* limit (eg 1/d) + Default of zero gives only the (000) reflection + """ + if dsmax == self.limit and self.peaks!=None: + return self.peaks + h=k=0 + l=1 # skip 0,0,0 + hs=ks=ls=1 + b=0 + peaks=[] + while abs(h)<20: # H + while abs(k)<20: # K + while abs(l)<20: #L + ds=self.ds([h,k,l]) + if ds < dsmax: + if not self.absent(h,k,l): + peaks.append([ds,(h,k,l)]) + else: + pass + b=0 + else: + if ls==1: + ls=-1 + l=0 + else: + ls=1 + l=0 + b=b+1 + break + l=l+ls + k=k+ks + # l is always zero here + if b>1: + if ks==1: + ks=-1 + k=-1 + else: + ks=1 + k=0 + b=b+1 + break + h=h+hs + if b>3: + if hs==1: + hs=-1 + h=-1 + else: + hs=1 + h=0 + break + + peaks.sort() + #for peak in peaks: + # print peak + self.peaks=peaks + self.limit=dsmax + return peaks + + def ds(self,h): + return math.sqrt(dot(h,matrixmultiply(self.gi,h))) # 1/d or d* + + + def makerings(self,limit,tol=0.001): + """ + Makes a list of computed powder rings + The tolerance is the difference in d* to decide + if two peaks overlap + """ + self.peaks=self.gethkls(limit) # [ ds, [hkl] ] + self.ringds=[] # a list of floats + self.ringhkls={} # a dict of lists of integer hkl + # Append first peak + peak = self.peaks[0] + self.ringds.append(peak[0]) + self.ringhkls[peak[0]] = [peak[1]] + for peak in self.peaks[1:]: + if abs(peak[0] - self.ringds[-1]) < tol: + self.ringhkls[self.ringds[-1]].append(peak[1]) + else: + self.ringds.append(peak[0]) + self.ringhkls[self.ringds[-1]]= [peak[1]] + + + def orient(self,ring1,g1,ring2,g2,verbose=0): + """ + Compute an orientation matrix using cell parameters and the indexing + of two reflections + + Orientation matrix: + A matrix such that h = A^1 g + define 3 vectors t1,t2,t3 + t1 is parallel to first peak (unit vector along g1) + t2 is in the plane of both (unit vector along g1x(g1xg2)) + t3 is perpendicular to both (unit vector along g1xg2) + """ + + costheta = dot(g1,g2)/sqrt(dot(g2,g2))/sqrt(dot(g1,g1)) + if verbose==1: print "observed costheta",costheta + best=5. + for ha in self.ringhkls[self.ringds[ring1]]: + for hb in self.ringhkls[self.ringds[ring2]]: + ca=self.anglehkls(ha,hb) + if abs(ca[1]-costheta) < best and ha!=hb: + h1=ha + h2=hb + best=ca[1] + if verbose==1: + print "Assigning h1",h1,g1,self.ds(h1),sqrt(dot(g1,g1)),self.ds(h1)-sqrt(dot(g1,g1)) + print "Assigning h2",h2,g2,self.ds(h2),sqrt(dot(g2,g2)),self.ds(h1)-sqrt(dot(g1,g1)) + print "Cos angle calc",self.anglehkls(h1,h2),"obs",costheta + h1c=matrixmultiply(self.B,h1) + h2c=matrixmultiply(self.B,h2) + t1c=unit(h1c) + t3c=unit(cross(h1c,h2c)) + t2c=unit(cross(h1c,t3c)) + t1g=unit(g1) + t3g=unit(cross(g1,g2)) + t2g=unit(cross(g1,t3g)) + T_g = transpose(array([t1g,t2g,t3g])) # Array are stored by rows and + T_c = transpose(array([t1c,t2c,t3c])) # these are columns + U=matrixmultiply(T_g , inverse(T_c)) + UB=matrixmultiply(U,self.B) + UBI=inverse(UB) + if verbose==1: + print "UBI" + print UBI + print "Grain gi" + print matrixmultiply(transpose(UB),UB) + print "Cell gi" + print self.gi + h=matrixmultiply(UBI,g1) + print "(%9.3f, %9.3f, %9.3f)"%(h[0],h[1],h[2]) + h=matrixmultiply(UBI,g2) + print "(%9.3f, %9.3f, %9.3f)"%(h[0],h[1],h[2]) + self.UBI=UBI + self.UB=UB + + + +if __name__=="__main__": + import sys,time + start=time.time() + cell = unitcell([float(x) for x in sys.argv[1:7]],sys.argv[7]) + p=cell.gethkls(float(sys.argv[8])) + n=len(p) + print "Generated",n,"peaks in",time.time()-start,"/s" + for k in p: + try: + print k[1],k[0],1/k[0],k[1] + except: + pass + + + + + diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..9051c247 --- /dev/null +++ b/LICENSE @@ -0,0 +1,325 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your + freedom to share and change it. By contrast, the GNU General Public + License is intended to guarantee your freedom to share and change free + software--to make sure the software is free for all its users. This + General Public License applies to most of the Free Software + Foundation's software and to any other program whose authors commit to + using it. (Some other Free Software Foundation software is covered by + the GNU Library General Public License instead.) You can apply it to + your programs, too. + + When we speak of free software, we are referring to freedom, not + price. Our General Public Licenses are designed to make sure that you + have the freedom to distribute copies of free software (and charge for + this service if you wish), that you receive source code or can get it + if you want it, that you can change the software or use pieces of it + in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid + anyone to deny you these rights or to ask you to surrender the rights. + These restrictions translate to certain responsibilities for you if you + distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether + gratis or for a fee, you must give the recipients all the rights that + you have. You must make sure that they, too, receive or can get the + source code. And you must show them these terms so they know their + rights. + + We protect your rights with two steps: (1) copyright the software, and + (2) offer you this license which gives you legal permission to copy, + distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain + that everyone understands that there is no warranty for this free + software. If the software is modified by someone else and passed on, we + want its recipients to know that what they have is not the original, so + that any problems introduced by others will not reflect on the original + authors' reputations. + + Finally, any free program is threatened constantly by software + patents. We wish to avoid the danger that redistributors of a free + program will individually obtain patent licenses, in effect making the + program proprietary. To prevent this, we have made it clear that any + patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and + modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains + a notice placed by the copyright holder saying it may be distributed + under the terms of this General Public License. The "Program", below, + refers to any such program or work, and a "work based on the Program" + means either the Program or any derivative work under copyright law: + that is to say, a work containing the Program or a portion of it, + either verbatim or with modifications and/or translated into another + language. (Hereinafter, translation is included without limitation in + the term "modification".) Each licensee is addressed as "you". + + Activities other than copying, distribution and modification are not + covered by this License; they are outside its scope. The act of + running the Program is not restricted, and the output from the Program + is covered only if its contents constitute a work based on the + Program (independent of having been made by running the Program). + Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's + source code as you receive it, in any medium, provided that you + conspicuously and appropriately publish on each copy an appropriate + copyright notice and disclaimer of warranty; keep intact all the + notices that refer to this License and to the absence of any warranty; + and give any other recipients of the Program a copy of this License + along with the Program. + + You may charge a fee for the physical act of transferring a copy, and + you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion + of it, thus forming a work based on the Program, and copy and + distribute such modifications or work under the terms of Section 1 + above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + + These requirements apply to the modified work as a whole. If + identifiable sections of that work are not derived from the Program, + and can be reasonably considered independent and separate works in + themselves, then this License, and its terms, do not apply to those + sections when you distribute them as separate works. But when you + distribute the same sections as part of a whole which is a work based + on the Program, the distribution of the whole must be on the terms of + this License, whose permissions for other licensees extend to the + entire whole, and thus to each and every part regardless of who wrote it. + + Thus, it is not the intent of this section to claim rights or contest + your rights to work written entirely by you; rather, the intent is to + exercise the right to control the distribution of derivative or + collective works based on the Program. + + In addition, mere aggregation of another work not based on the Program + with the Program (or with a work based on the Program) on a volume of + a storage or distribution medium does not bring the other work under + the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, + under Section 2) in object code or executable form under the terms of + Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + + The source code for a work means the preferred form of the work for + making modifications to it. For an executable work, complete source + code means all the source code for all modules it contains, plus any + associated interface definition files, plus the scripts used to + control compilation and installation of the executable. However, as a + special exception, the source code distributed need not include + anything that is normally distributed (in either source or binary + form) with the major components (compiler, kernel, and so on) of the + operating system on which the executable runs, unless that component + itself accompanies the executable. + + If distribution of executable or object code is made by offering + access to copy from a designated place, then offering equivalent + access to copy the source code from the same place counts as + distribution of the source code, even though third parties are not + compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program + except as expressly provided under this License. Any attempt + otherwise to copy, modify, sublicense or distribute the Program is + void, and will automatically terminate your rights under this License. + However, parties who have received copies, or rights, from you under + this License will not have their licenses terminated so long as such + parties remain in full compliance. + + 5. You are not required to accept this License, since you have not + signed it. However, nothing else grants you permission to modify or + distribute the Program or its derivative works. These actions are + prohibited by law if you do not accept this License. Therefore, by + modifying or distributing the Program (or any work based on the + Program), you indicate your acceptance of this License to do so, and + all its terms and conditions for copying, distributing or modifying + the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the + Program), the recipient automatically receives a license from the + original licensor to copy, distribute or modify the Program subject to + these terms and conditions. You may not impose any further + restrictions on the recipients' exercise of the rights granted herein. + You are not responsible for enforcing compliance by third parties to + this License. + + 7. If, as a consequence of a court judgment or allegation of patent + infringement or for any other reason (not limited to patent issues), + conditions are imposed on you (whether by court order, agreement or + otherwise) that contradict the conditions of this License, they do not + excuse you from the conditions of this License. If you cannot + distribute so as to satisfy simultaneously your obligations under this + License and any other pertinent obligations, then as a consequence you + may not distribute the Program at all. For example, if a patent + license would not permit royalty-free redistribution of the Program by + all those who receive copies directly or indirectly through you, then + the only way you could satisfy both it and this License would be to + refrain entirely from distribution of the Program. + + If any portion of this section is held invalid or unenforceable under + any particular circumstance, the balance of the section is intended to + apply and the section as a whole is intended to apply in other + circumstances. + + It is not the purpose of this section to induce you to infringe any + patents or other property right claims or to contest validity of any + such claims; this section has the sole purpose of protecting the + integrity of the free software distribution system, which is + implemented by public license practices. Many people have made + generous contributions to the wide range of software distributed + through that system in reliance on consistent application of that + system; it is up to the author/donor to decide if he or she is willing + to distribute software through any other system and a licensee cannot + impose that choice. + + This section is intended to make thoroughly clear what is believed to + be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in + certain countries either by patents or by copyrighted interfaces, the + original copyright holder who places the Program under this License + may add an explicit geographical distribution limitation excluding + those countries, so that distribution is permitted only in or among + countries not thus excluded. In such case, this License incorporates + the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions + of the General Public License from time to time. Such new versions will + be similar in spirit to the present version, but may differ in detail to + address new problems or concerns. + + Each version is given a distinguishing version number. If the Program + specifies a version number of this License which applies to it and "any + later version", you have the option of following the terms and conditions + either of that version or of any later version published by the Free + Software Foundation. If the Program does not specify a version number of + this License, you may choose any version ever published by the Free Software + Foundation. + + 10. If you wish to incorporate parts of the Program into other free + programs whose distribution conditions are different, write to the author + to ask for permission. For software which is copyrighted by the Free + Software Foundation, write to the Free Software Foundation; we sometimes + make exceptions for this. Our decision will be guided by the two goals + of preserving the free status of all derivatives of our free software and + of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY + FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN + OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES + PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED + OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS + TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE + PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, + REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING + WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR + REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, + INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING + OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED + TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY + YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER + PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE + POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest + possible use to the public, the best way to achieve this is to make it + free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest + to attach them to the start of each source file to most effectively + convey the exclusion of warranty; and each file should have at least + the "copyright" line and a pointer to where the full notice is found. + + ImageD11_v0.4 Software for beamline ID11 + Copyright (C) 2005 Jon Wright + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + Also add information on how to contact you by electronic and paper mail. + + If the program is interactive, make it output a short notice like this + when it starts in an interactive mode: + + ImageD11 version 0.4, Copyright (C) 2005 Jon Wright + ImageD11 comes with ABSOLUTELY NO WARRANTY; for details select help, license. + This is free software, and you are welcome to redistribute it + under certain conditions; + + This General Public License does not permit incorporating your program into + proprietary programs. If your program is a subroutine library, you may + consider it more useful to permit linking proprietary applications with the + library. If this is what you want to do, use the GNU Library General + Public License instead of this License. diff --git a/MANIFEST b/MANIFEST new file mode 100644 index 00000000..79d91d98 --- /dev/null +++ b/MANIFEST @@ -0,0 +1,30 @@ +setup.py +ImageD11\__init__.py +ImageD11\bisplev.py +ImageD11\blobcorrector.py +ImageD11\data.py +ImageD11\guiindexer.py +ImageD11\guimaker.py +ImageD11\guipeaksearch.py +ImageD11\guitransformer.py +ImageD11\indexing.py +ImageD11\license.py +ImageD11\listdialog.py +ImageD11\opendata.py +ImageD11\peakmerge.py +ImageD11\plot3d.py +ImageD11\simplex.py +ImageD11\transform.py +ImageD11\twodplot.py +ImageD11\unitcell.py +src\closest.c +src\connectedpixels.c +src\splines.c +src\dset.h +src\bispev.f +scripts\ImageD11_gui.py +scripts\peaksearch.py +scripts\recoveromega.py +README.TXT +MANIFEST +LICENSE diff --git a/PKG-INFO b/PKG-INFO new file mode 100644 index 00000000..798f3146 --- /dev/null +++ b/PKG-INFO @@ -0,0 +1,10 @@ +Metadata-Version: 1.0 +Name: ImageD11 +Version: 0.4 +Summary: ImageD11 +Home-page: UNKNOWN +Author: Jon Wright +Author-email: wright@esrf.fr +License: GPL +Description: UNKNOWN +Platform: UNKNOWN diff --git a/README.TXT b/README.TXT new file mode 100644 index 00000000..bf7ecb4b --- /dev/null +++ b/README.TXT @@ -0,0 +1,17 @@ +ImageD11 +Version 0.4 +Jon Wright +wright@esrf.fr + +You probably want to try something like: +python setup.py build +or +python setup.py build --compiler=mingw32 + +Followed by: +python setup.py install + + +Good luck! + + diff --git a/scripts/ImageD11_gui.py b/scripts/ImageD11_gui.py new file mode 100644 index 00000000..028f8a9a --- /dev/null +++ b/scripts/ImageD11_gui.py @@ -0,0 +1,136 @@ +from ImageD11.guimaker import GuiMaker +from Tkinter import * + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + + +if __name__=="__main__": + def help(): + from tkMessageBox import showinfo + showinfo("Help","Sorry, no help for you!\nPlease try harder") + ImageD11credits = """Thanks to: + + All of the Fable team, which includes at least: + Andy Goetz, Gavin Vaughan, + Soren Schmidt, Henning Poulsen, Larry Margulies + ...and others who should remind me to mention them + + John Hunter for the matplotlib plotting + + Anyone who tests me and gives useful feedback + + Jon Wright, for writing me! + """ + from ImageD11 import license + license = license.license + + + def credits(): + from tkMessageBox import showinfo + showinfo("Credits",ImageD11credits) + + def showlicense(): + from ScrolledText import ScrolledText + win = ScrolledText(Toplevel(),width=100) + win.insert(END,license) +# print dir(win) + win.pack(expand=1,fill=BOTH) + win.focus_set() +# showinfo("License",license) + + import tkFileDialog,os + + class TestGuiMaker(GuiMaker): + def start(self): + from tkMessageBox import showinfo + + + startmessage = """ +ImageD11 version 0.4, Copyright (C) 2005 Jon Wright +ImageD11 comes with ABSOLUTELY NO WARRANTY; for details select help, license. +This is free software, and you are welcome to redistribute it under certain conditions + +Please send useful feedback to wright@esrf.fr +""" + + startmessage += """ +Stuff to do: + + Implement the image peaksearching in the gui (maybe display opengl images??) + Separate the peak merging/reading from the gui for standalone scripting + Same for the transformations - once parameters are known gui is not needed + Tidy the mulitple threshold consequences + Implement those filters based in intensity etc + + Sort out a file format which tracks all information properly? +""" + showinfo("Welcome to ImageD11_v0.4",startmessage) +# from tkMessageBox import showinfo +# showinfo("Sorry","You had to say yes then if you want to use the program") +# sys.exit() + + from ImageD11 import guipeaksearch + self.peaksearcher = guipeaksearch.guipeaksearcher(self) + self.finalpeaks=None + from ImageD11 import guitransformer + self.unitcell=None + self.gv=None + self.transformer = guitransformer.guitransformer(self) + from ImageD11 import guiindexer + self.indexer = guiindexer.guiindexer(self) + import sys + self.menuBar = [ ( "File", 0, + [ ( "Print" , 0, self.printplot ), + ( "Exit", 1, sys.exit) ] ) , + self.peaksearcher.menuitems, + self.transformer.menuitems, + self.indexer.menuitems, + ( "Plotting", 0, + [ ( "Autoscale", 0, self.autoscaleplot), + ( "Clear plot",0, self.clearplot), + ] ) , + + ( "Help", 0, + [ ( "Help Me!", 0, help) , + ( "Credits" , 0, credits) , + ( "License" , 0, showlicense) + ] ) ] + + def printplot(self): self.twodplotter.printplot() + + def autoscaleplot(self): + self.twodplotter.autoscale() + def clearplot(self): + self.twodplotter.plotitems={} + self.twodplotter.replot() + + def makeWidgets(self): + # Get size of TopLevels window + self.opener=tkFileDialog.Open() + self.saver=tkFileDialog.SaveAs() + from ImageD11 import twodplot + self.twodplotter=twodplot.twodplot(self) + self.twodplotter.pack(side=RIGHT, expand=1, fill=BOTH) + + root = Tk() + root.wm_title("ImageD11") + TestGuiMaker() + root.mainloop() diff --git a/scripts/peaksearch.py b/scripts/peaksearch.py new file mode 100644 index 00000000..e7102a22 --- /dev/null +++ b/scripts/peaksearch.py @@ -0,0 +1,115 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +import time +reallystart=time.time() +from math import sqrt +import sys,glob,os.path +from ImageD11 import blobcorrector +from ImageD11 import opendata +from ImageD11 import connectedpixels +import Numeric + +def peaksearch(filename, outputfile, corrector, blobim , thresholds): + """ + Searches for peaks + """ + t0=time.time() + f=open(outputfile,"aq") + data_object = opendata.openedf(filename) + picture = data_object.data + print "%s"%(filename), + f.write("\n\n# File %s\n"%(filename)) + f.write("# Processed on %s\n"%(time.asctime())) + f.write("# Spatial correction from %s\n"%(corrector.splinefile)) + f.write("# SPLINE X-PIXEL-SIZE %f\n"%(corrector.xsize)) + f.write("# SPLINE Y-PIXEL-SIZE %f\n"%(corrector.ysize)) + for item in data_object.header.keys(): + try: + f.write("# %s = %s\n"%(item,data_object.header[item])) + except KeyError: + pass + if blobim.shape != data_object.data.shape: + raise "Incompatible blobimage buffer for file %s"%(filename) + for threshold in thresholds: + npks=0 + np=connectedpixels.connectedpixels(picture,blobim,threshold,verbose=0) + npi,sum,sumsq,com0,com1,com00,com01,com11=connectedpixels.blobproperties(picture,blobim,np,verbose=0) + f.write("\n#Threshold level %f\n"%(threshold)) + f.write( "# Number_of_pixels Average_counts x y xc yc sig_x sig_y cov_xy\n") + + for i in range(len(npi)): + if npi[i]>1: # Throw out one pixel peaks (div zero) + npks=npks+1 + n = npi[i] + avg = sum[i]/n # Average intensity + si = sqrt((sumsq[i] - n*avg*avg)/(n-1.)) # Standard dev on intensity + c0 = com0[i]/sum[i] # Centre of mass in index 0 + c1 = com1[i]/sum[i] # Centre of mass in index 1 + try: + c00 = sqrt((com00[i]/sum[i] - c0*c0)) + except: + c00 = 0. # this means a line of pixels and rounding errors + try: + c11 = sqrt((com11[i]/sum[i] - c1*c1)) + except: + c11 = 0. + try: + c01 = (com01[i]/sum[i] - c0*c1)/c00/c11 + except: + c01=0. + # Spatial corrections : + c0c, c1c = corrector.correct(c0,c1) + s = "%d %f %f %f %f %f %f %f %f\n"%(n, avg, c0, c1, c0c, c1c, c00, c11, c01) + f.write(s) + print "T=%-5d n=%-5d;"%(int(threshold),npks), + f.close() + print " time %f/s"%(time.time()-t0) + return npks # Number of peaks found + + +if __name__=="__main__": + + try: + stem = sys.argv[1] + outfile = sys.argv[2] + corrector=blobcorrector.correctorclass(sys.argv[5]) + thresholds = [float(t) for t in sys.argv[6:]] + first = int(sys.argv[3]) + last = int(sys.argv[4]) + files = ["%s%04d%s"%(stem,i,".edf") for i in range(first,last+1)] + blobim=Numeric.zeros(opendata.openedf(files[0]).data.shape,Numeric.Int) + if len(files)==0: + raise "No files found for stem %s"%(stem) + files.sort() + start = time.time() + print "File being treated in -> out, elapsed time" + nt=len(thresholds) + for filein in files: + peaksearch(filein, outfile, corrector, blobim, thresholds) + except: + print "Usage: %s filename outputfile first last spline threshold [threshold...]" + raise + sys.exit() + +end=time.time() +t=end-reallystart +print "Total time = %f /s, per image = %f /s"%(t,t/len(files)) diff --git a/scripts/recoveromega.py b/scripts/recoveromega.py new file mode 100644 index 00000000..a8401196 --- /dev/null +++ b/scripts/recoveromega.py @@ -0,0 +1,57 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +import sys + +try: + logfile = open(sys.argv[1],"r") + pksfile = open(sys.argv[2],"r") + newpksfile = open(sys.argv[3],"w") +except: + print "Usage: %s logfile pksfile newpksfile"%(sys.argv[0]) + sys.exit() + +# Read all of the scans into a dictionary of filenames/omega angles + +lookups = {} + +for line in logfile.readlines(): + if line.find(".edf")>0: + try: + om = line.split()[0] + fullname=line.split()[1] + name=fullname.split('/')[-1] + lookups[name]=om + except: + print line + raise + +logfile.close() + +for line in pksfile.readlines(): + newpksfile.write(line) + if line.find("File ")>0: + name=line.split()[-1] + + newpksfile.write("# Omega = %s\n"%(lookups[name]) ) + + +pksfile.close() +newpksfile.close() diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..27635d93 --- /dev/null +++ b/setup.py @@ -0,0 +1,56 @@ + + + +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + + + +from distutils.core import setup, Extension + +cl = Extension("closest", sources=['src/closest.c']) +cp = Extension("connectedpixels", sources = ['src/connectedpixels.c']) +# Fortran hack +import os +os.system("g77 -c src/bispev.f -o src/bispev.o") +os.system("ar -cr src/libsplines.a src/bispev.o") +bl = Extension("_splines", sources = ['src/splines.c'], + libraries = ["splines"], library_dirs = ["src"] ) + + +setup(name='ImageD11', + version=0.4, + author='Jon Wright', + author_email='wright@esrf.fr', + description='ImageD11', + license = "GPL", + ext_package = "ImageD11", + ext_modules = [cl,cp,bl], + packages = ["ImageD11"], + scripts = ["scripts/peaksearch.py", "scripts/ImageD11_gui.py"]) + + + +"""py_modules = ["bisplev","blobcorrector","data","opendata", + "peakmerge","peaksearch","simplex", + "transform","unitcell","indexing", + "guiindexer","guimaker","guipeaksearch", + "guitransformer","listdialog","plot3d","twodplot"] + )""" + diff --git a/src/bispev.f b/src/bispev.f new file mode 100644 index 00000000..d7e52d01 --- /dev/null +++ b/src/bispev.f @@ -0,0 +1,1880 @@ +C +C Routines from netlib +C Copyright by Paul Dierckx +C +C + + + subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, + * iwrk,kwrk,ier) +c subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,... +c ,my a bivariate spline s(x,y) of degrees kx and ky, given in the +c b-spline representation. +c +c calling sequence: +c call bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, +c * iwrk,kwrk,ier) +c +c input parameters: +c tx : real array, length nx, which contains the position of the +c knots in the x-direction. +c nx : integer, giving the total number of knots in the x-direction +c ty : real array, length ny, which contains the position of the +c knots in the y-direction. +c ny : integer, giving the total number of knots in the y-direction +c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the +c b-spline coefficients. +c kx,ky : integer values, giving the degrees of the spline. +c x : real array of dimension (mx). +c before entry x(i) must be set to the x co-ordinate of the +c i-th grid point along the x-axis. +c tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx. +c mx : on entry mx must specify the number of grid points along +c the x-axis. mx >=1. +c y : real array of dimension (my). +c before entry y(j) must be set to the y co-ordinate of the +c j-th grid point along the y-axis. +c ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my. +c my : on entry my must specify the number of grid points along +c the y-axis. my >=1. +c wrk : real array of dimension lwrk. used as workspace. +c lwrk : integer, specifying the dimension of wrk. +c lwrk >= mx*(kx+1)+my*(ky+1) +c iwrk : integer array of dimension kwrk. used as workspace. +c kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my. +c +c output parameters: +c z : real array of dimension (mx*my). +c on succesful exit z(my*(i-1)+j) contains the value of s(x,y) +c at the point (x(i),y(j)),i=1,...,mx;j=1,...,my. +c ier : integer error flag +c ier=0 : normal return +c ier=10: invalid input data (see restrictions) +c +c restrictions: +c mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my +c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx +c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my +c +c other subroutines required: +c fpbisp,fpbspl +c +c references : +c de boor c : on calculating with b-splines, j. approximation theory +c 6 (1972) 50-62. +c cox m.g. : the numerical evaluation of b-splines, j. inst. maths +c applics 10 (1972) 134-149. +c dierckx p. : curve and surface fitting with splines, monographs on +c numerical analysis, oxford university press, 1993. +c +c author : +c p.dierckx +c dept. computer science, k.u.leuven +c celestijnenlaan 200a, b-3001 heverlee, belgium. +c e-mail : Paul.Dierckx@cs.kuleuven.ac.be +c +c latest update : march 1987 +c +c ..scalar arguments.. + integer nx,ny,kx,ky,mx,my,lwrk,kwrk,ier +c ..array arguments.. + integer iwrk(kwrk) + real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my), + * wrk(lwrk) +c ..local scalars.. + integer i,iw,lwest +c .. +c before starting computations a data check is made. if the input data +c are invalid control is immediately repassed to the calling program. + ier = 10 + lwest = (kx+1)*mx+(ky+1)*my + if(lwrk.lt.lwest) go to 100 + if(kwrk.lt.(mx+my)) go to 100 + if(mx-1) 100,30,10 + 10 do 20 i=2,mx + if(x(i).lt.x(i-1)) go to 100 + 20 continue + 30 if(my-1) 100,60,40 + 40 do 50 i=2,my + if(y(i).lt.y(i-1)) go to 100 + 50 continue + 60 ier = 0 + iw = mx*(kx+1)+1 + call fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),wrk(iw), + * iwrk(1),iwrk(mx+1)) + 100 return + end + subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly) +c ..scalar arguments.. + integer nx,ny,kx,ky,mx,my +c ..array arguments.. + integer lx(mx),ly(my) + real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my), + * wx(mx,kx+1),wy(my,ky+1) +c ..local scalars.. + integer kx1,ky1,l,l1,l2,m,nkx1,nky1 + real*8 arg,sp,tb,te +c ..local arrays.. + real*8 h(6) +c ..subroutine references.. +c fpbspl +c .. + kx1 = kx+1 + nkx1 = nx-kx1 + tb = tx(kx1) + te = tx(nkx1+1) + l = kx1 + l1 = l+1 + do 40 i=1,mx + arg = x(i) + if(arg.lt.tb) arg = tb + if(arg.gt.te) arg = te + 10 if(arg.lt.tx(l1) .or. l.eq.nkx1) go to 20 + l = l1 + l1 = l+1 + go to 10 + 20 call fpbspl(tx,nx,kx,arg,l,h) + lx(i) = l-kx1 + do 30 j=1,kx1 + wx(i,j) = h(j) + 30 continue + 40 continue + ky1 = ky+1 + nky1 = ny-ky1 + tb = ty(ky1) + te = ty(nky1+1) + l = ky1 + l1 = l+1 + do 80 i=1,my + arg = y(i) + if(arg.lt.tb) arg = tb + if(arg.gt.te) arg = te + 50 if(arg.lt.ty(l1) .or. l.eq.nky1) go to 60 + l = l1 + l1 = l+1 + go to 50 + 60 call fpbspl(ty,ny,ky,arg,l,h) + ly(i) = l-ky1 + do 70 j=1,ky1 + wy(i,j) = h(j) + 70 continue + 80 continue + m = 0 + do 130 i=1,mx + l = lx(i)*nky1 + do 90 i1=1,kx1 + h(i1) = wx(i,i1) + 90 continue + do 120 j=1,my + l1 = l+ly(j) + sp = 0. + do 110 i1=1,kx1 + l2 = l1 + do 100 j1=1,ky1 + l2 = l2+1 + sp = sp+c(l2)*h(i1)*wy(j,j1) + 100 continue + l1 = l1+nky1 + 110 continue + m = m+1 + z(m) = sp + 120 continue + 130 continue + return + end + + subroutine fpbspl(t,n,k,x,l,h) +c subroutine fpbspl evaluates the (k+1) non-zero b-splines of +c degree k at t(l) <= x < t(l+1) using the stable recurrence +c relation of de boor and cox. +c .. +c ..scalar arguments.. + real*8 x + integer n,k,l +c ..array arguments.. + real*8 t(n),h(6) +c ..local scalars.. + real*8 f,one + integer i,j,li,lj +c ..local arrays.. + real*8 hh(5) +c .. + one = 0.1d+01 + h(1) = one + do 20 j=1,k + do 10 i=1,j + hh(i) = h(i) + 10 continue + h(1) = 0.0d0 + do 20 i=1,j + li = l+i + lj = li-j + f = hh(i)/(t(li)-t(lj)) + h(i) = h(i)+f*(t(li)-x) + h(i+1) = f*(x-t(lj)) + 20 continue + return + end + subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z, + * wrk,lwrk,iwrk,kwrk,ier) +c subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,... +c ,my the partial derivative ( order nux,nuy) of a bivariate spline +c s(x,y) of degrees kx and ky, given in the b-spline representation. +c +c calling sequence: +c call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk, +c * iwrk,kwrk,ier) +c +c input parameters: +c tx : real array, length nx, which contains the position of the +c knots in the x-direction. +c nx : integer, giving the total number of knots in the x-direction +c ty : real array, length ny, which contains the position of the +c knots in the y-direction. +c ny : integer, giving the total number of knots in the y-direction +c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the +c b-spline coefficients. +c kx,ky : integer values, giving the degrees of the spline. +c nux : integer values, specifying the order of the partial +c nuy derivative. 0<=nux=1. +c y : real array of dimension (my). +c before entry y(j) must be set to the y co-ordinate of the +c j-th grid point along the y-axis. +c ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my. +c my : on entry my must specify the number of grid points along +c the y-axis. my >=1. +c wrk : real array of dimension lwrk. used as workspace. +c lwrk : integer, specifying the dimension of wrk. +c lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1) +c iwrk : integer array of dimension kwrk. used as workspace. +c kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my. +c +c output parameters: +c z : real array of dimension (mx*my). +c on succesful exit z(my*(i-1)+j) contains the value of the +c specified partial derivative of s(x,y) at the point +c (x(i),y(j)),i=1,...,mx;j=1,...,my. +c ier : integer error flag +c ier=0 : normal return +c ier=10: invalid input data (see restrictions) +c +c restrictions: +c mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my +c lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1), +c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx +c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my +c +c other subroutines required: +c fpbisp,fpbspl +c +c references : +c de boor c : on calculating with b-splines, j. approximation theory +c 6 (1972) 50-62. +c dierckx p. : curve and surface fitting with splines, monographs on +c numerical analysis, oxford university press, 1993. +c +c author : +c p.dierckx +c dept. computer science, k.u.leuven +c celestijnenlaan 200a, b-3001 heverlee, belgium. +c e-mail : Paul.Dierckx@cs.kuleuven.ac.be +c +c latest update : march 1989 +c +c ..scalar arguments.. + integer nx,ny,kx,ky,nux,nuy,mx,my,lwrk,kwrk,ier +c ..array arguments.. + integer iwrk(kwrk) + real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my), + * wrk(lwrk) +c ..local scalars.. + integer i,iwx,iwy,j,kkx,kky,kx1,ky1,lx,ly,lwest,l1,l2,m,m0,m1, + * nc,nkx1,nky1,nxx,nyy + real*8 ak,fac +c .. +c before starting computations a data check is made. if the input data +c are invalid control is immediately repassed to the calling program. + ier = 10 + kx1 = kx+1 + ky1 = ky+1 + nkx1 = nx-kx1 + nky1 = ny-ky1 + nc = nkx1*nky1 + if(nux.lt.0 .or. nux.ge.kx) go to 400 + if(nuy.lt.0 .or. nuy.ge.ky) go to 400 + lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my + if(lwrk.lt.lwest) go to 400 + if(kwrk.lt.(mx+my)) go to 400 + if(mx-1) 400,30,10 + 10 do 20 i=2,mx + if(x(i).lt.x(i-1)) go to 400 + 20 continue + 30 if(my-1) 400,60,40 + 40 do 50 i=2,my + if(y(i).lt.y(i-1)) go to 400 + 50 continue + 60 ier = 0 + nxx = nkx1 + nyy = nky1 + kkx = kx + kky = ky +c the partial derivative of order (nux,nuy) of a bivariate spline of +c degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy. +c we calculate the b-spline coefficients of this spline + do 70 i=1,nc + wrk(i) = c(i) + 70 continue + if(nux.eq.0) go to 200 + lx = 1 + do 100 j=1,nux + ak = kkx + nxx = nxx-1 + l1 = lx + m0 = 1 + do 90 i=1,nxx + l1 = l1+1 + l2 = l1+kkx + fac = tx(l2)-tx(l1) + if(fac.le.0.) go to 90 + do 80 m=1,nyy + m1 = m0+nyy + wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac + m0 = m0+1 + 80 continue + 90 continue + lx = lx+1 + kkx = kkx-1 + 100 continue + 200 if(nuy.eq.0) go to 300 + ly = 1 + do 230 j=1,nuy + ak = kky + nyy = nyy-1 + l1 = ly + do 220 i=1,nyy + l1 = l1+1 + l2 = l1+kky + fac = ty(l2)-ty(l1) + if(fac.le.0.) go to 220 + m0 = i + do 210 m=1,nxx + m1 = m0+1 + wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac + m0 = m0+nky1 + 210 continue + 220 continue + ly = ly+1 + kky = kky-1 + 230 continue + m0 = nyy + m1 = nky1 + do 250 m=2,nxx + do 240 i=1,nyy + m0 = m0+1 + m1 = m1+1 + wrk(m0) = wrk(m1) + 240 continue + m1 = m1+nuy + 250 continue +c we partition the working space and evaluate the partial derivative + 300 iwx = 1+nxx*nyy + iwy = iwx+mx*(kx1-nux) + call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky, + * x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1)) + 400 return + end +c + subroutine surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, + * nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier) +c given the set of data points (x(i),y(i),z(i)) and the set of positive +c numbers w(i),i=1,...,m, subroutine surfit determines a smooth bivar- +c iate spline approximation s(x,y) of degrees kx and ky on the rect- +c angle xb <= x <= xe, yb <= y <= ye. +c if iopt = -1 surfit calculates the weighted least-squares spline +c according to a given set of knots. +c if iopt >= 0 the total numbers nx and ny of these knots and their +c position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic- +c ally by the routine. the smoothness of s(x,y) is then achieved by +c minimalizing the discontinuity jumps in the derivatives of s(x,y) +c across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1). +c the amounth of smoothness is determined by the condition that f(p) = +c sum ((w(i)*(z(i)-s(x(i),y(i))))**2) be <= s, with s a given non-neg- +c ative constant, called the smoothing factor. +c the fit is given in the b-spline representation (b-spline coefficients +c c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval- +c uated by means of subroutine bispev. +c +c calling sequence: +c call surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, +c * nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier) +c +c parameters: +c iopt : integer flag. on entry iopt must specify whether a weighted +c least-squares spline (iopt=-1) or a smoothing spline (iopt=0 +c or 1) must be determined. +c if iopt=0 the routine will start with an initial set of knots +c tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i= +c 1,...,ky+1. if iopt=1 the routine will continue with the set +c of knots found at the last call of the routine. +c attention: a call with iopt=1 must always be immediately pre- +c ceded by another call with iopt=1 or iopt=0. +c unchanged on exit. +c m : integer. on entry m must specify the number of data points. +c m >= (kx+1)*(ky+1). unchanged on exit. +c x : real array of dimension at least (m). +c y : real array of dimension at least (m). +c z : real array of dimension at least (m). +c before entry, x(i),y(i),z(i) must be set to the co-ordinates +c of the i-th data point, for i=1,...,m. the order of the data +c points is immaterial. unchanged on exit. +c w : real array of dimension at least (m). before entry, w(i) must +c be set to the i-th value in the set of weights. the w(i) must +c be strictly positive. unchanged on exit. +c xb,xe : real values. on entry xb,xe,yb and ye must specify the bound- +c yb,ye aries of the rectangular approximation domain. +c xb<=x(i)<=xe,yb<=y(i)<=ye,i=1,...,m. unchanged on exit. +c kx,ky : integer values. on entry kx and ky must specify the degrees +c of the spline. 1<=kx,ky<=5. it is recommended to use bicubic +c (kx=ky=3) splines. unchanged on exit. +c s : real. on entry (in case iopt>=0) s must specify the smoothing +c factor. s >=0. unchanged on exit. +c for advice on the choice of s see further comments +c nxest : integer. unchanged on exit. +c nyest : integer. unchanged on exit. +c on entry, nxest and nyest must specify an upper bound for the +c number of knots required in the x- and y-directions respect. +c these numbers will also determine the storage space needed by +c the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1). +c in most practical situation nxest = kx+1+sqrt(m/2), nyest = +c ky+1+sqrt(m/2) will be sufficient. see also further comments. +c nmax : integer. on entry nmax must specify the actual dimension of +c the arrays tx and ty. nmax >= nxest, nmax >=nyest. +c unchanged on exit. +c eps : real. +c on entry, eps must specify a threshold for determining the +c effective rank of an over-determined linear system of equat- +c ions. 0 < eps < 1. if the number of decimal digits in the +c computer representation of a real number is q, then 10**(-q) +c is a suitable value for eps in most practical applications. +c unchanged on exit. +c nx : integer. +c unless ier=10 (in case iopt >=0), nx will contain the total +c number of knots with respect to the x-variable, of the spline +c approximation returned. if the computation mode iopt=1 is +c used, the value of nx should be left unchanged between sub- +c sequent calls. +c in case iopt=-1, the value of nx should be specified on entry +c tx : real array of dimension nmax. +c on succesful exit, this array will contain the knots of the +c spline with respect to the x-variable, i.e. the position of +c the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the +c position of the additional knots tx(1)=...=tx(kx+1)=xb and +c tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat. +c if the computation mode iopt=1 is used, the values of tx(1), +c ...,tx(nx) should be left unchanged between subsequent calls. +c if the computation mode iopt=-1 is used, the values tx(kx+2), +c ...tx(nx-kx-1) must be supplied by the user, before entry. +c see also the restrictions (ier=10). +c ny : integer. +c unless ier=10 (in case iopt >=0), ny will contain the total +c number of knots with respect to the y-variable, of the spline +c approximation returned. if the computation mode iopt=1 is +c used, the value of ny should be left unchanged between sub- +c sequent calls. +c in case iopt=-1, the value of ny should be specified on entry +c ty : real array of dimension nmax. +c on succesful exit, this array will contain the knots of the +c spline with respect to the y-variable, i.e. the position of +c the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the +c position of the additional knots ty(1)=...=ty(ky+1)=yb and +c ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat. +c if the computation mode iopt=1 is used, the values of ty(1), +c ...,ty(ny) should be left unchanged between subsequent calls. +c if the computation mode iopt=-1 is used, the values ty(ky+2), +c ...ty(ny-ky-1) must be supplied by the user, before entry. +c see also the restrictions (ier=10). +c c : real array of dimension at least (nxest-kx-1)*(nyest-ky-1). +c on succesful exit, c contains the coefficients of the spline +c approximation s(x,y) +c fp : real. unless ier=10, fp contains the weighted sum of +c squared residuals of the spline approximation returned. +c wrk1 : real array of dimension (lwrk1). used as workspace. +c if the computation mode iopt=1 is used the value of wrk1(1) +c should be left unchanged between subsequent calls. +c on exit wrk1(2),wrk1(3),...,wrk1(1+(nx-kx-1)*(ny-ky-1)) will +c contain the values d(i)/max(d(i)),i=1,...,(nx-kx-1)*(ny-ky-1) +c with d(i) the i-th diagonal element of the reduced triangular +c matrix for calculating the b-spline coefficients. it includes +c those elements whose square is less than eps,which are treat- +c ed as 0 in the case of presumed rank deficiency (ier<-2). +c lwrk1 : integer. on entry lwrk1 must specify the actual dimension of +c the array wrk1 as declared in the calling (sub)program. +c lwrk1 must not be too small. let +c u = nxest-kx-1, v = nyest-ky-1, km = max(kx,ky)+1, +c ne = max(nxest,nyest), bx = kx*v+ky+1, by = ky*u+kx+1, +c if(bx.le.by) b1 = bx, b2 = b1+v-ky +c if(bx.gt.by) b1 = by, b2 = b1+u-kx then +c lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1 +c wrk2 : real array of dimension (lwrk2). used as workspace, but +c only in the case a rank deficient system is encountered. +c lwrk2 : integer. on entry lwrk2 must specify the actual dimension of +c the array wrk2 as declared in the calling (sub)program. +c lwrk2 > 0 . a save upper boundfor lwrk2 = u*v*(b2+1)+b2 +c where u,v and b2 are as above. if there are enough data +c points, scattered uniformly over the approximation domain +c and if the smoothing factor s is not too small, there is a +c good chance that this extra workspace is not needed. a lot +c of memory might therefore be saved by setting lwrk2=1. +c (see also ier > 10) +c iwrk : integer array of dimension (kwrk). used as workspace. +c kwrk : integer. on entry kwrk must specify the actual dimension of +c the array iwrk as declared in the calling (sub)program. +c kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1). +c ier : integer. unless the routine detects an error, ier contains a +c non-positive value on exit, i.e. +c ier=0 : normal return. the spline returned has a residual sum of +c squares fp such that abs(fp-s)/s <= tol with tol a relat- +c ive tolerance set to 0.001 by the program. +c ier=-1 : normal return. the spline returned is an interpolating +c spline (fp=0). +c ier=-2 : normal return. the spline returned is the weighted least- +c squares polynomial of degrees kx and ky. in this extreme +c case fp gives the upper bound for the smoothing factor s. +c ier<-2 : warning. the coefficients of the spline returned have been +c computed as the minimal norm least-squares solution of a +c (numerically) rank deficient system. (-ier) gives the rank. +c especially if the rank deficiency which can be computed as +c (nx-kx-1)*(ny-ky-1)+ier, is large the results may be inac- +c curate. they could also seriously depend on the value of +c eps. +c ier=1 : error. the required storage space exceeds the available +c storage space, as specified by the parameters nxest and +c nyest. +c probably causes : nxest or nyest too small. if these param- +c eters are already large, it may also indicate that s is +c too small +c the approximation returned is the weighted least-squares +c spline according to the current set of knots. +c the parameter fp gives the corresponding weighted sum of +c squared residuals (fp>s). +c ier=2 : error. a theoretically impossible result was found during +c the iteration proces for finding a smoothing spline with +c fp = s. probably causes : s too small or badly chosen eps. +c there is an approximation returned but the corresponding +c weighted sum of squared residuals does not satisfy the +c condition abs(fp-s)/s < tol. +c ier=3 : error. the maximal number of iterations maxit (set to 20 +c by the program) allowed for finding a smoothing spline +c with fp=s has been reached. probably causes : s too small +c there is an approximation returned but the corresponding +c weighted sum of squared residuals does not satisfy the +c condition abs(fp-s)/s < tol. +c ier=4 : error. no more knots can be added because the number of +c b-spline coefficients (nx-kx-1)*(ny-ky-1) already exceeds +c the number of data points m. +c probably causes : either s or m too small. +c the approximation returned is the weighted least-squares +c spline according to the current set of knots. +c the parameter fp gives the corresponding weighted sum of +c squared residuals (fp>s). +c ier=5 : error. no more knots can be added because the additional +c knot would (quasi) coincide with an old one. +c probably causes : s too small or too large a weight to an +c inaccurate data point. +c the approximation returned is the weighted least-squares +c spline according to the current set of knots. +c the parameter fp gives the corresponding weighted sum of +c squared residuals (fp>s). +c ier=10 : error. on entry, the input data are controlled on validity +c the following restrictions must be satisfied. +c -1<=iopt<=1, 1<=kx,ky<=5, m>=(kx+1)*(ky+1), nxest>=2*kx+2, +c nyest>=2*ky+2, 0=nxest, nmax>=nyest, +c xb<=x(i)<=xe, yb<=y(i)<=ye, w(i)>0, i=1,...,m +c lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1 +c kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1) +c if iopt=-1: 2*kx+2<=nx<=nxest +c xb=0: s>=0 +c if one of these conditions is found to be violated,control +c is immediately repassed to the calling program. in that +c case there is no approximation returned. +c ier>10 : error. lwrk2 is too small, i.e. there is not enough work- +c space for computing the minimal least-squares solution of +c a rank deficient system of linear equations. ier gives the +c requested value for lwrk2. there is no approximation re- +c turned but, having saved the information contained in nx, +c ny,tx,ty,wrk1, and having adjusted the value of lwrk2 and +c the dimension of the array wrk2 accordingly, the user can +c continue at the point the program was left, by calling +c surfit with iopt=1. +c +c further comments: +c by means of the parameter s, the user can control the tradeoff +c between closeness of fit and smoothness of fit of the approximation. +c if s is too large, the spline will be too smooth and signal will be +c lost ; if s is too small the spline will pick up too much noise. in +c the extreme cases the program will return an interpolating spline if +c s=0 and the weighted least-squares polynomial (degrees kx,ky)if s is +c very large. between these extremes, a properly chosen s will result +c in a good compromise between closeness of fit and smoothness of fit. +c to decide whether an approximation, corresponding to a certain s is +c satisfactory the user is highly recommended to inspect the fits +c graphically. +c recommended values for s depend on the weights w(i). if these are +c taken as 1/d(i) with d(i) an estimate of the standard deviation of +c z(i), a good s-value should be found in the range (m-sqrt(2*m),m+ +c sqrt(2*m)). if nothing is known about the statistical error in z(i) +c each w(i) can be set equal to one and s determined by trial and +c error, taking account of the comments above. the best is then to +c start with a very large value of s ( to determine the least-squares +c polynomial and the corresponding upper bound fp0 for s) and then to +c progressively decrease the value of s ( say by a factor 10 in the +c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the +c approximation shows more detail) to obtain closer fits. +c to choose s very small is strongly discouraged. this considerably +c increases computation time and memory requirements. it may also +c cause rank-deficiency (ier<-2) and endager numerical stability. +c to economize the search for a good s-value the program provides with +c different modes of computation. at the first call of the routine, or +c whenever he wants to restart with the initial set of knots the user +c must set iopt=0. +c if iopt=1 the program will continue with the set of knots found at +c the last call of the routine. this will save a lot of computation +c time if surfit is called repeatedly for different values of s. +c the number of knots of the spline returned and their location will +c depend on the value of s and on the complexity of the shape of the +c function underlying the data. if the computation mode iopt=1 +c is used, the knots returned may also depend on the s-values at +c previous calls (if these were smaller). therefore, if after a number +c of trials with different s-values and iopt=1, the user can finally +c accept a fit as satisfactory, it may be worthwhile for him to call +c surfit once more with the selected value for s but now with iopt=0. +c indeed, surfit may then return an approximation of the same quality +c of fit but with fewer knots and therefore better if data reduction +c is also an important objective for the user. +c the number of knots may also depend on the upper bounds nxest and +c nyest. indeed, if at a certain stage in surfit the number of knots +c in one direction (say nx) has reached the value of its upper bound +c (nxest), then from that moment on all subsequent knots are added +c in the other (y) direction. this may indicate that the value of +c nxest is too small. on the other hand, it gives the user the option +c of limiting the number of knots the routine locates in any direction +c for example, by setting nxest=2*kx+2 (the lowest allowable value for +c nxest), the user can indicate that he wants an approximation which +c is a simple polynomial of degree kx in the variable x. +c +c other subroutines required: +c fpback,fpbspl,fpsurf,fpdisc,fpgivs,fprank,fprati,fprota,fporde +c +c references: +c dierckx p. : an algorithm for surface fitting with spline functions +c ima j. numer. anal. 1 (1981) 267-283. +c dierckx p. : an algorithm for surface fitting with spline functions +c report tw50, dept. computer science,k.u.leuven, 1980. +c dierckx p. : curve and surface fitting with splines, monographs on +c numerical analysis, oxford university press, 1993. +c +c author: +c p.dierckx +c dept. computer science, k.u. leuven +c celestijnenlaan 200a, b-3001 heverlee, belgium. +c e-mail : Paul.Dierckx@cs.kuleuven.ac.be +c +c creation date : may 1979 +c latest update : march 1987 +c +c .. +c ..scalar arguments.. + real*8 xb,xe,yb,ye,s,eps,fp + integer iopt,m,kx,ky,nxest,nyest,nmax,nx,ny,lwrk1,lwrk2,kwrk,ier +c ..array arguments.. + real*8 x(m),y(m),z(m),w(m),tx(nmax),ty(nmax), + * c((nxest-kx-1)*(nyest-ky-1)),wrk1(lwrk1),wrk2(lwrk2) + integer iwrk(kwrk) +c ..local scalars.. + real*8 tol + integer i,ib1,ib3,jb1,ki,kmax,km1,km2,kn,kwest,kx1,ky1,la,lbx, + * lby,lco,lf,lff,lfp,lh,lq,lsx,lsy,lwest,maxit,ncest,nest,nek, + * nminx,nminy,nmx,nmy,nreg,nrint,nxk,nyk +c ..function references.. + integer max0 +c ..subroutine references.. +c fpsurf +c .. +c we set up the parameters tol and maxit. + maxit = 20 + tol = 0.1e-02 +c before starting computations a data check is made. if the input data +c are invalid,control is immediately repassed to the calling program. + ier = 10 + if(eps.le.0. .or. eps.ge.1.) go to 70 + if(kx.le.0 .or. kx.gt.5) go to 70 + kx1 = kx+1 + if(ky.le.0 .or. ky.gt.5) go to 70 + ky1 = ky+1 + kmax = max0(kx,ky) + km1 = kmax+1 + km2 = km1+1 + if(iopt.lt.(-1) .or. iopt.gt.1) go to 70 + if(m.lt.(kx1*ky1)) go to 70 + nminx = 2*kx1 + if(nxest.lt.nminx .or. nxest.gt.nmax) go to 70 + nminy = 2*ky1 + if(nyest.lt.nminy .or. nyest.gt.nmax) go to 70 + nest = max0(nxest,nyest) + nxk = nxest-kx1 + nyk = nyest-ky1 + ncest = nxk*nyk + nmx = nxest-nminx+1 + nmy = nyest-nminy+1 + nrint = nmx+nmy + nreg = nmx*nmy + ib1 = kx*nyk+ky1 + jb1 = ky*nxk+kx1 + ib3 = kx1*nyk+1 + if(ib1.le.jb1) go to 10 + ib1 = jb1 + ib3 = ky1*nxk+1 + 10 lwest = ncest*(2+ib1+ib3)+2*(nrint+nest*km2+m*km1)+ib3 + kwest = m+nreg + if(lwrk1.lt.lwest .or. kwrk.lt.kwest) go to 70 + if(xb.ge.xe .or. yb.ge.ye) go to 70 + do 20 i=1,m + if(w(i).le.0.) go to 70 + if(x(i).lt.xb .or. x(i).gt.xe) go to 70 + if(y(i).lt.yb .or. y(i).gt.ye) go to 70 + 20 continue + if(iopt.ge.0) go to 50 + if(nx.lt.nminx .or. nx.gt.nxest) go to 70 + nxk = nx-kx1 + tx(kx1) = xb + tx(nxk+1) = xe + do 30 i=kx1,nxk + if(tx(i+1).le.tx(i)) go to 70 + 30 continue + if(ny.lt.nminy .or. ny.gt.nyest) go to 70 + nyk = ny-ky1 + ty(ky1) = yb + ty(nyk+1) = ye + do 40 i=ky1,nyk + if(ty(i+1).le.ty(i)) go to 70 + 40 continue + go to 60 + 50 if(s.lt.0.) go to 70 + 60 ier = 0 +c we partition the working space and determine the spline approximation + kn = 1 + ki = kn+m + lq = 2 + la = lq+ncest*ib3 + lf = la+ncest*ib1 + lff = lf+ncest + lfp = lff+ncest + lco = lfp+nrint + lh = lco+nrint + lbx = lh+ib3 + nek = nest*km2 + lby = lbx+nek + lsx = lby+nek + lsy = lsx+m*km1 + call fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, + * eps,tol,maxit,nest,km1,km2,ib1,ib3,ncest,nrint,nreg,nx,tx, + * ny,ty,c,fp,wrk1(1),wrk1(lfp),wrk1(lco),wrk1(lf),wrk1(lff), + * wrk1(la),wrk1(lq),wrk1(lbx),wrk1(lby),wrk1(lsx),wrk1(lsy), + * wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier) + 70 return + end + subroutine fpback(a,z,n,k,c,nest) +c subroutine fpback calculates the solution of the system of +c equations a*c = z with a a n x n upper triangular matrix +c of bandwidth k. +c .. +c ..scalar arguments.. + integer n,k,nest +c ..array arguments.. + real*8 a(nest,k),z(n),c(n) +c ..local scalars.. + real*8 store + integer i,i1,j,k1,l,m +c .. + k1 = k-1 + c(n) = z(n)/a(n,1) + i = n-1 + if(i.eq.0) go to 30 + do 20 j=2,n + store = z(i) + i1 = k1 + if(j.le.k1) i1 = j-1 + m = i + do 10 l=1,i1 + m = m+1 + store = store-c(m)*a(i,l+1) + 10 continue + c(i) = store/a(i,1) + i = i-1 + 20 continue + 30 return + end + subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest, + * nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest, + * nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h, + * index,nummer,wrk,lwrk,ier) +c .. +c ..scalar arguments.. + real*8 xb,xe,yb,ye,s,eta,tol,fp,fp0 + integer iopt,m,kxx,kyy,nxest,nyest,maxit,nmax,km1,km2,ib1,ib3, + * nc,intest,nrest,nx0,ny0,lwrk,ier +c ..array arguments.. + real*8 x(m),y(m),z(m),w(m),tx(nmax),ty(nmax),c(nc),fpint(intest), + * coord(intest),f(nc),ff(nc),a(nc,ib1),q(nc,ib3),bx(nmax,km2), + * by(nmax,km2),spx(m,km1),spy(m,km1),h(ib3),wrk(lwrk) + integer index(nrest),nummer(m) +c ..local scalars.. + real*8 acc,arg,cos,dmax,fac1,fac2,fpmax,fpms,f1,f2,f3,hxi,p,pinv, + * piv,p1,p2,p3,sigma,sin,sq,store,wi,x0,x1,y0,y1,zi,eps, + * rn,one,con1,con9,con4,half,ten + integer i,iband,iband1,iband3,iband4,ibb,ichang,ich1,ich3,ii, + * in,irot,iter,i1,i2,i3,j,jrot,jxy,j1,kx,kx1,kx2,ky,ky1,ky2,l, + * la,lf,lh,lwest,lx,ly,l1,l2,n,ncof,nk1x,nk1y,nminx,nminy,nreg, + * nrint,num,num1,nx,nxe,nxx,ny,nye,nyy,n1,rank +c ..local arrays.. + real*8 hx(6),hy(6) +c ..function references.. + real*8 abs,fprati,sqrt + integer min0 +c ..subroutine references.. +c fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota +c .. +c set constants + one = 0.1e+01 + con1 = 0.1e0 + con9 = 0.9e0 + con4 = 0.4e-01 + half = 0.5e0 + ten = 0.1e+02 +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c part 1: determination of the number of knots and their position. c +c **************************************************************** c +c given a set of knots we compute the least-squares spline sinf(x,y), c +c and the corresponding weighted sum of squared residuals fp=f(p=inf). c +c if iopt=-1 sinf(x,y) is the requested approximation. c +c if iopt=0 or iopt=1 we check whether we can accept the knots: c +c if fp <=s we will continue with the current set of knots. c +c if fp > s we will increase the number of knots and compute the c +c corresponding least-squares spline until finally fp<=s. c +c the initial choice of knots depends on the value of s and iopt. c +c if iopt=0 we first compute the least-squares polynomial of degree c +c kx in x and ky in y; nx=nminx=2*kx+2 and ny=nminy=2*ky+2. c +c fp0=f(0) denotes the corresponding weighted sum of squared c +c residuals c +c if iopt=1 we start with the knots found at the last call of the c +c routine, except for the case that s>=fp0; then we can compute c +c the least-squares polynomial directly. c +c eventually the independent variables x and y (and the corresponding c +c parameters) will be switched if this can reduce the bandwidth of the c +c system to be solved. c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c ichang denotes whether(1) or not(-1) the directions have been inter- +c changed. + ichang = -1 + x0 = xb + x1 = xe + y0 = yb + y1 = ye + kx = kxx + ky = kyy + kx1 = kx+1 + ky1 = ky+1 + nxe = nxest + nye = nyest + eps = sqrt(eta) + if(iopt.lt.0) go to 20 +c calculation of acc, the absolute tolerance for the root of f(p)=s. + acc = tol*s + if(iopt.eq.0) go to 10 + if(fp0.gt.s) go to 20 +c initialization for the least-squares polynomial. + 10 nminx = 2*kx1 + nminy = 2*ky1 + nx = nminx + ny = nminy + ier = -2 + go to 30 + 20 nx = nx0 + ny = ny0 +c main loop for the different sets of knots. m is a save upper bound +c for the number of trials. + 30 do 420 iter=1,m +c find the position of the additional knots which are needed for the +c b-spline representation of s(x,y). + l = nx + do 40 i=1,kx1 + tx(i) = x0 + tx(l) = x1 + l = l-1 + 40 continue + l = ny + do 50 i=1,ky1 + ty(i) = y0 + ty(l) = y1 + l = l-1 + 50 continue +c find nrint, the total number of knot intervals and nreg, the number +c of panels in which the approximation domain is subdivided by the +c intersection of knots. + nxx = nx-2*kx1+1 + nyy = ny-2*ky1+1 + nrint = nxx+nyy + nreg = nxx*nyy +c find the bandwidth of the observation matrix a. +c if necessary, interchange the variables x and y, in order to obtain +c a minimal bandwidth. + iband1 = kx*(ny-ky1)+ky + l = ky*(nx-kx1)+kx + if(iband1.le.l) go to 130 + iband1 = l + ichang = -ichang + do 60 i=1,m + store = x(i) + x(i) = y(i) + y(i) = store + 60 continue + store = x0 + x0 = y0 + y0 = store + store = x1 + x1 = y1 + y1 = store + n = min0(nx,ny) + do 70 i=1,n + store = tx(i) + tx(i) = ty(i) + ty(i) = store + 70 continue + n1 = n+1 + if(nx-ny) 80,120,100 + 80 do 90 i=n1,ny + tx(i) = ty(i) + 90 continue + go to 120 + 100 do 110 i=n1,nx + ty(i) = tx(i) + 110 continue + 120 l = nx + nx = ny + ny = l + l = nxe + nxe = nye + nye = l + l = nxx + nxx = nyy + nyy = l + l = kx + kx = ky + ky = l + kx1 = kx+1 + ky1 = ky+1 + 130 iband = iband1+1 +c arrange the data points according to the panel they belong to. + call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg) +c find ncof, the number of b-spline coefficients. + nk1x = nx-kx1 + nk1y = ny-ky1 + ncof = nk1x*nk1y +c initialize the observation matrix a. + do 140 i=1,ncof + f(i) = 0. + do 140 j=1,iband + a(i,j) = 0. + 140 continue +c initialize the sum of squared residuals. + fp = 0. +c fetch the data points in the new order. main loop for the +c different panels. + do 250 num=1,nreg +c fix certain constants for the current panel; jrot records the column +c number of the first non-zero element in a row of the observation +c matrix according to a data point of the panel. + num1 = num-1 + lx = num1/nyy + l1 = lx+kx1 + ly = num1-lx*nyy + l2 = ly+ky1 + jrot = lx*nk1y+ly +c test whether there are still data points in the panel. + in = index(num) + 150 if(in.eq.0) go to 250 +c fetch a new data point. + wi = w(in) + zi = z(in)*wi +c evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in). + call fpbspl(tx,nx,kx,x(in),l1,hx) +c evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in). + call fpbspl(ty,ny,ky,y(in),l2,hy) +c store the value of these b-splines in spx and spy respectively. + do 160 i=1,kx1 + spx(in,i) = hx(i) + 160 continue + do 170 i=1,ky1 + spy(in,i) = hy(i) + 170 continue +c initialize the new row of observation matrix. + do 180 i=1,iband + h(i) = 0. + 180 continue +c calculate the non-zero elements of the new row by making the cross +c products of the non-zero b-splines in x- and y-direction. + i1 = 0 + do 200 i=1,kx1 + hxi = hx(i) + j1 = i1 + do 190 j=1,ky1 + j1 = j1+1 + h(j1) = hxi*hy(j)*wi + 190 continue + i1 = i1+nk1y + 200 continue +c rotate the row into triangle by givens transformations . + irot = jrot + do 220 i=1,iband + irot = irot+1 + piv = h(i) + if(piv.eq.0.) go to 220 +c calculate the parameters of the givens transformation. + call fpgivs(piv,a(irot,1),cos,sin) +c apply that transformation to the right hand side. + call fprota(cos,sin,zi,f(irot)) + if(i.eq.iband) go to 230 +c apply that transformation to the left hand side. + i2 = 1 + i3 = i+1 + do 210 j=i3,iband + i2 = i2+1 + call fprota(cos,sin,h(j),a(irot,i2)) + 210 continue + 220 continue +c add the contribution of the row to the sum of squares of residual +c right hand sides. + 230 fp = fp+zi**2 +c find the number of the next data point in the panel. + 240 in = nummer(in) + go to 150 + 250 continue +c find dmax, the maximum value for the diagonal elements in the reduced +c triangle. + dmax = 0. + do 260 i=1,ncof + if(a(i,1).le.dmax) go to 260 + dmax = a(i,1) + 260 continue +c check whether the observation matrix is rank deficient. + sigma = eps*dmax + do 270 i=1,ncof + if(a(i,1).le.sigma) go to 280 + 270 continue +c backward substitution in case of full rank. + call fpback(a,f,ncof,iband,c,nc) + rank = ncof + do 275 i=1,ncof + q(i,1) = a(i,1)/dmax + 275 continue + go to 300 +c in case of rank deficiency, find the minimum norm solution. +c check whether there is sufficient working space + 280 lwest = ncof*iband+ncof+iband + if(lwrk.lt.lwest) go to 780 + do 290 i=1,ncof + ff(i) = f(i) + do 290 j=1,iband + q(i,j) = a(i,j) + 290 continue + lf =1 + lh = lf+ncof + la = lh+iband + call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la), + * wrk(lf),wrk(lh)) + do 295 i=1,ncof + q(i,1) = q(i,1)/dmax + 295 continue +c add to the sum of squared residuals, the contribution of reducing +c the rank. + fp = fp+sq + 300 if(ier.eq.(-2)) fp0 = fp +c test whether the least-squares spline is an acceptable solution. + if(iopt.lt.0) go to 820 + fpms = fp-s + if(abs(fpms).le.acc) if(fp) 815,815,820 +c test whether we can accept the choice of knots. + if(fpms.lt.0.) go to 430 +c test whether we cannot further increase the number of knots. + if(ncof.gt.m) go to 790 + ier = 0 +c search where to add a new knot. +c find for each interval the sum of squared residuals fpint for the +c data points having the coordinate belonging to that knot interval. +c calculate also coord which is the same sum, weighted by the position +c of the data points considered. + 310 do 320 i=1,nrint + fpint(i) = 0. + coord(i) = 0. + 320 continue + do 360 num=1,nreg + num1 = num-1 + lx = num1/nyy + l1 = lx+1 + ly = num1-lx*nyy + l2 = ly+1+nxx + jrot = lx*nk1y+ly + in = index(num) + 330 if(in.eq.0) go to 360 + store = 0. + i1 = jrot + do 350 i=1,kx1 + hxi = spx(in,i) + j1 = i1 + do 340 j=1,ky1 + j1 = j1+1 + store = store+hxi*spy(in,j)*c(j1) + 340 continue + i1 = i1+nk1y + 350 continue + store = (w(in)*(z(in)-store))**2 + fpint(l1) = fpint(l1)+store + coord(l1) = coord(l1)+store*x(in) + fpint(l2) = fpint(l2)+store + coord(l2) = coord(l2)+store*y(in) + in = nummer(in) + go to 330 + 360 continue +c find the interval for which fpint is maximal on the condition that +c there still can be added a knot. + 370 l = 0 + fpmax = 0. + l1 = 1 + l2 = nrint + if(nx.eq.nxe) l1 = nxx+1 + if(ny.eq.nye) l2 = nxx + if(l1.gt.l2) go to 810 + do 380 i=l1,l2 + if(fpmax.ge.fpint(i)) go to 380 + l = i + fpmax = fpint(i) + 380 continue +c test whether we cannot further increase the number of knots. + if(l.eq.0) go to 785 +c calculate the position of the new knot. + arg = coord(l)/fpint(l) +c test in what direction the new knot is going to be added. + if(l.gt.nxx) go to 400 +c addition in the x-direction. + jxy = l+kx1 + fpint(l) = 0. + fac1 = tx(jxy)-arg + fac2 = arg-tx(jxy-1) + if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370 + j = nx + do 390 i=jxy,nx + tx(j+1) = tx(j) + j = j-1 + 390 continue + tx(jxy) = arg + nx = nx+1 + go to 420 +c addition in the y-direction. + 400 jxy = l+ky1-nxx + fpint(l) = 0. + fac1 = ty(jxy)-arg + fac2 = arg-ty(jxy-1) + if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370 + j = ny + do 410 i=jxy,ny + ty(j+1) = ty(j) + j = j-1 + 410 continue + ty(jxy) = arg + ny = ny+1 +c restart the computations with the new set of knots. + 420 continue +c test whether the least-squares polynomial is a solution of our +c approximation problem. + 430 if(ier.eq.(-2)) go to 830 +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c part 2: determination of the smoothing spline sp(x,y) c +c ***************************************************** c +c we have determined the number of knots and their position. we now c +c compute the b-spline coefficients of the smoothing spline sp(x,y). c +c the observation matrix a is extended by the rows of a matrix, c +c expressing that sp(x,y) must be a polynomial of degree kx in x and c +c ky in y. the corresponding weights of these additional rows are set c +c to 1./p. iteratively we than have to determine the value of p c +c such that f(p)=sum((w(i)*(z(i)-sp(x(i),y(i))))**2) be = s. c +c we already know that the least-squares polynomial corresponds to c +c p=0 and that the least-squares spline corresponds to p=infinity. c +c the iteration process which is proposed here makes use of rational c +c interpolation. since f(p) is a convex and strictly decreasing c +c function of p, it can be approximated by a rational function r(p)= c +c (u*p+v)/(p+w). three values of p(p1,p2,p3) with corresponding values c +c of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the c +c new value of p such that r(p)=s. convergence is guaranteed by taking c +c f1 > 0 and f3 < 0. c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + kx2 = kx1+1 +c test whether there are interior knots in the x-direction. + if(nk1x.eq.kx1) go to 440 +c evaluate the discotinuity jumps of the kx-th order derivative of +c the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1. + call fpdisc(tx,nx,kx2,bx,nmax) + 440 ky2 = ky1 + 1 +c test whether there are interior knots in the y-direction. + if(nk1y.eq.ky1) go to 450 +c evaluate the discontinuity jumps of the ky-th order derivative of +c the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1. + call fpdisc(ty,ny,ky2,by,nmax) +c initial value for p. + 450 p1 = 0. + f1 = fp0-s + p3 = -one + f3 = fpms + p = 0. + do 460 i=1,ncof + p = p+a(i,1) + 460 continue + rn = ncof + p = rn/p +c find the bandwidth of the extended observation matrix. + iband3 = kx1*nk1y + iband4 = iband3 +1 + ich1 = 0 + ich3 = 0 +c iteration process to find the root of f(p)=s. + do 770 iter=1,maxit + pinv = one/p +c store the triangularized observation matrix into q. + do 480 i=1,ncof + ff(i) = f(i) + do 470 j=1,iband + q(i,j) = a(i,j) + 470 continue + ibb = iband+1 + do 480 j=ibb,iband4 + q(i,j) = 0. + 480 continue + if(nk1y.eq.ky1) go to 560 +c extend the observation matrix with the rows of a matrix, expressing +c that for x=cst. sp(x,y) must be a polynomial in y of degree ky. + do 550 i=ky2,nk1y + ii = i-ky1 + do 550 j=1,nk1x +c initialize the new row. + do 490 l=1,iband + h(l) = 0. + 490 continue +c fill in the non-zero elements of the row. jrot records the column +c number of the first non-zero element in the row. + do 500 l=1,ky2 + h(l) = by(ii,l)*pinv + 500 continue + zi = 0. + jrot = (j-1)*nk1y+ii +c rotate the new row into triangle by givens transformations without +c square roots. + do 540 irot=jrot,ncof + piv = h(1) + i2 = min0(iband1,ncof-irot) + if(piv.eq.0.) if(i2) 550,550,520 +c calculate the parameters of the givens transformation. + call fpgivs(piv,q(irot,1),cos,sin) +c apply that givens transformation to the right hand side. + call fprota(cos,sin,zi,ff(irot)) + if(i2.eq.0) go to 550 +c apply that givens transformation to the left hand side. + do 510 l=1,i2 + l1 = l+1 + call fprota(cos,sin,h(l1),q(irot,l1)) + 510 continue + 520 do 530 l=1,i2 + h(l) = h(l+1) + 530 continue + h(i2+1) = 0. + 540 continue + 550 continue + 560 if(nk1x.eq.kx1) go to 640 +c extend the observation matrix with the rows of a matrix expressing +c that for y=cst. sp(x,y) must be a polynomial in x of degree kx. + do 630 i=kx2,nk1x + ii = i-kx1 + do 630 j=1,nk1y +c initialize the new row + do 570 l=1,iband4 + h(l) = 0. + 570 continue +c fill in the non-zero elements of the row. jrot records the column +c number of the first non-zero element in the row. + j1 = 1 + do 580 l=1,kx2 + h(j1) = bx(ii,l)*pinv + j1 = j1+nk1y + 580 continue + zi = 0. + jrot = (i-kx2)*nk1y+j +c rotate the new row into triangle by givens transformations . + do 620 irot=jrot,ncof + piv = h(1) + i2 = min0(iband3,ncof-irot) + if(piv.eq.0.) if(i2) 630,630,600 +c calculate the parameters of the givens transformation. + call fpgivs(piv,q(irot,1),cos,sin) +c apply that givens transformation to the right hand side. + call fprota(cos,sin,zi,ff(irot)) + if(i2.eq.0) go to 630 +c apply that givens transformation to the left hand side. + do 590 l=1,i2 + l1 = l+1 + call fprota(cos,sin,h(l1),q(irot,l1)) + 590 continue + 600 do 610 l=1,i2 + h(l) = h(l+1) + 610 continue + h(i2+1) = 0. + 620 continue + 630 continue +c find dmax, the maximum value for the diagonal elements in the +c reduced triangle. + 640 dmax = 0. + do 650 i=1,ncof + if(q(i,1).le.dmax) go to 650 + dmax = q(i,1) + 650 continue +c check whether the matrix is rank deficient. + sigma = eps*dmax + do 660 i=1,ncof + if(q(i,1).le.sigma) go to 670 + 660 continue +c backward substitution in case of full rank. + call fpback(q,ff,ncof,iband4,c,nc) + rank = ncof + go to 675 +c in case of rank deficiency, find the minimum norm solution. + 670 lwest = ncof*iband4+ncof+iband4 + if(lwrk.lt.lwest) go to 780 + lf = 1 + lh = lf+ncof + la = lh+iband4 + call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la), + * wrk(lf),wrk(lh)) + 675 do 680 i=1,ncof + q(i,1) = q(i,1)/dmax + 680 continue +c compute f(p). + fp = 0. + do 720 num = 1,nreg + num1 = num-1 + lx = num1/nyy + ly = num1-lx*nyy + jrot = lx*nk1y+ly + in = index(num) + 690 if(in.eq.0) go to 720 + store = 0. + i1 = jrot + do 710 i=1,kx1 + hxi = spx(in,i) + j1 = i1 + do 700 j=1,ky1 + j1 = j1+1 + store = store+hxi*spy(in,j)*c(j1) + 700 continue + i1 = i1+nk1y + 710 continue + fp = fp+(w(in)*(z(in)-store))**2 + in = nummer(in) + go to 690 + 720 continue +c test whether the approximation sp(x,y) is an acceptable solution. + fpms = fp-s + if(abs(fpms).le.acc) go to 820 +c test whether the maximum allowable number of iterations has been +c reached. + if(iter.eq.maxit) go to 795 +c carry out one more step of the iteration process. + p2 = p + f2 = fpms + if(ich3.ne.0) go to 740 + if((f2-f3).gt.acc) go to 730 +c our initial choice of p is too large. + p3 = p2 + f3 = f2 + p = p*con4 + if(p.le.p1) p = p1*con9 + p2*con1 + go to 770 + 730 if(f2.lt.0.) ich3 = 1 + 740 if(ich1.ne.0) go to 760 + if((f1-f2).gt.acc) go to 750 +c our initial choice of p is too small + p1 = p2 + f1 = f2 + p = p/con4 + if(p3.lt.0.) go to 770 + if(p.ge.p3) p = p2*con1 + p3*con9 + go to 770 + 750 if(f2.gt.0.) ich1 = 1 +c test whether the iteration process proceeds as theoretically +c expected. + 760 if(f2.ge.f1 .or. f2.le.f3) go to 800 +c find the new value of p. + p = fprati(p1,f1,p2,f2,p3,f3) + 770 continue +c error codes and messages. + 780 ier = lwest + go to 830 + 785 ier = 5 + go to 830 + 790 ier = 4 + go to 830 + 795 ier = 3 + go to 830 + 800 ier = 2 + go to 830 + 810 ier = 1 + go to 830 + 815 ier = -1 + fp = 0. + 820 if(ncof.ne.rank) ier = -rank +c test whether x and y are in the original order. + 830 if(ichang.lt.0) go to 930 +c if not, interchange x and y once more. + l1 = 1 + do 840 i=1,nk1x + l2 = i + do 840 j=1,nk1y + f(l2) = c(l1) + l1 = l1+1 + l2 = l2+nk1x + 840 continue + do 850 i=1,ncof + c(i) = f(i) + 850 continue + do 860 i=1,m + store = x(i) + x(i) = y(i) + y(i) = store + 860 continue + n = min0(nx,ny) + do 870 i=1,n + store = tx(i) + tx(i) = ty(i) + ty(i) = store + 870 continue + n1 = n+1 + if(nx-ny) 880,920,900 + 880 do 890 i=n1,ny + tx(i) = ty(i) + 890 continue + go to 920 + 900 do 910 i=n1,nx + ty(i) = tx(i) + 910 continue + 920 l = nx + nx = ny + ny = l + 930 if(iopt.lt.0) go to 940 + nx0 = nx + ny0 = ny + 940 return + end + + subroutine fpdisc(t,n,k2,b,nest) +c subroutine fpdisc calculates the discontinuity jumps of the kth +c derivative of the b-splines of degree k at the knots t(k+2)..t(n-k-1) +c ..scalar arguments.. + integer n,k2,nest +c ..array arguments.. + real*8 t(n),b(nest,k2) +c ..local scalars.. + real*8 an,fac,prod + integer i,ik,j,jk,k,k1,l,lj,lk,lmk,lp,nk1,nrint +c ..local array.. + real*8 h(12) +c .. + k1 = k2-1 + k = k1-1 + nk1 = n-k1 + nrint = nk1-k + an = nrint + fac = an/(t(nk1+1)-t(k1)) + do 40 l=k2,nk1 + lmk = l-k1 + do 10 j=1,k1 + ik = j+k1 + lj = l+j + lk = lj-k2 + h(j) = t(l)-t(lk) + h(ik) = t(l)-t(lj) + 10 continue + lp = lmk + do 30 j=1,k2 + jk = j + prod = h(j) + do 20 i=1,k + jk = jk+1 + prod = prod*h(jk)*fac + 20 continue + lk = lp+k1 + b(lmk,j) = (t(lk)-t(lp))/prod + lp = lp+1 + 30 continue + 40 continue + return + end + subroutine fpgivs(piv,ww,cos,sin) +c subroutine fpgivs calculates the parameters of a givens +c transformation . +c .. +c ..scalar arguments.. + real*8 piv,ww,cos,sin +c ..local scalars.. + real*8 dd,one,store +c ..function references.. + real*8 abs,sqrt +c .. + one = 0.1e+01 + store = abs(piv) + if(store.ge.ww) dd = store*sqrt(one+(ww/piv)**2) + if(store.lt.ww) dd = ww*sqrt(one+(piv/ww)**2) + cos = ww/dd + sin = piv/dd + ww = dd + return + end + subroutine fprank(a,f,n,m,na,tol,c,sq,rank,aa,ff,h) +c subroutine fprank finds the minimum norm solution of a least- +c squares problem in case of rank deficiency. +c +c input parameters: +c a : array, which contains the non-zero elements of the observation +c matrix after triangularization by givens transformations. +c f : array, which contains the transformed right hand side. +c n : integer,wich contains the dimension of a. +c m : integer, which denotes the bandwidth of a. +c tol : real value, giving a threshold to determine the rank of a. +c +c output parameters: +c c : array, which contains the minimum norm solution. +c sq : real value, giving the contribution of reducing the rank +c to the sum of squared residuals. +c rank : integer, which contains the rank of matrix a. +c +c ..scalar arguments.. + integer n,m,na,rank + real*8 tol,sq +c ..array arguments.. + real*8 a(na,m),f(n),c(n),aa(n,m),ff(n),h(m) +c ..local scalars.. + integer i,ii,ij,i1,i2,j,jj,j1,j2,j3,k,kk,m1,nl + real*8 cos,fac,piv,sin,yi + double precision store,stor1,stor2,stor3 +c ..function references.. + integer min0 +c ..subroutine references.. +c fpgivs,fprota +c .. + m1 = m-1 +c the rank deficiency nl is considered to be the number of sufficient +c small diagonal elements of a. + nl = 0 + sq = 0. + do 90 i=1,n + if(a(i,1).gt.tol) go to 90 +c if a sufficient small diagonal element is found, we put it to +c zero. the remainder of the row corresponding to that zero diagonal +c element is then rotated into triangle by givens rotations . +c the rank deficiency is increased by one. + nl = nl+1 + if(i.eq.n) go to 90 + yi = f(i) + do 10 j=1,m1 + h(j) = a(i,j+1) + 10 continue + h(m) = 0. + i1 = i+1 + do 60 ii=i1,n + i2 = min0(n-ii,m1) + piv = h(1) + if(piv.eq.0.) go to 30 + call fpgivs(piv,a(ii,1),cos,sin) + call fprota(cos,sin,yi,f(ii)) + if(i2.eq.0) go to 70 + do 20 j=1,i2 + j1 = j+1 + call fprota(cos,sin,h(j1),a(ii,j1)) + h(j) = h(j1) + 20 continue + go to 50 + 30 if(i2.eq.0) go to 70 + do 40 j=1,i2 + h(j) = h(j+1) + 40 continue + 50 h(i2+1) = 0. + 60 continue +c add to the sum of squared residuals the contribution of deleting +c the row with small diagonal element. + 70 sq = sq+yi**2 + 90 continue +c rank denotes the rank of a. + rank = n-nl +c let b denote the (rank*n) upper trapezoidal matrix which can be +c obtained from the (n*n) upper triangular matrix a by deleting +c the rows and interchanging the columns corresponding to a zero +c diagonal element. if this matrix is factorized using givens +c transformations as b = (r) (u) where +c r is a (rank*rank) upper triangular matrix, +c u is a (rank*n) orthonormal matrix +c then the minimal least-squares solution c is given by c = b' v, +c where v is the solution of the system (r) (r)' v = g and +c g denotes the vector obtained from the old right hand side f, by +c removing the elements corresponding to a zero diagonal element of a. +c initialization. + do 100 i=1,rank + do 100 j=1,m + aa(i,j) = 0. + 100 continue +c form in aa the upper triangular matrix obtained from a by +c removing rows and columns with zero diagonal elements. form in ff +c the new right hand side by removing the elements of the old right +c hand side corresponding to a deleted row. + ii = 0 + do 120 i=1,n + if(a(i,1).le.tol) go to 120 + ii = ii+1 + ff(ii) = f(i) + aa(ii,1) = a(i,1) + jj = ii + kk = 1 + j = i + j1 = min0(j-1,m1) + if(j1.eq.0) go to 120 + do 110 k=1,j1 + j = j-1 + if(a(j,1).le.tol) go to 110 + kk = kk+1 + jj = jj-1 + aa(jj,kk) = a(j,k+1) + 110 continue + 120 continue +c form successively in h the columns of a with a zero diagonal element. + ii = 0 + do 200 i=1,n + ii = ii+1 + if(a(i,1).gt.tol) go to 200 + ii = ii-1 + if(ii.eq.0) go to 200 + jj = 1 + j = i + j1 = min0(j-1,m1) + do 130 k=1,j1 + j = j-1 + if(a(j,1).le.tol) go to 130 + h(jj) = a(j,k+1) + jj = jj+1 + 130 continue + do 140 kk=jj,m + h(kk) = 0. + 140 continue +c rotate this column into aa by givens transformations. + jj = ii + do 190 i1=1,ii + j1 = min0(jj-1,m1) + piv = h(1) + if(piv.ne.0.) go to 160 + if(j1.eq.0) go to 200 + do 150 j2=1,j1 + j3 = j2+1 + h(j2) = h(j3) + 150 continue + go to 180 + 160 call fpgivs(piv,aa(jj,1),cos,sin) + if(j1.eq.0) go to 200 + kk = jj + do 170 j2=1,j1 + j3 = j2+1 + kk = kk-1 + call fprota(cos,sin,h(j3),aa(kk,j3)) + h(j2) = h(j3) + 170 continue + 180 jj = jj-1 + h(j3) = 0. + 190 continue + 200 continue +c solve the system (aa) (f1) = ff + ff(rank) = ff(rank)/aa(rank,1) + i = rank-1 + if(i.eq.0) go to 230 + do 220 j=2,rank + store = ff(i) + i1 = min0(j-1,m1) + k = i + do 210 ii=1,i1 + k = k+1 + stor1 = ff(k) + stor2 = aa(i,ii+1) + store = store-stor1*stor2 + 210 continue + stor1 = aa(i,1) + ff(i) = store/stor1 + i = i-1 + 220 continue +c solve the system (aa)' (f2) = f1 + 230 ff(1) = ff(1)/aa(1,1) + if(rank.eq.1) go to 260 + do 250 j=2,rank + store = ff(j) + i1 = min0(j-1,m1) + k = j + do 240 ii=1,i1 + k = k-1 + stor1 = ff(k) + stor2 = aa(k,ii+1) + store = store-stor1*stor2 + 240 continue + stor1 = aa(j,1) + ff(j) = store/stor1 + 250 continue +c premultiply f2 by the transpoze of a. + 260 k = 0 + do 280 i=1,n + store = 0. + if(a(i,1).gt.tol) k = k+1 + j1 = min0(i,m) + kk = k + ij = i+1 + do 270 j=1,j1 + ij = ij-1 + if(a(ij,1).le.tol) go to 270 + stor1 = a(ij,j) + stor2 = ff(kk) + store = store+stor1*stor2 + kk = kk-1 + 270 continue + c(i) = store + 280 continue +c add to the sum of squared residuals the contribution of putting +c to zero the small diagonal elements of matrix (a). + stor3 = 0. + do 310 i=1,n + if(a(i,1).gt.tol) go to 310 + store = f(i) + i1 = min0(n-i,m1) + if(i1.eq.0) go to 300 + do 290 j=1,i1 + ij = i+j + stor1 = c(ij) + stor2 = a(i,j+1) + store = store-stor1*stor2 + 290 continue + 300 fac = a(i,1)*c(i) + stor1 = a(i,1) + stor2 = c(i) + stor1 = stor1*stor2 + stor3 = stor3+stor1*(stor1-store-store) + 310 continue + fac = stor3 + sq = sq+fac + return + end + + real*8 function fprati(p1,f1,p2,f2,p3,f3) +c given three points (p1,f1),(p2,f2) and (p3,f3), function fprati +c gives the value of p such that the rational interpolating function +c of the form r(p) = (u*p+v)/(p+w) equals zero at p. +c .. +c ..scalar arguments.. + real*8 p1,f1,p2,f2,p3,f3 +c ..local scalars.. + real*8 h1,h2,h3,p +c .. + if(p3.gt.0.) go to 10 +c value of p in case p3 = infinity. + p = (p1*(f1-f3)*f2-p2*(f2-f3)*f1)/((f1-f2)*f3) + go to 20 +c value of p in case p3 ^= infinity. + 10 h1 = f1*(f2-f3) + h2 = f2*(f3-f1) + h3 = f3*(f1-f2) + p = -(p1*p2*h3+p2*p3*h1+p3*p1*h2)/(p1*h1+p2*h2+p3*h3) +c adjust the value of p1,f1,p3 and f3 such that f1 > 0 and f3 < 0. + 20 if(f2.lt.0.) go to 30 + p1 = p2 + f1 = f2 + go to 40 + 30 p3 = p2 + f3 = f2 + 40 fprati = p + return + end + subroutine fprota(cos,sin,a,b) +c subroutine fprota applies a givens rotation to a and b. +c .. +c ..scalar arguments.. + real*8 cos,sin,a,b +c ..local scalars.. + real*8 stor1,stor2 +c .. + stor1 = a + stor2 = b + b = cos*stor2+sin*stor1 + a = cos*stor1-sin*stor2 + return + end + subroutine fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg) +c subroutine fporde sorts the data points (x(i),y(i)),i=1,2,...,m +c according to the panel tx(l)<=x /* To talk to python */ +#include "Numeric/arrayobject.h" /* Access to Numeric */ + +inline int conv_double_to_int_fast(double); +inline int conv_double_to_int_safe(double); + + +static PyObject *closest( PyObject *self, PyObject *args, PyObject *keywds){ + PyArrayObject *ar=NULL, *vals=NULL; + double *x, *v; + int ibest,i,j,nj; + double best; + if(!PyArg_ParseTuple(args,"O!O!", + &PyArray_Type, &ar, /* array args */ + &PyArray_Type, &vals)) /* array args */ + return NULL; + + if(ar->nd != 1 || ar->descr->type_num!=PyArray_DOUBLE){ + PyErr_SetString(PyExc_ValueError, + "First array must be 1d and double"); + return NULL; + } + if(vals->nd != 1 || vals->descr->type_num!=PyArray_DOUBLE){ + PyErr_SetString(PyExc_ValueError, + "Second array must be 1d and double"); + return NULL; + } + + if(ar->strides[0] != sizeof(double) || vals->strides[0] != sizeof(double)){ + PyErr_SetString(PyExc_ValueError, + "Arrays must be flat (strides == sizeof(double))"); + return NULL; + } + + x = (double*) ar->data; + v = (double*) vals->data; + + best=99.; + ibest=0; + nj=vals->dimensions[0]; + for(i=0; idimensions[0];i++){ + for(j=0; jnd != 2 || ubi->descr->type_num!=PyArray_DOUBLE){ + printf("first arg nd %d\n",ubi->nd); + printf("first arg type %d\n",ubi->descr->type_num); + PyErr_SetString(PyExc_ValueError, + "First array must be 3x3 2d and double"); + return NULL; + } + if(gv->nd != 2 || gv->descr->type_num!=PyArray_DOUBLE){ + PyErr_SetString(PyExc_ValueError, + "Second array must be 3xn 2d and double"); + return NULL; + } + +/* if(gv->strides[0] != sizeof(double) || ubi->strides[0] != sizeof(double)){ + PyErr_SetString(PyExc_ValueError, + "Arrays must be flat (strides == sizeof(double))"); + return NULL; + } */ + + if(ubi->dimensions[0] != 3 || ubi->dimensions[1] != 3){ + PyErr_SetString(PyExc_ValueError, + "First array must be 3x3 2d, dimensions problem"); + return NULL; + } + + if(gv->dimensions[1] != 3){ + PyErr_SetString(PyExc_ValueError, + "Second array must be 3xn 2d and double"); + return NULL; + } +/* if(gv->strides[1] != sizeof(double)){ + PyErr_SetString(PyExc_ValueError, + "Second array must be 3xn 2d and double, and you want it aligned!!!"); + return NULL; + } +*/ + + /* for(i=0;i<3;i++) + for(j=0;j<3;j++) + u[i][j] = * (double *) (ubi->data + i*ubi->strides[0] + j*ubi->strides[1]); */ + + u00=* (double *) (ubi->data + 0*ubi->strides[0] + 0*ubi->strides[1]); + u01=* (double *) (ubi->data + 0*ubi->strides[0] + 1*ubi->strides[1]); + u02=* (double *) (ubi->data + 0*ubi->strides[0] + 2*ubi->strides[1]); + u10=* (double *) (ubi->data + 1*ubi->strides[0] + 0*ubi->strides[1]); + u11=* (double *) (ubi->data + 1*ubi->strides[0] + 1*ubi->strides[1]); + u12=* (double *) (ubi->data + 1*ubi->strides[0] + 2*ubi->strides[1]); + u20=* (double *) (ubi->data + 2*ubi->strides[0] + 0*ubi->strides[1]); + u21=* (double *) (ubi->data + 2*ubi->strides[0] + 1*ubi->strides[1]); + u22=* (double *) (ubi->data + 2*ubi->strides[0] + 2*ubi->strides[1]); + + n=0; + tol=tol*tol; + + for(k=0;kdimensions[0];k++){ /* Loop over observed peaks */ + /* Compute hkls of this peak as h = UBI g */ + g0 = * (double *) (gv->data + k*gv->strides[0] + 0*gv->strides[1]); + g1 = * (double *) (gv->data + k*gv->strides[0] + 1*gv->strides[1]); + g2 = * (double *) (gv->data + k*gv->strides[0] + 2*gv->strides[1]); + h0 = u00*g0 + u01*g1 + u02*g2; + h1 = u10*g0 + u11*g1 + u12*g2; + h2 = u20*g0 + u21*g1 + u22*g2; + t0=h0-conv_double_to_int_fast(h0); + t1=h1-conv_double_to_int_fast(h1); + t2=h2-conv_double_to_int_fast(h2); +/* assert(t[0]<0.501); + assert(t[1]<0.501); + assert(t[2]<0.501); */ + sumsq = t0*t0+t1*t1+t2*t2; + if (sumsq < tol) + n=n+1; + /* + printf("k=%d g=(",k); + printf(" %6.3f ",g[0]); + printf(" %6.3f ",g[1]); + printf(" %6.3f ) h=(",g[2]); + printf(" %6.3f ",h[0]); + printf(" %6.3f ",h[1]); + printf(" %6.3f ) ",h[2]); + printf("int(h)=(%3d %3d %3d) sumsq=%g n=%d\n",conv_double_to_int_fast(h[0]), + conv_double_to_int_fast(h[1]),conv_double_to_int_fast(h[2]),sumsq,n); + */ + } + + + return Py_BuildValue("i",n); +} + +inline int conv_double_to_int_safe(double x){ + int a; + a = floor(x+0.5); + return a; +} + + + +inline int conv_double_to_int_fast(double x){ + /* This was benched as about ten times faster than the safe mode!! */ + const int p=52; + const double c_p1 = (1L << (p/2)); + const double c_p2 = (1L << (p-p/2)); + const double c_mul = c_p1 * c_p2; + const double cs = 1.5*c_mul; + x+=cs; + const int a = *(int *)(&x); + return (a); +} + +/* Make a UBI matrix from two peaks and score it ?? */ + + + +static PyMethodDef closestMethods[] = { + { "closest", (PyCFunction) closest, METH_VARARGS, + "int,float closest(x,v)\nFind element in x closest to one in v\nx,v 1D flat Float arrays"}, + { "score", (PyCFunction) score, METH_VARARGS, + "int score(ubi,gv,tol)\nFind number of gv which are integer h within tol"}, + { NULL, NULL, 0, NULL } +}; + +void initclosest(void) +{ + PyObject *m, *d; + m=Py_InitModule("closest",closestMethods); + import_array(); + d=PyModule_GetDict(m); +} + + diff --git a/src/connectedpixels.c b/src/connectedpixels.c new file mode 100644 index 00000000..cc3c374a --- /dev/null +++ b/src/connectedpixels.c @@ -0,0 +1,455 @@ + + + + +/* +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +*/ + + +#include /* To talk to python */ +#include "Numeric/arrayobject.h" /* Access to Numeric */ +#include "dset.h" +#include +#include +#include +/* +int * dset_initialise(int size); +int * dset_new(int ** S); +void dset_makeunion(int * S, int r1, int r2); +void dset_link(int *S, int r1, int r2); +int dset_find(int x, int * S); +*/ + +#undef PARANOID + +int getval(char *p,int type); + +int getval(char *p,int type){ + switch (type){ + case PyArray_CHAR : return *(char *)p*1; + case PyArray_SBYTE : return *(signed char *)p*1; + case PyArray_SHORT : return *(short *)p*1; + case PyArray_INT : return *(int *)p*1; + case PyArray_LONG : return *(long *)p*1; + case PyArray_FLOAT : return *(float *)p*1; + case PyArray_DOUBLE : return *(double *)p*1; +#ifdef PyArray_UNSIGNED_TYPES + case PyArray_UBYTE : return *(unsigned char *)p*1; + case PyArray_USHORT : return *(unsigned short*)p*1; + case PyArray_UINT : return *(unsigned int *)p*1; +#endif + } + printf("Oh bugger in getval - unrecongnised numeric type\n"); + exit(1); + return 0; +} + +void ptype(int type){ + printf("Your input type was "); + switch (type){ + case PyArray_CHAR : printf("PyArray_CHAR *(char *)\n");break; + case PyArray_SBYTE : printf("PyArray_SBYTE *(signed char *)\n");break; + case PyArray_SHORT : printf("PyArray_SHORT *(short *)\n");break; + case PyArray_INT : printf("PyArray_INT *(int *)\n");break; + case PyArray_LONG : printf("PyArray_LONG *(long *)\n");break; + case PyArray_FLOAT : printf("PyArray_FLOAT *(float *)\n");break; + case PyArray_DOUBLE : printf("PyArray_DOUBLE *(double *)\n");break; +#ifdef PyArray_UNSIGNED_TYPES + case PyArray_UBYTE : printf("PyArray_UBYTE *(unsigned char *)\n");break; + case PyArray_USHORT : printf("PyArray_USHORT *(unsigned short*)\n");break; + case PyArray_UINT : printf("PyArray_UINT *(unsigned int *)\n");break; +#endif + } +} + + +/* make an image of peak assignement for pixels */ +static PyObject * connectedpixels (PyObject *self, PyObject *args, PyObject *keywds) +{ + PyArrayObject *dataarray=NULL,*results=NULL; /* in (not modified) and out */ + + int i,j,k,l,ival,f,s,np; + int *T; /* for the disjoint set copy */ + int verbose=0,type; /* whether to print stuff and the type of the input array */ + int percent,npover; /* for the progress indicator in verbose mode */ + static char *kwlist[] = {"data","results","threshold","verbose", NULL}; + clock_t tv1,tv1point5,tv2,tv3,tv4; + double time; + float threshold; + if(!PyArg_ParseTupleAndKeywords(args,keywds, "O!O!f|i",kwlist, + &PyArray_Type, &dataarray, /* array args */ + &PyArray_Type, &results, /* array args */ + &threshold, &verbose)) /* threshold and optional verbosity */ + return NULL; + if(verbose!=0){ + printf("\n\nHello there from connected pixels, you wanted the verbose output...\n"); + } + + tv1=clock(); + /* Check array is two dimensional and float */ + if(dataarray->nd != 2){ + PyErr_SetString(PyExc_ValueError, + "Array must be 2d, first arg problem"); + return NULL; + } + type=dataarray->descr->type_num; + if(verbose!=0)ptype(type); + if(verbose!=0)printf("Thresholding at level %f\n",threshold); + /* make an array to hold the integer peak assignments */ + /* PyArray_FromDims(int n_dimensions, int dimensions[n_dimensions], int type_num) */ + /* Make a new array to hold the results... how to do optional arguments?? */ + /* results=(PyArrayObject*)PyArray_FromDims(dataarray->nd, dataarray->dimensions, PyArray_INT); */ + /* if(results==NULL){*/ + /* printf("Could not make a results array\nGiving up\n");*/ + /* exit(2); } */ + /* if(verbose!=0)printf("Made results array\n");*/ + + if(results->nd != 2){ + PyErr_SetString(PyExc_ValueError, + "Array must be 2d, second arg problem"); + return NULL; + } + if( results->dimensions[0] != dataarray->dimensions[0] || + results->dimensions[1] != dataarray->dimensions[1] ){ + PyErr_SetString(PyExc_ValueError, + "Arrays must have same shape"); + return NULL; + } + + S = dset_initialise(16384); /* Default number before reallocation */ + if(verbose!=0)printf("Initialised the disjoint set\n"); + npover=0; /* number of pixels over threshold */ + /* Decide on fast/slow loop - inner versus outer */ + if(dataarray->strides[0] > dataarray->strides[1]) { + f=1; s=0;} + else { + f=0; s=1; + } + if (verbose!=0){ + printf("Fast index is %d, slow index is %d, ",f,s); + printf("strides[0]=%d, strides[1]=%d\n",dataarray->strides[0],dataarray->strides[1]); + } + percent=(results->dimensions[s]-1)/80.0; + if(percent < 1)percent=1; + tv1point5=clock(); + if(verbose!=0){ + time=(tv1point5-tv1)/CLOCKS_PER_SEC; + printf("Ready to scan image, setup took %g and clock is only good to ~0.015 seconds\n",time); + } + + if(verbose!=0)printf("Scanning image\n"); + for( i = 0 ; i < (results->dimensions[s]) ; i++ ){ /* i,j is looping along the indices data array */ + if(verbose!=0 && (i%percent == 0) )printf("."); + + + for( j = 0 ; j < (results->dimensions[f]) ; j++ ){ + + /* Set result for this pixel to zero */ + + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = 0; + + ival= (* (unsigned short *) (dataarray->data + i*dataarray->strides[s] + j*dataarray->strides[f])) ; + /* val=getval((dataarray->data + i*dataarray->strides[s] + j*dataarray->strides[f]),type); */ + if( ival > threshold) { + npover++; + k=0;l=0; + /* peak needs to be assigned */ + /* i-1, j-1, assign same index */ + if(i!=0 && j!=0) + if( ( k=*(int *)(results->data + (i-1)*results->strides[s] + (j-1)*results->strides[f] )) >0 ) { + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = k; + l=k; /* l holds the assignment of the current pixel */ + } + /* i-1,j */ + if(i!=0) + if( ( k=*(int *)(results->data + (i-1)*results->strides[s] + j*results->strides[f] )) >0 ) { + if(l==0){ + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = k; + l=k; /* l holds the assignment of the current pixel */ + } else { + if(l!=k)dset_makeunion(S,k,l); + } + } + /* i-1, j+1 */ + if(i!=0 && j!=(results->dimensions[f]-1)) + if( ( k=*(int *)(results->data + (i-1)*results->strides[s] + (j+1)*results->strides[f] )) >0 ) { + if(l==0){ + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = k; + l=k; /* l holds the assignment of the current pixel */ + } else { + if(l!=k)dset_makeunion(S,k,l); + } + } + /* i, j-1 */ + if(j!=0) + if( ( k=*(int *)(results->data + i*results->strides[s] + (j-1)*results->strides[f] )) >0 ) { + if(l==0){ + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = k; + l=k; /* l holds the assignment of the current pixel */ + } else { + if(l!=k)dset_makeunion(S,k,l); + } + } + if(l==0){ /* pixel has no neighbours thus far */ + S = dset_new(&S); + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = S[S[0]-1]; + } + } /* active pixel */ + else { /* inactive pixel - zero in results */ + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f])) = 0; + } + } /* j */ + } /* i */ + tv2=clock(); + if(verbose!=0){ + printf("\nFinished scanning image, now make unique list and sums\n"); + printf("Number of pixels over threshold was : %d\n",npover); + time=(tv2-tv1)/CLOCKS_PER_SEC; + printf("That took %f seconds and %f in total being aware that the clock is only good to ~0.015 seconds\n",(tv2-tv1point5)/CLOCKS_PER_SEC,time); + } + /* First, count independent peaks */ + + /* count peaks */ + np=S[S[0]-1]; + /* Now make each T[i] contain the unique ascending integer for the set */ + T=NULL; + if(verbose!=0)printf("Now I want to malloc T=%p at size np=%d and sizeof(int)=%d\n",T,np,sizeof(int)); + if( ( T=(int *) ( malloc( (np+3)*sizeof(int) ) ) ) == NULL){ + printf("Memory allocation error in connected pixels\n"); + return 0; + } + if(verbose!=0)printf("Now T should be allocated, it is at %p\n",T); + np=0; + for(i=1;i=i && verbose){ + printf("Oh dear - there was a problem compressing the disjoint set, j=%d, i=%d \n",j,i); + } + if(S[j]!=j && verbose!=0){ + printf("Oh dear - the disjoint set is squiff,S[j]=%d,j=%d\n",S[j],j); + } + } + } + if(verbose!=0)printf("\n"); + tv3=clock(); + if(verbose!=0){ + printf("Compressing list, time=%g, total time=%g\n",(tv3-tv2)/CLOCKS_PER_SEC,(tv3-tv1)/CLOCKS_PER_SEC); + } + + if(verbose)printf("Found %d peaks\nNow assigning unique identifiers\n",np); + + + /* Now loop through the results image again, assigning the corrected numbers to each pixel */ + for( i = 0 ; i < (results->dimensions[s]) ; i++ ){ /* i,j is looping along the indices data array */ + for( j = 0 ; j < (results->dimensions[f]) ; j++ ){ + ival=*(int *)(results->data + i*results->strides[s] + j*results->strides[f]); + if( ival > 0){ + if(T[ival]==0){ + printf("Something buggered up, ival=%d, T[ival]=%d, dset_find(ival,S)=%d\n",ival,T[ival],dset_find(ival,S)); + /* free(S); */ + exit(2); + return NULL; + } + (*(int *)(results->data + i*results->strides[s] + j*results->strides[f]))=T[ival]; + } + } + } + + + + + if(verbose!=0)printf("Freeing S=%p and T=%p\n",S,T); + free(S); + free(T); + if(verbose!=0){ + printf("Freed S=%p and T=%p\n",S,T); + printf("Finished assigning unique values and tidied up, now returning"); + tv4=clock(); + time=(tv4-tv3)/CLOCKS_PER_SEC; + printf("\nThat took %g seconds, being aware that the clock is only good to ~0.015 seconds\n",time); + printf("The total time in this routine was about %f\n",(tv4-tv1)/CLOCKS_PER_SEC); + } + +/* return PyArray_Return(results); */ +/* Py_DECREF(results);*/ + /* Py_DECREF(dataarray);*/ + return Py_BuildValue("i", np); +} + + +/* BLOB PROPERTIES BIT - things to measure are: */ +/* Centre of mass (x,y) - optionally use x,y image */ +/* Amount of intensity in blob */ +/* std dev of intensity in blob */ +/* ...? */ + + +static PyObject * blobproperties (PyObject *self, PyObject *args, PyObject *keywds) +{ + PyArrayObject *dataarray=NULL,*blobarray=NULL; /* in (not modified) */ + static char *kwlist[] = {"data","blob","npeaks","verbose", NULL}; + PyArrayObject *sum,*sumsq,*com0,*com1,*npix, *com00, *com11, *com01; + double *asum,*asumsq,*acom0,*acom1, *acom00, *acom11, *acom01 ; + double fval; + int np,verbose=0,type,f,s,peak,bad; + int i,j,safelyneed,*anpix, percent; + if(!PyArg_ParseTupleAndKeywords(args,keywds, "O!O!i|i",kwlist, + &PyArray_Type, &dataarray, /* array args - data */ + &PyArray_Type, &blobarray, /* blobs */ + &np, /* Number of peaks to treat */ + &verbose)) /* threshold and optional verbosity */ + return NULL; + + /* Check array is two dimensional */ + if(dataarray->nd != 2){ + PyErr_SetString(PyExc_ValueError, + "data array must be 2d, first arg problem"); + return NULL; + } + if(verbose!=0)printf("Welcome to blobproperties\n"); + /* Check array is two dimensional and float */ + if(blobarray->nd != 2 && blobarray->descr->type_num != PyArray_INT){ + PyErr_SetString(PyExc_ValueError, + "Blob array must be 2d and integer, second arg problem"); + return NULL; + } + type=dataarray->descr->type_num; + /* Decide on fast/slow loop - inner versus outer */ + if(blobarray->strides[0] > blobarray->strides[1]) { + f=1; s=0;} + else { + f=0; s=1; + } + if (verbose!=0){ + printf("Fast index is %d, slow index is %d, ",f,s); + printf("strides[0]=%d, strides[1]=%d\n",dataarray->strides[0],dataarray->strides[1]); + } + /* results arrays */ + safelyneed=np+3; + npix = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_INT); + sumsq = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + sum = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + com0 = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + com1 = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + com00 = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + com01 = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + com11 = (PyArrayObject *)PyArray_FromDims(1,&safelyneed,PyArray_DOUBLE); + + if (npix == NULL || sumsq == NULL || sum == NULL || com0 == NULL || com1 == NULL) goto fail; + anpix = (int *) npix->data; + asumsq = (double *) sumsq->data; + asum = (double *) sum->data; + acom0 = (double *) com0->data; + acom1 = (double *) com1->data; + acom00 = (double *) com00->data; + acom01 = (double *) com01->data; + acom11 = (double *) com11->data; + + for ( i=0 ; idimensions[s]-1)/80.0; + if(percent < 1)percent=1; + bad=0; + if(verbose!=0)printf("Scanning image\n"); + for( i = 0 ; i <= (blobarray->dimensions[s]-1) ; i++ ){ /* i,j is looping along the indices data array */ + if(verbose!=0 && (i%percent == 0) )printf("."); + for( j = 0 ; j <= (blobarray->dimensions[f]-1) ; j++ ){ + peak=* (int *) (blobarray->data + i*blobarray->strides[s] + j*blobarray->strides[f]); + if( peak > 0 && peak <=np ) { + fval=(double)(getval((dataarray->data + i*dataarray->strides[s] + j*dataarray->strides[f]),type)); + asum[peak] =asum[peak]+fval; + asumsq[peak]=asumsq[peak]+fval*fval; + acom0[peak] =acom0[peak]+i*fval; + acom1[peak] =acom1[peak]+j*fval; + acom00[peak] =acom00[peak]+i*i*fval; + acom01[peak] =acom01[peak]+i*j*fval; + acom11[peak] =acom11[peak]+j*j*fval; + anpix[peak] =anpix[peak]++; +#ifdef PARANOID + if(verbose!=0){ + printf("peak=%d i=%d j=%d val=%f acom00[peak]=%f 01=%f 11=%f\n",peak,i,j,fval,acom00[peak], + acom01[peak],acom11[peak]);} +#endif + } + else{ + if(peak!=0){ + bad++; + if(verbose!=0 && bad<10){printf("Found %d in your blob image at i=%d, j=%d\n",peak,i,j);} + } + } + } /* j */ + } /* i */ + if(verbose){ + printf("\nFound %d bad pixels in the blob image\n",bad); + } + return Py_BuildValue("OOOOOOOO", PyArray_Return(npix), + PyArray_Return(sum),PyArray_Return(sumsq), + PyArray_Return(com0),PyArray_Return(com1),PyArray_Return(com00), + PyArray_Return(com01),PyArray_Return(com11) ); + fail: + Py_XDECREF(npix); + Py_XDECREF(sum); + Py_XDECREF(sumsq); + Py_XDECREF(com0); + Py_XDECREF(com1); + Py_XDECREF(com00); + Py_XDECREF(com01); + Py_XDECREF(com11); + return NULL; +} + + + + +static PyMethodDef connectedpixelsMethods[] = { + {"connectedpixels", (PyCFunction) connectedpixels, METH_VARARGS | METH_KEYWORDS, + "Assign connected pixels in image (first arg) above threshold (float second arg) to return image of peak number.Optional arguement verbose!=0 prints some info as it runs"}, + {"blobproperties", (PyCFunction) blobproperties, METH_VARARGS | METH_KEYWORDS, + "Computes various properties of a blob image (created by connecetedpixels)"}, + + {NULL, NULL, 0, NULL} /* setinel */ +}; + +void +initconnectedpixels(void) +{ + PyObject *m, *d; + + m=Py_InitModule("connectedpixels", connectedpixelsMethods); + import_array(); + d=PyModule_GetDict(m); +} diff --git a/src/dset.h b/src/dset.h new file mode 100644 index 00000000..536f3f47 --- /dev/null +++ b/src/dset.h @@ -0,0 +1,98 @@ + + +/* +# ImageD11_v0.4 Software for beamline ID11 +# Copyright (C) 2005 Jon Wright +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + + +#ifndef _dset_h +#define _dset_h +int * dset_initialise(int size); /* array to hold real values of each */ +int * dset_new(int ** S); +void dset_makeunion(int * S, int r1, int r2); +void dset_link(int *S, int r1, int r2); +int dset_find(int x, int * S); + +int *S; + +int * dset_initialise(int size){ + int i; + /* printf("Initialising the disjoint set at size %d, S=%x\n",size,S); */ + S=(int *) ( malloc(size*sizeof(int) ) ); + if (S==NULL){ + printf("Memory allocation error in dset_initialise\n"); + exit(1); + } + for(i=0;ilength){ + S=(int *)(realloc(S,length*2*sizeof(int))); + if(S==NULL){ + printf("Memory allocation error in dset_new\n"); + exit(1); + } +/* printf("Realloced S to %d in dset_new\n",length*2); */ + (int *) pS = S; + S[0]=length*2; + S[length-1]=0; + for(i=length-1;ir2){S[r1]=r2;} + if(r1 +#include "Numeric/arrayobject.h" + +#define BISPEV bispev_ +#define PARDER parder_ +#define SURFIT surfit_ + +void BISPEV(double*,int*,double*,int*,double*,int*,int*,double*,int*, + double*,int*,double*,double*,int*,int*,int*,int*); +void PARDER(double*,int*,double*,int*,double*,int*,int*,int*,int*,double*, + int*,double*,int*,double*,double*,int*,int*,int*,int*); +void SURFIT(int*,int*,double*,double*,double*,double*,double*,double*,double*,double*, + int*,int*,double*,int*,int*,int*,double*,int*,double*,int*,double*,double*, + double*,double*,int*,double*,int*,int*,int*,int*); + +static char doc_bispev[] = " [z,ier] = _bispev(tx,ty,c,kx,ky,x,y,nux,nuy)"; + +static PyObject *_bispev(PyObject *dummy, PyObject *args) { + int nx,ny,kx,ky,mx,my,lwrk,*iwrk,kwrk,ier,lwa,mxy,nux,nuy; + double *tx,*ty,*c,*x,*y,*z,*wrk,*wa = NULL; + PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_z = NULL,*ap_tx = NULL,\ + *ap_ty = NULL,*ap_c = NULL; + PyObject *x_py = NULL,*y_py = NULL,*c_py = NULL,*tx_py = NULL,*ty_py = NULL; + if (!PyArg_ParseTuple(args, "OOOiiOOii",&tx_py,&ty_py,&c_py,&kx,&ky, + &x_py,&y_py,&nux,&nuy)) + return NULL; + ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1); + ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, PyArray_DOUBLE, 0, 1); + ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1); + ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, PyArray_DOUBLE, 0, 1); + ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, PyArray_DOUBLE, 0, 1); + if (ap_x == NULL || ap_y == NULL || ap_c == NULL || ap_tx == NULL \ + || ap_ty == NULL) goto fail; + x = (double *) ap_x->data; + y = (double *) ap_y->data; + c = (double *) ap_c->data; + tx = (double *) ap_tx->data; + ty = (double *) ap_ty->data; + nx = ap_tx->dimensions[0]; + ny = ap_ty->dimensions[0]; + mx = ap_x->dimensions[0]; + my = ap_y->dimensions[0]; + mxy = mx*my; + ap_z = (PyArrayObject *)PyArray_FromDims(1,&mxy,PyArray_DOUBLE); + z = (double *) ap_z->data; + if (nux || nuy) + lwrk = mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1); + else + lwrk = mx*(kx+1)+my*(ky+1); + kwrk = mx+my; + lwa = lwrk+kwrk; + if ((wa = (double *)malloc(lwa*sizeof(double)))==NULL) { + PyErr_NoMemory(); + goto fail; + } + wrk = wa; + iwrk = (int *)(wrk+lwrk); + if (nux || nuy) + PARDER(tx,&nx,ty,&ny,c,&kx,&ky,&nux,&nuy,x,&mx,y,&my,z,wrk,&lwrk,iwrk,&kwrk,&ier); + else + BISPEV(tx,&nx,ty,&ny,c,&kx,&ky,x,&mx,y,&my,z,wrk,&lwrk,iwrk,&kwrk,&ier); + if (wa) free(wa); + Py_DECREF(ap_x); + Py_DECREF(ap_y); + Py_DECREF(ap_c); + Py_DECREF(ap_tx); + Py_DECREF(ap_ty); + return Py_BuildValue("Ni",PyArray_Return(ap_z),ier); + fail: + if (wa) free(wa); + Py_XDECREF(ap_x); + Py_XDECREF(ap_y); + Py_XDECREF(ap_z); + Py_XDECREF(ap_c); + Py_XDECREF(ap_tx); + Py_XDECREF(ap_ty); + return NULL; +} + +static char doc_surfit[] = + " [tx,ty,c,o] = " + "_surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,iopt,s,eps,tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)"; + +static PyObject *_surfit(PyObject *dummy, PyObject *args) { + int iopt,m,kx,ky,nxest,nyest,nx,ny,lwrk1,lwrk2,*iwrk,kwrk,ier,lwa,nxo,nyo,\ + i,lc,lcest,nmax; + double *x,*y,*z,*w,xb,xe,yb,ye,s,*tx,*ty,*c,fp,*wrk1,*wrk2,*wa = NULL,eps; + PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_z,*ap_w = NULL,\ + *ap_tx = NULL,*ap_ty = NULL,*ap_c = NULL; + PyArrayObject *ap_wrk = NULL; + PyObject *x_py = NULL,*y_py = NULL,*z_py = NULL,*w_py = NULL,\ + *tx_py = NULL,*ty_py = NULL; + PyObject *wrk_py=NULL; + nx=ny=ier=nxo=nyo=0; + if (!PyArg_ParseTuple(args, "OOOOddddiiiddOOiiOii",\ + &x_py,&y_py,&z_py,&w_py,&xb,&xe,\ + &yb,&ye,&kx,&ky,&iopt,&s,&eps,&tx_py,&ty_py,&nxest,&nyest,\ + &wrk_py,&lwrk1,&lwrk2)) return NULL; + ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1); + ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, PyArray_DOUBLE, 0, 1); + ap_z = (PyArrayObject *)PyArray_ContiguousFromObject(z_py, PyArray_DOUBLE, 0, 1); + ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, PyArray_DOUBLE, 0, 1); + ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, PyArray_DOUBLE, 0, 1); + /*ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, PyArray_INT, 0, 1);*/ + if (ap_x == NULL || ap_y == NULL || ap_z == NULL || ap_w == NULL \ + || ap_wrk == NULL) goto fail; + x = (double *) ap_x->data; + y = (double *) ap_y->data; + z = (double *) ap_z->data; + w = (double *) ap_w->data; + m = ap_x->dimensions[0]; + nmax=nxest; + if (nmaxdimensions[0]; + ny = nyo = ap_ty->dimensions[0]; + memcpy(tx,ap_tx->data,nx*sizeof(double)); + memcpy(ty,ap_ty->data,ny*sizeof(double)); + } + if (iopt==1) { + lc = (nx-kx-1)*(ny-ky-1); + memcpy(wrk1,ap_wrk->data,lc*sizeof(double)); + /*memcpy(iwrk,ap_iwrk->data,n*sizeof(int));*/ + } + SURFIT(&iopt,&m,x,y,z,w,&xb,&xe,&yb,&ye,&kx,&ky,&s,&nxest,&nyest,&nmax, + &eps,&nx,tx,&ny,ty,c,&fp,wrk1,&lwrk1,wrk2,&lwrk2,iwrk,&kwrk,&ier); + i=0; + while ((ier>10) && (i++<5)) { + lwrk2=ier; + if ((wrk2 = (double *)malloc(lwrk2*sizeof(double)))==NULL) { + PyErr_NoMemory(); + goto fail; + } + SURFIT(&iopt,&m,x,y,z,w,&xb,&xe,&yb,&ye,&kx,&ky,&s,&nxest,&nyest,&nmax, + &eps,&nx,tx,&ny,ty,c,&fp,wrk1,&lwrk1,wrk2,&lwrk2,iwrk,&kwrk,&ier); + if (wrk2) free(wrk2); + } + if (ier==10) goto fail; + lc = (nx-kx-1)*(ny-ky-1); + ap_tx = (PyArrayObject *)PyArray_FromDims(1,&nx,PyArray_DOUBLE); + ap_ty = (PyArrayObject *)PyArray_FromDims(1,&ny,PyArray_DOUBLE); + ap_c = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE); + if (ap_tx == NULL || ap_ty == NULL || ap_c == NULL) goto fail; + if ((iopt==0)||(nx>nxo)||(ny>nyo)) { + ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE); + if (ap_wrk == NULL) goto fail; + /*ap_iwrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_INT);*/ + } + memcpy(ap_tx->data,tx,nx*sizeof(double)); + memcpy(ap_ty->data,ty,ny*sizeof(double)); + memcpy(ap_c->data,c,lc*sizeof(double)); + memcpy(ap_wrk->data,wrk1,lc*sizeof(double)); + /*memcpy(ap_iwrk->data,iwrk,n*sizeof(int));*/ + if (wa) free(wa); + Py_DECREF(ap_x); + Py_DECREF(ap_y); + Py_DECREF(ap_z); + Py_DECREF(ap_w); + return Py_BuildValue("NNN{s:N,s:i,s:d}",PyArray_Return(ap_tx),\ + PyArray_Return(ap_ty),PyArray_Return(ap_c),\ + "wrk",PyArray_Return(ap_wrk),\ + "ier",ier,"fp",fp); + fail: + if (wa) free(wa); + Py_XDECREF(ap_x); + Py_XDECREF(ap_y); + Py_XDECREF(ap_z); + Py_XDECREF(ap_w); + Py_XDECREF(ap_tx); + Py_XDECREF(ap_ty); + Py_XDECREF(ap_wrk); + /*Py_XDECREF(ap_iwrk);*/ + return NULL; +} + + + +static struct PyMethodDef splines_methods[] = { + {"_bispev", _bispev, METH_VARARGS, doc_bispev}, + {"_surfit", _surfit, METH_VARARGS, doc_surfit}, + + {NULL, NULL, 0, NULL} + }; + +void init_splines(void) { + PyObject *m, *d, *s; + m = Py_InitModule("_splines", splines_methods); + import_array(); + d = PyModule_GetDict(m); + s = PyString_FromString(" 0.1 "); + PyDict_SetItemString(d, "__version__", s); + Py_DECREF(s); + if (PyErr_Occurred()) + Py_FatalError("can't initialize module splines"); +}