Initial revision

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<=dx<kx): raise ValueError,"0<=dx=%d<kx=%d must hold"%(dx,kx)
+   if not (0<=dy<ky): raise ValueError,"0<=dy=%d<ky=%d must hold"%(dy,ky)
+   x,y=map(myasarray,[x,y])
+   if (len(x.shape) != 1) or (len(y.shape) != 1):
+       raise ValueError, "First two entries should be rank-1 arrays."
+   z,ier=_splines._bispev(tx,ty,c,kx,ky,x,y,dx,dy)
+   if ier==10: raise ValueError,"Invalid input data"
+   if ier: raise TypeError,"An error occurred"
+   z.shape=len(x),len(y)
+   if len(z)>1: 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
+

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)
+
+