Source code for fabio.brukerimage

#!/usr/bin/env python
#coding: utf8


"""

Authors: Henning O. Sorensen & Erik Knudsen
         Center for Fundamental Research: Metal Structures in Four Dimensions
         Risoe National Laboratory
         Frederiksborgvej 399
         DK-4000 Roskilde
         email:erik.knudsen@risoe.dk

Based on: openbruker,readbruker, readbrukerheader functions in the opendata
         module of ImageD11 written by Jon Wright, ESRF, Grenoble, France

Writer by Jérôme Kieffer, ESRF, Grenoble, France

"""
# get ready for python3
from __future__ import with_statement, print_function

__authors__ = ["Henning O. Sorensen" , "Erik Knudsen", "Jon Wright", "Jérôme Kieffer"]
__date__ = "20130502"
__status__ = "development"
__copyright__ = "2007-2009 Risoe National Laboratory; 2010-2013 ESRF"
__licence__ = "GPL"

import numpy, logging, sys
from math import ceil
import os, getpass, time
logger = logging.getLogger("brukerimage")
from .fabioimage import fabioimage
from .fabioutils import pad
from types import StringTypes
if sys.version_info[0] < 3:
    bytes = str

[docs]class brukerimage(fabioimage): """ Read and eventually write ID11 bruker (eg smart6500) images TODO: int32 -> float32 conversion according to the "linear" keyword. This is done and works but we need to check with other program that we are appliing the right formula and not the reciprocal one. """ bpp_to_numpy = {1:numpy.uint8, 2:numpy.uint16, 4:numpy.uint32} # needed if you feel like writing - see ImageD11/scripts/edf2bruker.py SPACER = "\x1a\x04" #this is CTRL-Z CTRL-D HEADERS_KEYS = ["FORMAT", #Frame format. Always “86” or "100" for Bruker-format frames. "VERSION", #Header version #, such as: 1 to 17 (6 is obsolete). "HDRBLKS", #Header size in 512-byte blocks, such as 10 or 15. Determines where the image block begins. "TYPE", #String indicating kind of data in the frame. Used to determine if a spatial correction table was applied to the frame imag "SITE", #Site name "MODEL", #Diffractometer model "USER", #Username "SAMPLE", #Sample ID, "SETNAME", #Basic data set name "RUN", #Run number within the data set, usually starts at 0, but 1 for APEX2. "SAMPNUM", #Specimen number within the data set "TITLE", #User comments (8 lines) "NCOUNTS", #Total frame counts "NOVERFL", #Number of overflows when compression frame. "MINIMUM", #Minimum counts in a pixel (uncompressed value) "MAXIMUM", #Maximum counts in a pixel (uncompressed value) "NONTIME", #Number of on-time events "NLATE", #Number of late events. Always zero for many detectors. "FILENAM", #(Original) frame filename "CREATED", #Date and time of creation "CUMULAT", #Accumulated frame exposure time in seconds "ELAPSDR", #Requested time for last exposure in seconds "ELAPSDA", #Actual time for last exposure in seconds. "OSCILLA", #Nonzero if acquired by oscillation "NSTEPS", #steps or oscillations in this frame "RANGE", #Scan range in decimal degrees (unsigned) "START", #Starting scan angle value, decimal degrees "INCREME", #Scan angle increment between frames (signed) "NUMBER", #Sequence number of this frame in series, usually starts at 0, but 1 for APEX2 "NFRAMES", #Total number of frames in the series "ANGLES", #Diffractometer angles in Eulerian space ( 2T, OM, PH, CH). "NOVER64", #Number of pixels > 64K (actually LinearThreshold value) "NPIXELB", #Number of bytes/pixel, such as 1, 2, or 4. "NROWS", #Number of rasters in frame, such as 512, 1024, 2048, or 4096 "NCOLS", #Number of pixels/raster, such as 512, 1024, 2048 or 4096 "WORDORD", #Order of bytes in word (0=LSB first) "LONGORD", #Order of words in a longword (0=LSW first) "TARGET" , #X-ray target material: Cu, Mo, Ag, Fe, Cr, Co, Ni, W, Mn, or other. "SOURCEK", #X-ray source voltage in kV "SOURCEM", #X-ray source current in mA "FILTER" , #Filter/monochromator setting: Such as: Parallel, graphite, Ni Filter, C Filter, Zr Filter,Cross coupled Goebel Mirrors ... "CELL" , #Unit cell A,B,C,ALPHA,BETA,GAMMA "MATRIX" , #9R Orientation matrix (P3 conventions) "LOWTEMP", #Low temp flag. "TEMP", #set temperature "HITEMP", #Acquired at high temperature "ZOOM" , #Zoom: Xc, Yc, Mag used for HI-STAR detectors: 0.5 0.5 1.0 "CENTER" , #X, Y of direct beam at 2-theta = 0. These are raw center for raw frames and unwarped center for unwarped frames. "DISTANC", #Sample-detector distance, cm (see CmToGrid value) Adds: Sample-detector grid/phosphor distance, cm "TRAILER", #Byte pointer to trailer info "COMPRES", #Compression scheme ID, if any. Such as: NONE, LINEAR (Linear scale, offset for pixel values, typically 1.0, 0.0). "LINEAR", #Linear scale (1.0 0.0 for no change; 0.1 0 for divided by 10...) "PHD" , # Discriminator: Pulse height settings. X100 and X1000 only. Stores CCD phosphor efficiency (first field). "PREAMP" , #Preamp gain setting. X100 and X1000 only. SMART: Stores Roper CCD gain table index value. "CORRECT", #Flood table correction filename, UNKNOWN or LINEAR. "WARPFIL", #Brass plate correction filename, UNKNOWN or LINEAR. Note: A filename here does NOT mean that spatial correction was performed. See TYPE and string “UNWARP” to determine that. "WAVELEN", #Wavelengths (average, a1, a2) "MAXXY", #X,Y pixel # of maximum counts (from lower corner of 0,0) "AXIS", #Scan axis ib Eulerian space (1-4 for 2-theta, omega, phi, chi) (0 =none, 2 = default). "ENDING" , #Actual goniometer angles at end of frame in Eulerian space. "DETPAR" , #Detector position corrections (dX,dY,dDist,Pitch,Roll,Yaw) "LUT", #Recommended display lookup table "DISPLIM", #Recommended display limits "PROGRAM", #Name and version of program writing frame, such as: "ROTATE", #Non zero if acquired by rotation of phi during scan (or oscilate) "BITMASK", #File name of active pixel mask associated with this frame or $NULL "OCTMASK", #Octagon mask parameters to use if BITMASK=$null. Min X, Min X+Y, Min Y, Max X-Y, Max X, Max X+Y, Max Y, Max Y-X. "ESDCELL", #Unit cell parameter standard deviations "DETTYPE", #Detector or CCD chip type (as displayed on CEU). Default is MULTIWIRE but UNKNOWN is advised, can contain PIXPERCM: CMTOGRID: "NEXP", #Number of exposures: 1=single, 2=correlated sum.32 for most ccds, and 64 for 2K ccds. "CCDPARM", #CCD parameters: readnoise, e/ADU, e/photon, bias, full scale "BIS", #Potential full linear scale if rescan and attenuator used. "CHEM", #Chemical formula in CIFTAB string, such as “?” "MORPH", #Crystal morphology in CIFTAB string, such as “?” "CCOLOR", #Crystal color in CIFTAB string, such as “?” "CSIZE", #Crystal dimensions (3 ea) in CIFTAB string, such as “?” "DNSMET", #Density measurement method in CIFTAB string, such as “?” "DARK", #Name of dark current correction or NONE. "AUTORNG", #Auto-ranging: gain, high-speed time, scale, offset, full linear scale Note: If full linear scale is zero, then CCDPARM full scale is the full linear scale (BIS frames). "ZEROADJ", #Goniometer zero corrections (refined in least squares) "XTRANS", #Crystal XYZ translations (refined in least squares) "HKL&XY", #HKL and pixel XY for reciprocal space scan. GADDS only. "AXES2", #Diffractometer setting linear axes (4 ea). (X, Y, Z, Aux) "ENDING2", #Actual goniometer linear axes @ end of frame. (X, Y, Z, Aux) "FILTER2", #Monochromator 2-theta angle and monochromator roll angle. v15: Adds beam tilt angle and attenuator factor. "LEPTOS", # String for LEPTOS. "CFR", #Only in 21CFRPart11 mode, writes the checksum for header and image (2str).] ] def __init__(self, data=None , header=None): fabioimage.__init__(self, data, header) self.__bpp_file = None self.version = 86 self.__headerstring__ = "" def _readheader(self, infile): """ The bruker format uses 80 char lines in key : value format In the first 512*5 bytes of the header there should be a HDRBLKS key, whose value denotes how many 512 byte blocks are in the total header. The header is always n*5*512 bytes, otherwise it wont contain whole key: value pairs """ line = 80 blocksize = 512 nhdrblks = 5 #by default we always read 5 blocks of 512 self.__headerstring__ = infile.read(blocksize * nhdrblks) self.header = {} for i in range(0, nhdrblks * blocksize, line): if self.__headerstring__[i: i + line].find(":") > 0: key, val = self.__headerstring__[i: i + line].split(":", 1) key = key.strip() # remove the whitespace (why?) val = val.strip() if key in self.header: # append lines if key already there self.header[key] = self.header[key] + os.linesep + val else: self.header[key] = val self.header_keys.append(key) # we must have read this in the first 5*512 bytes. nhdrblks = int(self.header['HDRBLKS']) # Now read in the rest of the header blocks, appending self.__headerstring__ += infile.read(blocksize * (nhdrblks - 5)) for i in range(5 * blocksize, nhdrblks * blocksize, line): if self.__headerstring__[i: i + line].find(":") > 0: # as for first 512 bytes of header key, val = self.__headerstring__[i: i + line].split(":", 1) key = key.strip() val = val.strip() if key in self.header: self.header[key] = self.header[key] + os.linesep + val else: self.header[key] = val self.header_keys.append(key) # make a (new) header item called "datastart" self.header['datastart'] = blocksize * nhdrblks #set the image dimensions self.dim1 = int(self.header['NROWS']) self.dim2 = int(self.header['NCOLS'])
[docs] def read(self, fname, frame=None): """ Read in and unpack the pixels (including overflow table """ infile = self._open(fname, "rb") try: self._readheader(infile) except: raise rows = self.dim1 cols = self.dim2 try: # you had to read the Bruker docs to know this! npixelb = int(self.header['NPIXELB']) except Exception: errmsg = "length " + str(len(self.header['NPIXELB'])) + "\n" for byt in self.header['NPIXELB']: errmsg += "char: " + str(byt) + " " + str(ord(byt)) + "\n" logger.warning(errmsg) raise RuntimeError(errmsg) data = numpy.fromstring(infile.read(rows * cols * npixelb), dtype=self.bpp_to_numpy[npixelb]) #handle overflows nov = int(self.header['NOVERFL']) if nov > 0: # Read in the overflows # need at least int32 sized data I guess - can reach 2^21 data = data.astype(numpy.uint32) # 16 character overflows: # 9 characters of intensity # 7 character position for i in range(nov): ovfl = infile.read(16) intensity = int(ovfl[0: 9]) position = int(ovfl[9: 16]) data[position] = intensity infile.close() # Handle Float images ... if "LINEAR" in self.header: try: slope, offset = self.header["LINEAR"].split(None, 1) slope = float(slope) offset = float(offset) except Exception: logger.warning("Error in converting to float data with linear parameter: %s" % self.header["LINEAR"]) self.data = data else: if slope == 1 and offset == 0: self.data = data else: #TODO: check that the formula is OK, not reverted. logger.warning("performing correction with slope=%s, offset=%s (LINEAR=%s)" % (slope, offset, self.header["LINEAR"])) self.data = (data * slope + offset).astype(numpy.float32) else: self.data = data self.data.shape = self.dim1, self.dim2 self.resetvals() self.pilimage = None return self
[docs] def write(self, fname): """ Write a bruker image """ if numpy.issubdtype(self.data.dtype, float): if "LINEAR" in self.header: try: slope, offset = self.header["LINEAR"].split(None, 1) slope = float(slope) offset = float(offset) except Exception: logger.warning("Error in converting to float data with linear parameter: %s" % self.header["LINEAR"]) slope, offset = 1.0, 0.0 else: offset = self.data.min() max_data = self.data.max() max_range = 2 ** 24 - 1 #similar to the mantissa of a float32 if max_data > offset: slope = (max_data - offset) / float(max_range) else: slope = 1.0 tmp_data = numpy.round(((self.data - offset) / slope)).astype(numpy.uint32) self.header["LINEAR"] = "%s %s" % (slope, offset) else: tmp_data = self.data bpp = self.calc_bpp(tmp_data) self.basic_translate(fname) limit = 2 ** (8 * bpp) - 1 data = tmp_data.astype(self.bpp_to_numpy[bpp]) reset = numpy.where(tmp_data >= limit) data[reset] = limit data = data.newbyteorder("<") #Bruker enforces little endian with self._open(fname, "wb") as bruker: bruker.write(self.gen_header()) bruker.write(data.tostring()) bruker.write(self.gen_overflow())
[docs] def calc_bpp(self, data=None, max_entry=4096): """ Calculate the number of byte per pixel to get an optimal overflow table. @return: byte per pixel """ if data is None: data = self.data if self.__bpp_file is None: for i in [1, 2]: overflown = (data >= (2 ** (8 * i) - 1)) if overflown.sum() < max_entry: self.__bpp_file = i break else: self.__bpp_file = 4 return self.__bpp_file
[docs] def gen_header(self): """ Generate headers (with some magic and guesses) @param format can be 86 or 100 """ headers = [] for key in self.HEADERS_KEYS: if key in self.header: value = self.header[key] line = key.ljust(7) + ":" if type(value) in StringTypes: if os.linesep in value: lines = value.split(os.linesep) for i in lines[:-1] : headers.append((line + bytes(i)).ljust(80, " ")) line = key.ljust(7) + ":" line += bytes(lines[-1]) elif len(value) < 72: line += bytes(value) else: for i in range(len(value) // 72): headers.append((line + bytes(value[72 * i:72 * (i + 1)]))) line = key.ljust(7) + ":" line += value[72 * (i + 1):] elif "__len__" in dir(value): f = "\%.%is" % 72 // len(value) - 1 line += " ".join([f % i for i in value]) else: line += bytes(value) headers.append(line.ljust(80, " ")) header = "".join(headers) if len(header) > 512 * self.header["HDRBLKS"]: tmp = ceil(len(header) / 512.0) self.header["HDRBLKS"] = int(ceil(tmp / 5.0) * 5.0) for i in range(len(headers)): if headers[i].startswith("HDRBLKS"): headers[i] = headers.append(("HDRBLKS:%s" % self.header["HDRBLKS"]).ljust(80, " ")) res = pad("".join(headers), self.SPACER + "."*78, 512 * int(self.header["HDRBLKS"])) return res
[docs] def gen_overflow(self): """ Generate an overflow table """ limit = 2 ** (8 * self.calc_bpp()) - 1 flat = self.data.ravel() #flat memory view overflow_pos = numpy.where(flat >= limit)[0] #list of indexes overflow_val = flat[overflow_pos] overflow = "".join(["%09i%07i" % (val, pos) for pos, val in zip(overflow_pos, overflow_val)]) return pad(overflow, ".", 512)
[docs] def basic_translate(self, fname=None): """ Does some basic population of the headers so that the writing is possible """ if not "FORMAT" in self.header: self.header["FORMAT"] = "86" if not "HDRBLKS" in self.header: self.header["HDRBLKS"] = 5 if not "TYPE" in self.header: self.header["TYPE"] = "UNWARPED" if not "USER" in self.header: self.header["USER"] = getpass.getuser() if not "FILENAM" in self.header: self.header["FILENAM"] = "%s" % fname if not "CREATED" in self.header: self.header["CREATED"] = time.ctime() if not "NOVERFL" in self.header: self.header["NOVERFL"] = "0" # if not "NPIXELB" in self.header: self.header["NPIXELB"] = self.calc_bpp() #if not "NROWS" in self.header: self.header["NROWS"] = self.data.shape[0] #if not "NCOLS" in self.header: self.header["NCOLS"] = self.data.shape[1] if not "WORDORD" in self.header: self.header["WORDORD"] = "0" if not "LONGORD" in self.header: self.header["LONGORD"] = "0"
[docs]def test(): """ a testcase """ import sys, time img = brukerimage() start = time.clock() for filename in sys.argv[1:]: img.read(filename) res = img.toPIL16() img.rebin(2, 2) print(filename + (": max=%d, min=%d, mean=%.2e, stddev=%.2e") % ( img.getmax(), img.getmin(), img.getmean(), img.getstddev())) print('integrated intensity (%d %d %d %d) =%.3f' % ( 10, 20, 20, 40, img.integrate_area((10, 20, 20, 40)))) end = time.clock() print (end - start)
if __name__ == '__main__': test()