Break the DIALS/xfel circular dependency (#627)

DIALS and dxtbx have had dependencies on cctbx_project/xfel, and vice versa. This PR is one of 8 that will break the cyclic dependency.

Notably for dxtbx, this specifically removes:
pytest.importorskip("xfel")

For more detail, see cctbx/cctbx_project#872 and dials/dials#2404

Commits:
* Time code moved from xfel to serialtbx
* Move in radial_average c++ code from xfel
* dpack moved to serialtbx
* Jiffy function iterate_detector_at_level moved from xfel to serialtbx
* Utility functions and constants from cspad_tbx moved to serialtbx
* basis_from_geo and basis moved into serialtbx
* Use dxtbx's cbf_wrapper by adding in the only function in the xfel subclass of cbf_wrapper
* Move in add_frame_specific_cbf_tables
* xfel.util.jungfrau moved to serialtbx.detector
* cspad constants and data reading moved to serialtbx
* Rayonix toolbox methods and constants moved to serialtbx
* Legacy image pickle metrology code moved to serialtbx.detector.legacy_metrology
* Remove xfel import tests
* Better handle if SIT_* variables are missing for psana (IE don't crash dxtbx in the import step for FormatXTC)
* Fix leftover conflict
* Remove xfel reference that snuck in with nxmx_writer
* Use CCTBX nightly until serialtbx hits stable cctbx
* Remove importer skip for xfel from test_nxmx_writer

---------

Co-authored-by: Nicholas Devenish <[email protected]>

.azure-pipelines/ci-conda-env.txt

@@ -2,7 +2,7 @@ conda-forge::boost
 conda-forge::boost-cpp
 conda-forge::bzip2
 conda-forge::c-compiler<1.5
-conda-forge::cctbx-base==2023.7
+cctbx-nightly::cctbx-base
 conda-forge::conda
 conda-forge::cxx-compiler<1.5
 conda-forge::python-dateutil

newsfragments/627.feature

@@ -0,0 +1 @@
+Remove circular dependencies between dxtbx and ``cctbx.xfel`` by using the new ``serialtbx``.

src/dxtbx/boost_python/ext.cpp

@@ -2,6 +2,8 @@
 #include <boost/python/def.hpp>
 #include <scitbx/array_family/shared.h>
 #include <scitbx/array_family/flex_types.h>
+#include <scitbx/constants.h>
+#include <scitbx/vec2.h>
 #include <boost_adaptbx/python_streambuf.h>
 #include <cstdint>
 #include <cctype>
@@ -193,6 +195,41 @@ namespace dxtbx { namespace boost_python {
     return PyBytes_FromStringAndSize(&*packed.begin(), packed.size());
   }
 
+  double distance_between_points(scitbx::vec2<int> const& a, scitbx::vec2<int> const& b) {
+    return std::sqrt((std::pow(double(b[0]-a[0]),2)+std::pow(double(b[1]-a[1]),2)));
+  }
+
+  void radial_average(scitbx::af::versa<double, scitbx::af::flex_grid<> > & data,
+  scitbx::af::versa<bool, scitbx::af::flex_grid<> > & mask,
+  scitbx::vec2<int> const& beam_center,
+  scitbx::af::shared<double> sums,
+  scitbx::af::shared<double> sums_sq,
+  scitbx::af::shared<int> counts,
+  double pixel_size, double distance,
+  scitbx::vec2<int> const& upper_left,
+  scitbx::vec2<int> const& lower_right) {
+    std::size_t extent = sums.size();
+    double extent_in_mm = extent * pixel_size;
+    double extent_two_theta = std::atan(extent_in_mm/distance)*180/scitbx::constants::pi;
+
+    for(std::size_t y = upper_left[1]; y < lower_right[1]; y++) {
+      for(std::size_t x = upper_left[0]; x < lower_right[0]; x++) {
+        double val = data(x,y);
+        if(val > 0 && mask(x,y)) {
+          scitbx::vec2<int> point((int)x,(int)y);
+          double d_in_mm = distance_between_points(point,beam_center) * pixel_size;
+          double twotheta = std::atan(d_in_mm/distance)*180/scitbx::constants::pi;
+          std::size_t bin = (std::size_t)std::floor(twotheta*extent/extent_two_theta);
+          if (bin >= extent)
+            continue;
+          sums[bin] += val;
+          sums_sq[bin] += val*val;
+          counts[bin]++;
+        }
+      }
+    }
+  }
+
   // Python entry point to decompress Rigaku Oxford Diffractometer TY6 compression
   scitbx::af::flex_int uncompress_rod_TY6(const boost::python::object &data,
                                           const boost::python::object &offsets,
@@ -224,6 +261,11 @@ namespace dxtbx { namespace boost_python {
     def("is_big_endian", is_big_endian);
     def("uncompress", &uncompress, (arg_("packed"), arg_("slow"), arg_("fast")));
     def("compress", &compress);
+    def("radial_average", &radial_average,
+      (arg("data"), arg("beam_center"), arg("sums"), arg("sums_sq"), arg("counts"),
+       arg("pixel_size"), arg("distance"),
+       arg("upper_left"), arg("lower_right")))
+    ;
     def("uncompress_rod_TY6",
         &uncompress_rod_TY6,
         (arg_("data"), arg_("offsets"), arg_("slow"), arg_("fast")));

src/dxtbx/command_line/detector_superpose.py

@@ -8,7 +8,7 @@
 from scitbx.array_family import flex
 from scitbx.math.superpose import least_squares_fit
 from scitbx.matrix import col
-from xfel.command_line.cspad_detector_congruence import iterate_detector_at_level
+from serialtbx.detector import iterate_detector_at_level
 
 import dials.util
 from dials.util.options import OptionParser

src/dxtbx/command_line/image2pickle.py

@@ -16,8 +16,9 @@
 import libtbx.option_parser
 from libtbx import easy_pickle
 from libtbx.utils import Usage
+import serialtbx.util
 from scitbx.array_family import flex
-from xfel.cxi.cspad_ana.cspad_tbx import dpack, evt_timestamp
+import serialtbx.detector.cspad
 
 import dxtbx.util
 
@@ -265,7 +266,7 @@ def run(args=None):
             timestamp = None
         else:
             msec, sec = math.modf(scan.get_epochs()[0])
-            timestamp = evt_timestamp((sec, msec))
+            timestamp = serialtbx.util.timestamp((sec, msec))
 
         if is_multi_image:
             for i in range(img.get_num_images()):
@@ -329,7 +330,7 @@ def save_image(
             print("Skipping %s, file exists" % imgpath)
             return
 
-    data = dpack(
+    data = serialtbx.detector.cspad.dpack(
         data=raw_data,
         distance=distance,
         pixel_size=pixel_size,

src/dxtbx/command_line/radial_average.py

@@ -16,7 +16,7 @@
 from libtbx.utils import Sorry, Usage
 from scitbx.array_family import flex
 from scitbx.matrix import col
-from xfel import radial_average
+from dxtbx.ext import radial_average
 
 import dxtbx
 import dxtbx.util

src/dxtbx/format/FormatCBFMultiTile.py

@@ -9,6 +9,7 @@
 import pycbf
 
 from scitbx.array_family import flex
+from scitbx import matrix
 
 from dxtbx.format.FormatCBF import FormatCBF
 from dxtbx.format.FormatCBFFull import FormatCBFFull
@@ -17,6 +18,13 @@
 from dxtbx.model.detector_helpers import find_undefined_value, find_underload_value
 
 
+def angle_and_axis(basis):
+    """Normalize a quaternion and return the angle and axis
+    @param params metrology object"""
+    q = matrix.col(basis.orientation).normalize()
+    return q.unit_quaternion_as_axis_and_angle(deg=True)
+
+
 class cbf_wrapper(pycbf.cbf_handle_struct):
     """Wrapper class that provides convenience functions for working with cbflib"""
 
@@ -56,6 +64,36 @@ def has_sections(self):
                 return False
             raise e
 
+    def add_frame_shift(self, basis, axis_settings):
+        """Add an axis representing a frame shift (a rotation axis with an offset)"""
+        angle, axis = angle_and_axis(basis)
+
+        if angle == 0:
+            axis = (0, 0, 1)
+
+        if basis.include_translation:
+            translation = basis.translation
+        else:
+            translation = (0, 0, 0)
+
+        self.add_row(
+            [
+                basis.axis_name,
+                "rotation",
+                "detector",
+                basis.depends_on,
+                str(axis[0]),
+                str(axis[1]),
+                str(axis[2]),
+                str(translation[0]),
+                str(translation[1]),
+                str(translation[2]),
+                basis.equipment_component,
+            ]
+        )
+
+        axis_settings.append([basis.axis_name, "FRAME1", str(angle), "0"])
+
 
 class FormatCBFMultiTile(FormatCBFFull):
     """An image reading class multi-tile CBF files"""

src/dxtbx/format/FormatNexusJungfrauExt.py

@@ -7,10 +7,7 @@
 
 from scitbx.array_family import flex
 
-try:
-    from xfel.util.jungfrau import pad_stacked_format
-except ImportError:
-    pass
+from serialtbx.detector.jungfrau import pad_stacked_format
 
 
 from dxtbx.format.FormatNexus import FormatNexus

src/dxtbx/format/FormatPYmultitile.py

@@ -8,10 +8,7 @@
 from iotbx.detectors.npy import image_dict_to_unicode
 from scitbx.matrix import col
 
-try:
-    from xfel.cftbx.detector.cspad_detector import CSPadDetector
-except ImportError:
-    CSPadDetector = None
+from serialtbx.detector.legacy_metrology.cspad_detector import CSPadDetector
 
 from dxtbx.format.FormatPY import FormatPY
 from dxtbx.model import Detector
@@ -20,8 +17,6 @@
 class FormatPYmultitile(FormatPY):
     @staticmethod
     def understand(image_file):
-        if not CSPadDetector:
-            return False
         try:
             with FormatPYmultitile.open_file(image_file, "rb") as fh:
                 data = pickle.load(fh, encoding="bytes")
@@ -60,13 +55,14 @@ def _detector(self):
         ASIC:s are considered, since DXTBX metrology is not concerned
         with hierarchies.
 
-        Merged from xfel.cftbx.detector.cspad_detector.readHeader() and
-        xfel.cftbx.detector.metrology.metrology_as_dxtbx_vectors().
+        Merged from serialtbx.detector.legacy_metrology.cspad_detector.readHeader()
+        and serialtbx.detector.legacy_metrology.metrology.metrology_as_dxtbx_vectors().
         """
 
-        # XXX Introduces dependency on cctbx.xfel!  Should probably be
-        # merged into the code here!
-        from xfel.cftbx.detector.metrology import _transform, get_projection_matrix
+        from serialtbx.detector.legacy_metrology.metrology import (
+            _transform,
+            get_projection_matrix,
+        )
 
         # Apply the detector distance to the translation of the root
         # detector object.

src/dxtbx/format/FormatXTC.py

@@ -15,14 +15,21 @@
 from dxtbx.format.FormatStill import FormatStill
 from dxtbx.model import Spectrum
 from dxtbx.util.rotate_and_average import rotate_and_average
+import serialtbx.util
+import serialtbx.detector.xtc
 
 try:
     import psana
-
-    from xfel.cxi.cspad_ana import cspad_tbx
 except ImportError:
     psana = None
-    cspad_tbx = None
+except TypeError:
+    # Check if SIT_* environment variables are set
+    import os
+
+    if os.environ.get("SIT_ROOT"):
+        # Variables are present, so must have been another error
+        raise
+    psana = None
 
 locator_str = """
   experiment = None
@@ -152,7 +159,7 @@ def understand(image_file):
         If PSANA fails to read it, then input may not be an xtc/smd file. If success, then OK.
         If detector_address is not provided, a command line promp will try to get the address
         from the user"""
-        if not psana or not cspad_tbx:
+        if not psana:
             return False
         try:
             params = FormatXTC.params_from_phil(locator_scope, image_file)
@@ -375,7 +382,7 @@ def get_psana_timestamp(self, index):
         sec = time[0]
         nsec = time[1]
 
-        return cspad_tbx.evt_timestamp((sec, nsec / 1e6))
+        return serialtbx.util.timestamp((sec, nsec / 1e6))
 
     def get_num_images(self):
         return self.n_images
@@ -397,7 +404,7 @@ def _beam(self, index=None):
             if spectrum:
                 wavelength = spectrum.get_weighted_wavelength()
             else:
-                wavelength = cspad_tbx.evt_wavelength(
+                wavelength = serialtbx.detector.xtc.evt_wavelength(
                     evt, delta_k=self.params.wavelength_delta_k
                 )
             if wavelength is None:

src/dxtbx/format/FormatXTCCspad.py

@@ -8,19 +8,13 @@
 from libtbx.phil import parse
 from scitbx.array_family import flex
 from scitbx.matrix import col
-from xfel.cftbx.detector.cspad_cbf_tbx import read_slac_metrology
-from xfel.cxi.cspad_ana.cspad_tbx import env_distance
+from serialtbx.detector import cspad
+from serialtbx.detector.xtc import env_distance
+import serialtbx.detector.cspad
 
 from dxtbx.format.FormatXTC import FormatXTC, locator_str
 from dxtbx.model import Detector, ParallaxCorrectedPxMmStrategy
 
-try:
-    from xfel.cftbx.detector import cspad_cbf_tbx
-    from xfel.cxi.cspad_ana import cspad_tbx
-except ImportError:
-    # xfel not configured
-    pass
-
 cspad_locator_str = """
   cspad {
     detz_offset = None
@@ -92,7 +86,7 @@ def get_raw_data(self, index=None):
         run_number = event.run()
         run = self._psana_runs[run_number]
         det = self._get_psana_detector(run)
-        data = cspad_cbf_tbx.get_psana_corrected_data(
+        data = cspad.get_psana_corrected_data(
             det,
             event,
             use_default=self.params.cspad.use_psana_calib,
@@ -131,7 +125,7 @@ def _detector(self, index=None):
         run = self.get_run_from_index(index)
         det = self._get_psana_detector(run)
         geom = det.pyda.geoaccess(run.run())
-        cob = read_slac_metrology(geometry=geom, include_asic_offset=True)
+        cob = cspad.read_slac_metrology(geometry=geom, include_asic_offset=True)
         distance = env_distance(
             self.params.detector_address[0], run.env(), self.params.cspad.detz_offset
         )
@@ -195,14 +189,12 @@ def _detector(self, index=None):
                         - origin
                     )
                     p.set_local_frame(fast.elems, slow.elems, origin.elems)
-                    p.set_pixel_size(
-                        (cspad_cbf_tbx.pixel_size, cspad_cbf_tbx.pixel_size)
-                    )
-                    p.set_image_size(cspad_cbf_tbx.asic_dimension)
+                    p.set_pixel_size((cspad.pixel_size, cspad.pixel_size))
+                    p.set_image_size(cspad.asic_dimension)
                     p.set_trusted_range(
                         (
-                            cspad_tbx.cspad_min_trusted_value,
-                            cspad_tbx.cspad_saturated_value,
+                            serialtbx.detector.cspad.cspad_min_trusted_value,
+                            serialtbx.detector.cspad.cspad_saturated_value,
                         )
                     )
                     p.set_name(val)

src/dxtbx/format/FormatXTCEpix.py

@@ -63,7 +63,7 @@ def get_detector(self, index=None):
     def _detector(self, index=None):
         from PSCalib.SegGeometryStore import sgs
 
-        from xfel.cftbx.detector.cspad_cbf_tbx import basis_from_geo
+        from serialtbx.detector.xtc import basis_from_geo
 
         run = self.get_run_from_index(index)
         if run.run() in self._cached_detector:

src/dxtbx/format/FormatXTCJungfrau.py

@@ -4,7 +4,17 @@
 import sys
 
 import numpy as np
-import psana
+
+try:
+    import psana
+except ImportError:
+    psana = None
+except TypeError:
+    # Check if SIT_* environment variables are set
+    if os.environ.get("SIT_ROOT"):
+        # Variables are present, so must have been another error
+        raise
+    psana = None
 
 from cctbx import factor_kev_angstrom
 from cctbx.eltbx import attenuation_coefficient
@@ -54,7 +64,7 @@ def understand(image_file):
         return any(["jungfrau" in src.lower() for src in params.detector_address])
 
     def get_raw_data(self, index=None):
-        from xfel.util import jungfrau
+        from serialtbx.detector.util import jungfrau
 
         if index is None:
             index = 0
@@ -94,7 +104,7 @@ def get_detector(self, index=None):
     def _detector(self, index=None):
         from PSCalib.SegGeometryStore import sgs
 
-        from xfel.cftbx.detector.cspad_cbf_tbx import basis_from_geo
+        from serialtbx.detector.xtc import basis_from_geo
 
         run = self.get_run_from_index(index)
         if run.run() in self._cached_detector: