From 4093109a9667e3ba8c5307e4031a688c1a99ccc5 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Mon, 1 Jul 2024 14:28:56 +1000 Subject: [PATCH 1/5] First attempt at Nonius Kappa CCD format. --- src/dxtbx/format/FormatNoniusKappaCCD.py | 262 +++++++++++++++++++++++ 1 file changed, 262 insertions(+) create mode 100644 src/dxtbx/format/FormatNoniusKappaCCD.py diff --git a/src/dxtbx/format/FormatNoniusKappaCCD.py b/src/dxtbx/format/FormatNoniusKappaCCD.py new file mode 100644 index 000000000..7c5a5ab2f --- /dev/null +++ b/src/dxtbx/format/FormatNoniusKappaCCD.py @@ -0,0 +1,262 @@ +from __future__ import annotations + +import io +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 read_uint16_bs +from dxtbx.format.Format import Format + +MINIMUM_KEYS = [ + # 'ByteOrder', Assume little by default + "Data type", + "X dimension", + "Y dimension", + "Number of readouts", +] + + +class FormatNoniusKappaCCD(Format): + """This format produced by 1990s era Nonius Kappa CCD instruments. + Consists of a header containing lines of the form XXX = YYY, followed by an image. + We reuse Fabio code. + """ + + @staticmethod + def read_header_lines(image_path): + """ + Adapted from Fabio KcdImage _readheader + """ + + hdr_lines = [] + end_of_headers = False + with Format.open_file(image_path, "rb") as f: + while not end_of_headers: + one_line = f.readline() + 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: + 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): + """ + Try to find some 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 + 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"]) + + # Kappa axis points into the kappa head as clockwise rotation when + # viewed from above. The datum position has the kappa head under + # the incoming beam. All primary axes rotate clockwise when viewed + # from above, so the kappa axis has positive x, z components, + # corresponding to "+z" when the calculation in in make_kappa_goniometer + # is reviewed. + + 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. CHECK""" + + 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"]), + ] + + self._detector_factory.two_theta( + sensor="CCD", + distance=float(self.headers["Dx start"]), + beam_centre=[image_size[0] / 2.0, image_size[1] / 2], + fast_direction="+y", # Increasing to the right + slow_direction="+x", # Down, same as rotation axis + two_theta_direction=[1, 0, 0], # Same sense as omega and phi + two_theta_angle=float(self.headers["Theta start"]) * 2.0, + pixel_size=pixel_size, + image_size=image_size, + identifier=self.headers["Detector ID"], + ) + + 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. Polarisation is that + obtained after single reflection from graphite 002 + monochromator.""" + a1 = float(self.headers["Alpha1"]) + a2 = float(self.headers["Alpha2"]) + wt = float(self.headers["Alpha1/Alpha2 ratio"]) + wavelength = wt * a1 + a2 + direction = [0.0, 0.0, 1.0] # imgCIF standard + polarisation, pol_dir = self.get_beam_polarisation() + return self._beam_factory.complex( + sample_to_source=direction, + wavelength=wavelength, + polarization=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.header["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. Hard-coded angles for MO and CU + 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 + + pol_frac = (1 + math.cos(two_theta * math.pi / 180) ** 2) / 2 + pol_dir = [1.0, 0.0, 0.0] # To be confirmed later. + return pol_frac, pol_dir + + def get_raw_data(self): + # Compute image size: always little-endian UInt16 = 2 bytes + + 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 + + # Now seek to beginning counting from the end + + with self.open_file(self._image_file, "rb") as infile: + try: + infile.seek(-expected_size, io.SEEK_END) + except Exception: + print( + "Warning: Seeking from end not implemented for file %s" + % self._image_file + ) + if hasattr(infile, "measure_size"): + fileSize = infile.measure_size() + elif hasattr(infile, "size"): + fileSize = infile.size + elif hasattr(infile, "getSize"): + fileSize = infile.getSize() + else: + print("Unable to guess the file-size of %s" % self._image_file) + fileSize = os.stat(self._image_file)[6] + infile.seek(fileSize - expected_size - infile.tell(), 1) + + # Read data into array. Data are little endian based on Fabio approach + # Data may be double-read, presumably in case some CCD bins had something + # left after the first read. + + for i in range(nbReadOut): + raw_data = read_uint16_bs(streambuf(infile), dim1 * dim2) + if i == 0: + data = raw_data + else: + data += raw_data + + data.reshape(flex.grid(dim2, dim1)) # Not sure about correct order of dims + return data + + +if __name__ == "__main__": + for arg in sys.argv[1:]: + print(FormatNoniusKappaCCD.understand(arg)) From f415890faafcf3fb3d6fed42f8b970d6bfa10f25 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Wed, 3 Jul 2024 17:50:58 +1000 Subject: [PATCH 2/5] Added and adjusted geometry to allow dials processing. KappaCCD files now display correctly and indexing for individual scans is correct. Kappa movement is not correct. --- src/dxtbx/format/FormatNoniusKappaCCD.py | 87 +++++++++++++----------- 1 file changed, 46 insertions(+), 41 deletions(-) diff --git a/src/dxtbx/format/FormatNoniusKappaCCD.py b/src/dxtbx/format/FormatNoniusKappaCCD.py index 7c5a5ab2f..10cbcb0eb 100644 --- a/src/dxtbx/format/FormatNoniusKappaCCD.py +++ b/src/dxtbx/format/FormatNoniusKappaCCD.py @@ -1,6 +1,5 @@ from __future__ import annotations -import io import math import os import re @@ -11,7 +10,7 @@ from scitbx.array_family import flex from dxtbx import IncorrectFormatError -from dxtbx.ext import read_uint16_bs +from dxtbx.ext import is_big_endian, read_uint16, read_uint16_bs from dxtbx.format.Format import Format MINIMUM_KEYS = [ @@ -37,9 +36,11 @@ def read_header_lines(image_path): 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: @@ -50,10 +51,11 @@ def read_header_lines(image_path): if not end_of_headers: if one_line.strip() == "Binned mode": one_line = "Mode = Binned" - if "=" not in one_line: + 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 @@ -90,6 +92,8 @@ def understand(image_file): 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): @@ -132,12 +136,13 @@ def _detector(self): account for binning, which doubles the size. CHECK""" pix_fac = 0.001 - # if self.headers["Mode"] == "Binned": - # pix_fac = 0.002 + 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, + # float(self.headers["pixel X-size (um)"]) * pix_fac #For testing + float(self.headers["pixel Y-size (um)"]) * pix_fac, # True ] image_size = [ @@ -145,17 +150,27 @@ def _detector(self): float(self.headers["Y dimension"]), ] - self._detector_factory.two_theta( + 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=[image_size[0] / 2.0, image_size[1] / 2], + beam_centre=beam_centre, fast_direction="+y", # Increasing to the right slow_direction="+x", # Down, same as rotation axis - two_theta_direction=[1, 0, 0], # Same sense as omega and phi + 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): @@ -167,13 +182,13 @@ def _beam(self): a1 = float(self.headers["Alpha1"]) a2 = float(self.headers["Alpha2"]) wt = float(self.headers["Alpha1/Alpha2 ratio"]) - wavelength = wt * a1 + a2 + wavelength = (wt * a1 + a2) / (wt + 1.0) direction = [0.0, 0.0, 1.0] # imgCIF standard polarisation, pol_dir = self.get_beam_polarisation() return self._beam_factory.complex( sample_to_source=direction, wavelength=wavelength, - polarization=pol_dir, + polarization_plane_normal=pol_dir, polarization_fraction=polarisation, ) @@ -190,7 +205,7 @@ def _scan(self): osc_range = float(self.headers["%s scan range" % rotax]) exposure_time = float(self.headers["Exposure time"]) - simpletime = self.header["Date"] + simpletime = self.headers["Date"] epoch = time.mktime(time.strptime(simpletime, "%a %m/%d/%y %I:%M:%S %p")) return self._scan_factory.single_file( @@ -209,7 +224,8 @@ def get_beam_polarisation(self): elif self.headers["Target material"] == "AG": two_theta = 9.62 # Calculated - pol_frac = (1 + math.cos(two_theta * math.pi / 180) ** 2) / 2 + # Check pol frac against Azaroff paper + pol_frac = 0.5 * (1 + math.cos(two_theta * math.pi / 180) ** 2) / 2 pol_dir = [1.0, 0.0, 0.0] # To be confirmed later. return pol_frac, pol_dir @@ -224,36 +240,25 @@ def get_raw_data(self): # Now seek to beginning counting from the end with self.open_file(self._image_file, "rb") as infile: - try: - infile.seek(-expected_size, io.SEEK_END) - except Exception: - print( - "Warning: Seeking from end not implemented for file %s" - % self._image_file - ) - if hasattr(infile, "measure_size"): - fileSize = infile.measure_size() - elif hasattr(infile, "size"): - fileSize = infile.size - elif hasattr(infile, "getSize"): - fileSize = infile.getSize() + fileSize = os.stat(self._image_file)[6] + infile.seek(fileSize - expected_size, os.SEEK_SET) + + # Read data into array. Data are little endian based on Fabio approach + # Data may be double-read, presumably in case some CCD bins had something + # left after the first read. + + for i in range(nbReadOut): + if is_big_endian(): + raw_data = read_uint16_bs(streambuf(infile), dim1 * dim2) else: - print("Unable to guess the file-size of %s" % self._image_file) - fileSize = os.stat(self._image_file)[6] - infile.seek(fileSize - expected_size - infile.tell(), 1) - - # Read data into array. Data are little endian based on Fabio approach - # Data may be double-read, presumably in case some CCD bins had something - # left after the first read. - - for i in range(nbReadOut): - raw_data = read_uint16_bs(streambuf(infile), dim1 * dim2) - if i == 0: - data = raw_data - else: - data += raw_data + raw_data = read_uint16(streambuf(infile), dim1 * dim2) + + if i == 0: + data = raw_data + else: + data += raw_data - data.reshape(flex.grid(dim2, dim1)) # Not sure about correct order of dims + data.reshape(flex.grid(dim2, dim1)) # These dims work out return data From 9d8ecbf179ce6e4b0a36532927f562b6a1d9fed0 Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Thu, 4 Jul 2024 14:11:27 +1000 Subject: [PATCH 3/5] Fix kappa axis direction and polarisation calculation, clean up comments. --- src/dxtbx/format/FormatNoniusKappaCCD.py | 68 +++++++++--------------- 1 file changed, 25 insertions(+), 43 deletions(-) diff --git a/src/dxtbx/format/FormatNoniusKappaCCD.py b/src/dxtbx/format/FormatNoniusKappaCCD.py index 10cbcb0eb..9de12dbb7 100644 --- a/src/dxtbx/format/FormatNoniusKappaCCD.py +++ b/src/dxtbx/format/FormatNoniusKappaCCD.py @@ -13,20 +13,9 @@ from dxtbx.ext import is_big_endian, read_uint16, read_uint16_bs from dxtbx.format.Format import Format -MINIMUM_KEYS = [ - # 'ByteOrder', Assume little by default - "Data type", - "X dimension", - "Y dimension", - "Number of readouts", -] - class FormatNoniusKappaCCD(Format): - """This format produced by 1990s era Nonius Kappa CCD instruments. - Consists of a header containing lines of the form XXX = YYY, followed by an image. - We reuse Fabio code. - """ + """A class for reading files produced by 1990s era Nonius Kappa CCD instruments.""" @staticmethod def read_header_lines(image_path): @@ -78,7 +67,7 @@ def parse_header(header_lines): @staticmethod def understand(image_file): """ - Try to find some characteristic character sequences in the + Look for characteristic character sequences in the header. """ try: @@ -112,15 +101,7 @@ def _goniometer(self): omega = float(self.headers["Omega start"]) kappa = float(self.headers["Kappa start"]) phi = float(self.headers["Phi start"]) - - # Kappa axis points into the kappa head as clockwise rotation when - # viewed from above. The datum position has the kappa head under - # the incoming beam. All primary axes rotate clockwise when viewed - # from above, so the kappa axis has positive x, z components, - # corresponding to "+z" when the calculation in in make_kappa_goniometer - # is reviewed. - - direction = "+z" # + direction = "-z" # if "Phi scan range" in self.headers: scan_axis = "phi" @@ -133,7 +114,7 @@ def _goniometer(self): def _detector(self): """It appears that pixel size reported in the header does not - account for binning, which doubles the size. CHECK""" + account for binning, which doubles the size in both dimensions.""" pix_fac = 0.001 if self.headers["Mode"] == "Binned": @@ -141,8 +122,7 @@ def _detector(self): pixel_size = [ float(self.headers["pixel X-size (um)"]) * pix_fac, - # float(self.headers["pixel X-size (um)"]) * pix_fac #For testing - float(self.headers["pixel Y-size (um)"]) * pix_fac, # True + float(self.headers["pixel Y-size (um)"]) * pix_fac, ] image_size = [ @@ -176,14 +156,13 @@ def _detector(self): 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. Polarisation is that - obtained after single reflection from graphite 002 - monochromator.""" + 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] # imgCIF standard + direction = [0.0, 0.0, 1.0] polarisation, pol_dir = self.get_beam_polarisation() return self._beam_factory.complex( sample_to_source=direction, @@ -213,10 +192,10 @@ def _scan(self): ) def get_beam_polarisation(self): - """Polarisation for single reflection off graphite - 002 monochromator. Hard-coded angles for MO and CU - targets - """ + """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": @@ -224,29 +203,32 @@ def get_beam_polarisation(self): elif self.headers["Target material"] == "AG": two_theta = 9.62 # Calculated - # Check pol frac against Azaroff paper - pol_frac = 0.5 * (1 + math.cos(two_theta * math.pi / 180) ** 2) / 2 - pol_dir = [1.0, 0.0, 0.0] # To be confirmed later. + # 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): - # Compute image size: always little-endian UInt16 = 2 bytes + """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 - # Now seek to beginning counting from the end + # 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) - # Read data into array. Data are little endian based on Fabio approach - # Data may be double-read, presumably in case some CCD bins had something - # left after the first read. - for i in range(nbReadOut): if is_big_endian(): raw_data = read_uint16_bs(streambuf(infile), dim1 * dim2) @@ -258,7 +240,7 @@ def get_raw_data(self): else: data += raw_data - data.reshape(flex.grid(dim2, dim1)) # These dims work out + data.reshape(flex.grid(dim2, dim1)) return data From 7950b3669cc5d3579842cc63cb31557eeefe14df Mon Sep 17 00:00:00 2001 From: "James.Hester" Date: Thu, 4 Jul 2024 14:18:12 +1000 Subject: [PATCH 4/5] Add newsfragment for kappa CCD update. --- newsfragments/xxx.feature | 1 + 1 file changed, 1 insertion(+) create mode 100644 newsfragments/xxx.feature diff --git a/newsfragments/xxx.feature b/newsfragments/xxx.feature new file mode 100644 index 000000000..5c3023f73 --- /dev/null +++ b/newsfragments/xxx.feature @@ -0,0 +1 @@ +Add Nonius KappaCCD format. From f6e027157e5c4cfb8864f2448084c52dd56051db Mon Sep 17 00:00:00 2001 From: DiamondLightSource-build-server Date: Thu, 4 Jul 2024 04:35:28 +0000 Subject: [PATCH 5/5] Rename newsfragments/xxx.feature to newsfragments/741.feature --- newsfragments/{xxx.feature => 741.feature} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename newsfragments/{xxx.feature => 741.feature} (100%) diff --git a/newsfragments/xxx.feature b/newsfragments/741.feature similarity index 100% rename from newsfragments/xxx.feature rename to newsfragments/741.feature