FormatROD: Use a multiaxis goniometer model to allow simultaneous

indexing of multiple scans

src/dxtbx/format/FormatROD.py

@@ -23,6 +23,7 @@
 import numpy as np
 
 from scitbx.array_family import flex
+from scitbx.math import r3_rotation_axis_and_angle_as_matrix
 
 from dxtbx.format.Format import Format
 
@@ -153,7 +154,7 @@ def _read_binary_header(
             f.seek(offset + general_nbytes + special_nbytes + 640)
             # detector rotation in degrees along e1, e2, e3
             detector_rotns = struct.unpack("<ddd", f.read(24))
-            # FIXME: direct beam position when all angles are zero?
+            # direct beam position when all angles are zero?
             origin_px_x, origin_px_y = struct.unpack("<dd", f.read(16))
             # alpha and beta are angles between KAPPA(=CHI) and THETA, and e3.
             angles_in_deg = struct.unpack(
@@ -198,6 +199,25 @@ def _start(self):
         self._txt_header = self._read_ascii_header(self._image_file)
         self._bin_header = self._read_binary_header(self._image_file)
 
+        self._gonio_start_angles = (
+            np.array(self._bin_header["start_angles_steps"])
+            * np.array(self._bin_header["step_to_rad"])
+            / np.pi
+            * 180
+        )
+        self._gonio_end_angles = (
+            np.array(self._bin_header["end_angles_steps"])
+            * np.array(self._bin_header["step_to_rad"])
+            / np.pi
+            * 180
+        )
+
+        self._scan_axis = -1
+        for axis in [0, 3]:  # 0 - OMEGA, 3 - PHI; the default is omega scan
+            if self._gonio_start_angles[axis] != self._gonio_end_angles[axis]:
+                self._scan_axis = axis
+                break
+
         return
 
     # Rigaku/Oxford coordinate system:
@@ -229,12 +249,36 @@ def get_goniometer(self, index=None):
         return Format.get_goniometer(self)
 
     def _goniometer(self):
+        if self._scan_axis == -1:
+            raise NotImplementedError("Still shots not implemented yet.")
+        elif self._scan_axis == 0:  # OMEGA
+            dxtbx_scan_axis = 2
+        elif self._scan_axis == 3:  # PHI
+            dxtbx_scan_axis = 0
+        else:
+            pass  # should not happen
+
         # FIXME: sometimes XDS.INP generated by CrysAlisPro has
         #  tiny deviations from (0, 1, 0). I don't know how to calculate it.
-        direction = (0.0, -1.0, 0.0)
+        alpha_rad = self._bin_header["angles_in_deg"][0] * np.pi / 180
+        axes = flex.vec3_double(
+            ((0, -1, 0), (0, -np.cos(alpha_rad), np.sin(alpha_rad)), (0, -1, 0))
+        )
 
-        # FIXME: represent kappa axis in fixed_rotation
-        return self._goniometer_factory.known_axis(direction)
+        # angles[self._scan_axis] is not used anyway
+        # angles are in degrees!
+        angles = flex.double(
+            (
+                self._gonio_start_angles[3],
+                self._gonio_start_angles[2],
+                self._gonio_start_angles[0],
+            )
+        )
+        names = flex.std_string(("PHI", "KAPPA=CHI", "OMEGA"))
+
+        return self._goniometer_factory.make_multi_axis_goniometer(
+            axes, angles, names, dxtbx_scan_axis
+        )
 
     def get_beam(self, index=None):
         return Format.get_beam(self)
@@ -251,53 +295,27 @@ def get_detector(self, index=None):
         return Format.get_detector(self)
 
     def _detector(self):
-        gonio_angles = np.array(self._bin_header["start_angles_steps"]) * np.array(
-            self._bin_header["step_to_rad"]
-        )
+        theta_rad = self._gonio_start_angles[1] / 180 * np.pi
         detector_rotns_rad = np.array(self._bin_header["detector_rotns"]) / 180 * np.pi
 
         # I don't know why only rot_e1 is clockwise but this matches
         # DIRECTION_OF_DETECTOR_X-AXIS/Y-AXIS in XDS.INP from CrysAlisPro.
         # Note that XDS.INP's directions of Y and Z are opposite from ours.
-        rot_e1 = np.array(  # clockwise along e1 = Z
-            [
-                np.cos(detector_rotns_rad[0]),
-                -np.sin(detector_rotns_rad[0]),
-                0,
-                np.sin(detector_rotns_rad[0]),
-                np.cos(detector_rotns_rad[0]),
-                0,
-                0,
-                0,
-                1,
-            ]
-        ).reshape(3, 3)
-        rot_e2 = np.array(  # ANTI-clockwise along e2 = X
-            [
-                1,
-                0,
-                0,
-                0,
-                np.cos(detector_rotns_rad[1]),
-                np.sin(detector_rotns_rad[1]),
-                0,
-                -np.sin(detector_rotns_rad[1]),
-                np.cos(detector_rotns_rad[1]),
-            ]
-        ).reshape(3, 3)
-        rot_theta = np.array(  # ANTI-clockwise along e3 = Y
-            [
-                np.cos(gonio_angles[1]),
-                0,
-                -np.sin(gonio_angles[1]),
-                0,
-                1,
-                0,
-                np.sin(gonio_angles[1]),
-                0,
-                np.cos(gonio_angles[1]),
-            ]
-        ).reshape(3, 3)
+        rot_e1 = np.array(
+            r3_rotation_axis_and_angle_as_matrix([0, 0, 1], detector_rotns_rad[0])
+        ).reshape(
+            3, 3
+        )  # clockwise along e1 = Z
+        rot_e2 = np.array(
+            r3_rotation_axis_and_angle_as_matrix([-1, 0, 0], detector_rotns_rad[1])
+        ).reshape(
+            3, 3
+        )  # ANTI-clockwise along e2 = X
+        rot_theta = np.array(
+            r3_rotation_axis_and_angle_as_matrix([0, -1, 0], theta_rad)
+        ).reshape(
+            3, 3
+        )  # ANTI-clockwise along e3 = Y
         detector_axes = rot_theta.dot(rot_e2.dot(rot_e1))
 
         pixel_size_x = self._bin_header["real_px_size_x"]
@@ -309,9 +327,6 @@ def _detector(self):
                 -self._bin_header["distance_mm"],
             ]
         )
-
-        # XDS exporter in IPR's CrysAlisPro seems broken. It writes the same ORGX/Y
-        # regardless of the theta angles.
         origin = detector_axes.dot(origin_at_zero)
 
         detector = self._detector_factory.make_detector(
@@ -332,23 +347,8 @@ def get_scan(self, index=None):
     def _scan(self):
         """Return the scan information for this image."""
 
-        for axis in [0, 3]:  # 0 - OMEGA, 3 - PHI; the default is omega scan
-            start_angle = (
-                self._bin_header["start_angles_steps"][axis]
-                * self._bin_header["step_to_rad"][axis]
-                / np.pi
-                * 180
-            )
-            end_angle = (
-                self._bin_header["end_angles_steps"][axis]
-                * self._bin_header["step_to_rad"][axis]
-                / np.pi
-                * 180
-            )
-            if start_angle != end_angle:
-                if axis == 3:  # PHI scan
-                    pass  # TODO: FIXME: the goniometer axis is not 0,-1,0 in this case!
-                break
+        start_angle = self._gonio_start_angles[self._scan_axis]
+        end_angle = self._gonio_end_angles[self._scan_axis]
 
         return self._scan_factory.single_file(
             filename=self._image_file,
@@ -481,18 +481,8 @@ def decode_TY6_oneline(self, linedata, w):
             print(
                 "Starting angles in degrees (OMEGA, THETA, KAPPA=CHI, PHI, OMEGA PRIME, THETA PRIME)"
             )
-            print(
-                np.array(reader._bin_header["start_angles_steps"])
-                * np.array(reader._bin_header["step_to_rad"])
-                * 180
-                / np.pi
-            )
+            print(reader._gonio_start_angles)
             print("Ending angles in degrees")
-            print(
-                np.array(reader._bin_header["end_angles_steps"])
-                * np.array(reader._bin_header["step_to_rad"])
-                * 180
-                / np.pi
-            )
+            print(reader._gonio_end_angles)
         else:
             print("Unsupported format")