First commit of the reciprocal space reconstruction code

ImageD11/ImageD11_file_series.py

@@ -0,0 +1,142 @@
+
+"""
+To be moved to fabio sometime
+"""
+
+import fabio.file_series
+import fabio.fabioimage
+import fabio.openimage
+import numpy
+
+def get_options(parser):
+
+    parser.add_option("-5","--hdf5",action="store", type="string",
+                      dest = "hdf5", default = None,
+                      help = "hdf file containing input image series")
+    # or, eventually:
+    # stem, first, last, format, (omegas better be in the headers)
+    parser.add_option("-n","--stem",action="store", type="string",
+                      dest = "stem", default = None,
+                      help = "stem name for input image series")
+    parser.add_option("-f","--first",action="store", type="int",
+                      dest = "first", default = None,
+                      help = "first number for input image series")
+    parser.add_option("-l","--last",action="store", type="int",
+                      dest = "last", default = None,
+                      help = "last number for input image series")
+    parser.add_option("--ndigits", action="store", type="int",
+                      dest = "ndigits", default = 4,
+                      help = "Number of digits in file numbering [4]")
+    parser.add_option("-P", "--padding", action="store",
+                      type="choice", choices=["Y","N"], 
+                      default="Y", dest="padding",
+                      help="Is the image number to padded Y|N, e.g. "\
+                          "should 1 be 0001 or just 1 in image name, default=Y")
+    parser.add_option("-F","--format",action="store", type="string",
+                      dest = "format", default = ".edf",
+                      help = "format [.edf] for input image series")
+
+    parser.add_option("-O", "--flood", action="store", type="string",
+                      dest = "flood", default = None,
+                      help = "Flood")
+
+    parser.add_option("-d", "--dark", action="store", type="string",
+                      dest = "dark", default = None,
+                      help = "Dark image")
+    return parser
+
+
+def get_series_from_hdf( hdf_file, dark = None, flood = None ):
+    groups = hdf_file.listnames()
+    for group in groups:
+        imagenames = hdf_file[group].listnames()
+        for image in imagenames:
+            im = hdf_file[group][image]
+            om = float(im.attrs['Omega'])
+            data = im[:,:]
+            if (dark, flood) is not (None, None):
+                data = data.astype(numpy.float32)
+            if dark != None:
+                numpy.subtract( data, dark, data )
+            if flood != None:
+                numpy.divide( data, flood, data )
+            yield fabio.fabioimage.fabioimage( data = data,
+                                               header = {
+                    'Omega': om } )
+
+def series_from_fabioseries( fabioseries, dark=None, flood=None ):
+    for filename in fabioseries:
+        fim = fabio.openimage.openimage(filename)
+        if (dark, flood) is not (None, None):
+            fim.data = fim.data.astype(numpy.float32)
+        if dark != None:
+            numpy.subtract( fim.data, dark, fim.data )
+        if flood != None:
+            numpy.divide( fim.data, flood, fim.data )
+        fim.header['Omega'] = float(fim.header['Omega'])
+        yield fim
+
+
+
+def get_series_from_stemnum( options, args, dark = None, flood = None ):
+    """
+    Returns a file series thing - not a fabio one
+    """
+    if options.format in ['bruker', 'BRUKER', 'Bruker']:
+        extn = ""
+    elif options.format == 'GE':
+        extn = ""
+    else:
+        extn = options.format
+
+    fso = fabio.file_series.numbered_file_series(
+        options.stem,
+        options.first,
+        options.last,
+        extn,
+        digits = options.ndigits,
+        padding = options.padding )
+    return series_from_fabioseries( fso , dark, flood )
+
+
+def get_series_from_options( options, args ):
+    """
+    Returns a file series thing - not a fabio one
+
+    This gives back a fabioimage object with dark and flood
+    corrected data
+    """
+
+    try:
+        if options.dark is not None:
+            dark = fabio.openimage( options.dark ).data
+        else:
+            dark = None
+    except:
+        print "Problem with your dark",options.dark
+        raise
+
+    try:
+        if options.flood is not None:
+            flood = fabio.openimage( options.flood ).data
+        else:
+            flood = None
+    except:
+        print "Problem with your flood",options.flood
+        raise
+
+
+    if len(args) > 0 :
+        # We assume unlabelled arguments are filenames 
+        fso = fabio.file_series.file_series(args)
+        return series_from_fabioseries( fso, dark, flood )
+
+    if options.hdf5 is not None:
+        hf = h5py.File(options.hdf5)
+        # print "Getting images from",options.hdf5
+        return get_series_from_hdf( hf, dark, flood )
+
+    return get_series_from_stemnum( options, args,
+                                     dark, flood) 
+
+

ImageD11/rsv.py

@@ -0,0 +1,132 @@
+
+"""
+Reciprocal Space Volume
+
+This class is just for holding a volume and loading/saving it to disk
+Another class (rsv_mapper) will add images into it.
+"""
+
+
+# ImageD11_v1.0 Software for beamline ID11
+# Copyright (C) 2005-2007  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 FITESS 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  0211-1307  USA
+
+
+import logging, numpy, h5py
+
+class rsv(object):
+    """
+    A reciprocal space volume
+    """
+    def __init__(self, dimensions , **kwds ):
+        """
+        dimensions = NX*NY*NZ grid for the space
+        uspace = a 3x3 matrix describing the grid
+         eg:
+            pixel at vol[i,j,k] comes from:
+               [i,j,k] = (uspace).gvec
+          gvec is the scattering vector in reciprocal space
+          so uspace are vectors in real space 
+        uorigin = a 3 vector giving the position of the [0,0,0] pixel
+        """
+        assert len(dimensions) == 3
+        self.SIG = None # signal
+        self.MON = None # monitor
+        self.NR =  dimensions # dimensions      
+        self.NORMED = None
+        self.metadata = kwds
+        self.allocate_vol()
+
+    def allocate_vol( self ):
+        """ 
+        Allocates memory for a volume data
+        """
+        if self.NR is None:
+            raise Exception("Cannot allocate rsv")
+        total = self.NR[0]*self.NR[1]*self.NR[2]
+        logging.info("rsv: memory used = %.2f MB"%(total*8.0/1024/1024))
+        self.SIG = numpy.zeros( total, numpy.float32 )
+        self.MON = numpy.zeros( total, numpy.float32 )
+        self.metadata = {}
+
+    def normalise( self ):
+        """
+        Return the normalised but avoid divide by zero
+        """
+        self.NORMED = numpy.where( self.MON > 0.1, 
+                              self.SIG/(self.MON+1e-32), 
+                              0.0)
+
+
+def writevol(vol, filename):
+    """
+    Write volume in vol to filename
+    
+    Compress -1 is for the zeros, which there might be a lot of
+    """
+    if not isinstance( vol, rsv ):
+        raise Exception("First arg to writevol should be an rsv object")
+    if None in [vol.NR, vol.SIG, vol.MON]:
+        raise Exception("Cannot save rsv, has not data in it")
+    volout = h5py.File( filename,"w")
+    if vol.SIG.dtype != numpy.float32:
+        logging.warning("rsv SIG was not float32, converting")
+        vol.SIG = vol.SIG.astype(numpy.float32)
+    volout.create_dataset( "signal",
+                           (vol.NR[0],vol.NR[1],vol.NR[2]),
+                           vol.SIG.dtype,
+                           data = vol.SIG,
+                           compression = 'gzip', 
+                           compression_opts = 1)
+    if vol.MON.dtype != numpy.float32:
+        logging.warning("rsv MON was not float32, converting")
+        vol.MON = vol.MON.astype(numpy.float32)
+    volout.create_dataset( "monitor",
+                           (vol.NR[0],vol.NR[1],vol.NR[2]),
+                           vol.MON.dtype,
+                           data = vol.MON,
+                           compression = 'gzip',
+                           compression_opts = 1)
+    for key, value in vol.metadata.iteritems():
+        volout.attrs[key]=value
+    volout.flush()
+    volout.close()
+
+
+def readvol(filename):
+    """
+    Read volume from a file
+    returns an rsv object
+    """
+    volf = h5py.File(filename)
+    try:
+        mon = volf['monitor']
+        sig = volf['signal']
+    except:
+        raise Exception("Your file %s is not an rsv"%(filename))
+    mona= mon[:,:,:]
+    siga= sig[:,:,:]
+    assert mona.shape == siga.shape
+    assert len(mona.shape)==3
+    vol = rsv( mona.shape )
+    for name, value in volf.attrs.iteritems():
+        vol.metadata[name] = value
+    volf.close()
+    return vol
+
+
+
+