Add Nonius Kappa CCD format module

newsfragments/741.feature

@@ -0,0 +1 @@
+Add Nonius KappaCCD format.

src/dxtbx/format/FormatNoniusKappaCCD.py

@@ -0,0 +1,249 @@
+from __future__ import annotations
+
+import math
+import os
+import re
+import sys
+import time
+
+from boost_adaptbx.boost.python import streambuf
+from scitbx.array_family import flex
+
+from dxtbx import IncorrectFormatError
+from dxtbx.ext import is_big_endian, read_uint16, read_uint16_bs
+from dxtbx.format.Format import Format
+
+
+class FormatNoniusKappaCCD(Format):
+ """A class for reading files produced by 1990s era Nonius Kappa CCD instruments."""
+
+ @staticmethod
+ def read_header_lines(image_path):
+ """
+ Adapted from Fabio KcdImage _readheader
+ """
+
+ hdr_lines = []
+ end_of_headers = False
+ linect = 0
+ with Format.open_file(image_path, "rb") as f:
+ while not end_of_headers:
+ one_line = f.readline()
+ linect += 1
+ try:
+ one_line = one_line.decode("ASCII")
+ except UnicodeDecodeError:
+ end_of_headers = True
+ else:
+ if len(one_line) > 100:
+ end_of_headers = True
+ if not end_of_headers:
+ if one_line.strip() == "Binned mode":
+ one_line = "Mode = Binned"
+ if "=" not in one_line and linect > 1:
+ end_of_headers = True
+ if not end_of_headers:
+ hdr_lines.append(one_line)
+
+ return hdr_lines
+
+ @staticmethod
+ def parse_header(header_lines):
+ header_dic = {}
+
+ for l in header_lines:
+ separator = l.find("=")
+ if separator == -1:
+ continue
+ k = l[:separator].strip()
+ v = l[separator + 1 :].strip()
+ if k in header_dic:
+ header_dic[k] = header_dic[k] + "\n" + v
+ else:
+ header_dic[k] = v
+
+ return header_dic
+
+ @staticmethod
+ def understand(image_file):
+ """
+ Look for characteristic character sequences in the
+ header.
+ """
+ try:
+ with Format.open_file(image_file, "rb") as fh:
+ tag = fh.read(1024).decode("latin-1", "replace")
+ except OSError:
+ return False
+ if tag[0:2] != "No":
+ return False
+ if re.search("^Kappa-support angle", tag, re.MULTILINE) is None:
+ return False
+ if re.search("^Binned mode", tag, re.MULTILINE) is None:
+ return False
+ if re.search("^Shutter = Closed", tag, re.MULTILINE) is not None:
+ return False # ignore dark current measurement
+ return True
+
+ def __init__(self, image_file, **kwargs):
+ """Initialise the image structure from the given file."""
+
+ if not self.understand(image_file):
+ raise IncorrectFormatError(self, image_file)
+ super().__init__(str(image_file), **kwargs)
+
+ def _start(self):
+ self.headers = self.parse_header(self.read_header_lines(self._image_file))
+
+ def _goniometer(self):
+ """Return model for a kappa goniometer"""
+ alpha = float(self.headers["Kappa-support angle"])
+ omega = float(self.headers["Omega start"])
+ kappa = float(self.headers["Kappa start"])
+ phi = float(self.headers["Phi start"])
+ direction = "-z" #
+
+ if "Phi scan range" in self.headers:
+ scan_axis = "phi"
+ else:
+ scan_axis = "omega"
+
+ return self._goniometer_factory.make_kappa_goniometer(
+ alpha, omega, kappa, phi, direction, scan_axis
+ )
+
+ def _detector(self):
+ """It appears that pixel size reported in the header does not
+ account for binning, which doubles the size in both dimensions."""
+
+ pix_fac = 0.001
+ if self.headers["Mode"] == "Binned":
+ pix_fac = 0.002
+
+ pixel_size = [
+ float(self.headers["pixel X-size (um)"]) * pix_fac,
+ float(self.headers["pixel Y-size (um)"]) * pix_fac,
+ ]
+
+ image_size = [
+ float(self.headers["X dimension"]),
+ float(self.headers["Y dimension"]),
+ ]
+
+ beam_centre = (
+ image_size[0] * pixel_size[0] * 0.5,
+ image_size[1] * pixel_size[1] * 0.5,
+ )
+
+ gain = self.headers.get("Detector gain (ADU/count)", "1.0")
+ gain = float(gain)
+
+ return self._detector_factory.two_theta(
+ sensor="CCD",
+ distance=float(self.headers["Dx start"]),
+ beam_centre=beam_centre,
+ fast_direction="+y", # Increasing to the right
+ slow_direction="+x", # Down, same as rotation axis
+ two_theta_direction="+x", # Same sense as omega and phi
+ two_theta_angle=float(self.headers["Theta start"]) * 2.0,
+ pixel_size=pixel_size,
+ image_size=image_size,
+ gain=gain,
+ identifier=self.headers["Detector ID"],
+ trusted_range=(0.0, 131070.0), # May readout UInt16 twice
+ )
+
+ def _beam(self):
+ """Return a simple model for a lab-based beam. As dxtbx
+ has no laboratory target model, take the weighted values
+ of the alpha1/alpha2 wavelengths."""
+
+ a1 = float(self.headers["Alpha1"])
+ a2 = float(self.headers["Alpha2"])
+ wt = float(self.headers["Alpha1/Alpha2 ratio"])
+ wavelength = (wt * a1 + a2) / (wt + 1.0)
+ direction = [0.0, 0.0, 1.0]
+ polarisation, pol_dir = self.get_beam_polarisation()
+ return self._beam_factory.complex(
+ sample_to_source=direction,
+ wavelength=wavelength,
+ polarization_plane_normal=pol_dir,
+ polarization_fraction=polarisation,
+ )
+
+ def _scan(self):
+ """All scan dates will be in the range 1969-2068 as per
+ Python strptime. Please retire your CCD before 2069."""
+
+ if "Phi scan range" in self.headers:
+ rotax = "Phi"
+ else:
+ rotax = "Omega"
+
+ osc_start = float(self.headers["%s start" % rotax])
+ osc_range = float(self.headers["%s scan range" % rotax])
+ exposure_time = float(self.headers["Exposure time"])
+
+ simpletime = self.headers["Date"]
+ epoch = time.mktime(time.strptime(simpletime, "%a %m/%d/%y %I:%M:%S %p"))
+
+ return self._scan_factory.single_file(
+ self._image_file, exposure_time, osc_start, osc_range, epoch
+ )
+
+ def get_beam_polarisation(self):
+ """Polarisation for single reflection off graphite 002
+ monochromator as per manual. Hard-coded angles for Mo, Cu and
+ Ag targets."""
+
+ if self.headers["Target material"] == "MO":
+ two_theta = 12.2 # From manual
+ elif self.headers["Target material"] == "CU":
+ two_theta = 26.6 # From manual
+ elif self.headers["Target material"] == "AG":
+ two_theta = 9.62 # Calculated
+
+ # Assume an ideally imperfect monochromator, with the axis
+ # of rotation of the monochromator in the horizontal direction
+ # in a laboratory setup.
+
+ pol_frac = 1.0 / (1 + math.cos(two_theta * math.pi / 180) ** 2)
+ pol_dir = [0.0, 1.0, 0.0]
+ return pol_frac, pol_dir
+
+ def get_raw_data(self):
+ """Return raw data from a Nonius CCD frame."""
+
+ # Frame file contains a series of readouts, each pixel is
+ # an unsigned little-endian 16-bit integer
+
+ dim1 = int(self.headers["X dimension"])
+ dim2 = int(self.headers["Y dimension"])
+ nbReadOut = int(self.headers["Number of readouts"])
+ expected_size = dim1 * dim2 * 2 * nbReadOut
+
+ # Not clear if the image offset is present in the header,
+ # therefore we seek from the end of the file.
+
+ with self.open_file(self._image_file, "rb") as infile:
+ fileSize = os.stat(self._image_file)[6]
+ infile.seek(fileSize - expected_size, os.SEEK_SET)
+
+ for i in range(nbReadOut):
+ if is_big_endian():
+ raw_data = read_uint16_bs(streambuf(infile), dim1 * dim2)
+ else:
+ raw_data = read_uint16(streambuf(infile), dim1 * dim2)
+
+ if i == 0:
+ data = raw_data
+ else:
+ data += raw_data
+
+ data.reshape(flex.grid(dim2, dim1))
+ return data
+
+
+if __name__ == "__main__":
+ for arg in sys.argv[1:]:
+ print(FormatNoniusKappaCCD.understand(arg))