Make use of ringextraction

doc/source/usage/tutorial/Recalib/Recalib_notebook.ipynb

@@ -46,7 +46,7 @@
     "import fabio\n",
     "from pyFAI.test.utilstest import UtilsTest\n",
     "from pyFAI.calibrant import CALIBRANT_FACTORY\n",
-    "from pyFAI.goniometer import SingleGeometry\n",
+    "from pyFAI.single_geometry import SingleGeometry\n",
     "print(f\"Using pyFAI version: {pyFAI.version}\")\n",
     "start_time = time.perf_counter()"
    ]

src/pyFAI/goniometer.py

@@ -44,12 +44,7 @@
 import numpy
 from collections import OrderedDict, namedtuple
 from scipy.optimize import minimize
-from silx.image import marchingsquares
-from .massif import Massif
-from .control_points import ControlPoints
 from .detectors import detector_factory, Detector
-from .geometry import Geometry
-from .geometryRefinement import GeometryRefinement
 from .azimuthalIntegrator import AzimuthalIntegrator
 from .utils import StringTypes
 from .multi_geometry import MultiGeometry
@@ -581,168 +576,6 @@ def sload(cls, filename):
         return gonio
 
 
-class SingleGeometry(object):
-    """This class represents a single geometry of a detector position on a
-    goniometer arm
-    """
-
-    def __init__(self, label, image=None, metadata=None, pos_function=None,
-                 control_points=None, calibrant=None, detector=None, geometry=None):
-        """Constructor of the SingleGeometry class, used for calibrating a
-        multi-geometry setup with a moving detector.
-
-        :param label: name of the geometry, a string or anything unmutable
-        :param image: image with Debye-Scherrer rings as 2d numpy array
-        :param metadata: anything which contains the goniometer position
-        :param pos_function: a function which takes the metadata as input
-                                 and returns the goniometer arm position
-        :param control_points: a pyFAI.control_points.ControlPoints instance
-            (optional parameter)
-        :param calibrant: a pyFAI.calibrant.Calibrant instance.
-                        Contains the wavelength to be used (optional parameter)
-        :param detector: a pyFAI.detectors.Detector instance or something like
-                        that Contains the mask to be used (optional parameter)
-        :param geometry: an azimuthal integrator or a ponifile
-                        (or a dict with the geometry) (optional parameter)
-        """
-        dict_geo = {}
-        self.label = label
-        self.image = image
-        self.metadata = metadata  # may be anything
-        self.calibrant = calibrant
-        if control_points is None or isinstance(control_points, ControlPoints):
-            self.control_points = control_points
-        else:
-            # Probaly a NPT file
-            self.control_points = ControlPoints(control_points, calibrant=calibrant)
-
-        if detector is not None:
-            self.detector = detector_factory(detector)
-        else:
-            self.detector = None
-        if isinstance(geometry, Geometry):
-            dict_geo = geometry.getPyFAI()
-        elif isinstance(geometry, StringTypes) and os.path.exists(geometry):
-            dict_geo = Geometry.sload(geometry).getPyFAI()
-        elif isinstance(geometry, dict):
-            dict_geo = geometry
-
-        if self.detector is not None:
-            dict_geo["detector"] = self.detector
-        if self.control_points is not None:
-            dict_geo["data"] = self.control_points.getList()
-        if self.calibrant is not None:
-            dict_geo["calibrant"] = self.calibrant
-            if self.calibrant.wavelength:
-                dict_geo["wavelength"] = self.calibrant.wavelength
-        if "max_shape" in dict_geo:
-            # not used in constructor
-            dict_geo.pop("max_shape")
-        self.geometry_refinement = GeometryRefinement(**dict_geo)
-        if self.detector is None:
-            self.detector = self.geometry_refinement.detector
-        self.pos_function = pos_function
-        self.massif = None
-
-    def get_position(self):
-        """This method  is in charge of calculating the motor position from metadata/label/..."""
-        return self.pos_function(self.metadata)
-
-    def extract_cp(self, max_rings=None, pts_per_deg=1.0, Imin=0):
-        """Performs an automatic keypoint extraction and update the geometry refinement part
-
-        :param max_ring: extract at most N rings from the image
-        :param pts_per_deg: number of control points per azimuthal degree (increase for better precision)
-        """
-        if self.massif is None:
-            if self.detector:
-                mask = self.detector.dynamic_mask(self.image)
-            else:
-                mask = None
-            self.massif = Massif(self.image, mask)
-
-        tth = numpy.array([i for i in self.calibrant.get_2th() if i is not None])
-        tth = numpy.unique(tth)
-        tth_min = numpy.zeros_like(tth)
-        tth_max = numpy.zeros_like(tth)
-        delta = (tth[1:] - tth[:-1]) / 4.0
-        tth_max[:-1] = delta
-        tth_max[-1] = delta[-1]
-        tth_min[1:] = -delta
-        tth_min[0] = -delta[0]
-        tth_max += tth
-        tth_min += tth
-        shape = self.image.shape
-        ttha = self.geometry_refinement.twoThetaArray(shape)
-        chia = self.geometry_refinement.chiArray(shape)
-        rings = 0
-        cp = ControlPoints(calibrant=self.calibrant)
-        if max_rings is None:
-            max_rings = tth.size
-
-        qmask, count = build_qmask(ttha, tth_min, tth_max, self.geometry_refinement.detector.mask)
-        mask2 = numpy.empty(qmask.shape, dtype=bool)
-        ms = marchingsquares.MarchingSquaresMergeImpl(ttha,
-                                                      mask=self.geometry_refinement.detector.mask,
-                                                      use_minmax_cache=True)
-        for i in range(tth.size):
-            if rings >= max_rings:
-                break
-            if count[i]:
-                rings += 1
-                mask = qmask == i
-                sub_data = self.image[mask]
-                mean = sub_data.mean(dtype=numpy.float64)
-                std = sub_data.std(dtype=numpy.float64)
-                upper_limit = mean + std
-                numpy.logical_and(self.image > upper_limit, mask, out=mask2)
-                size2 = mask2.sum(dtype=int)
-                if size2 < 1000:
-                    upper_limit = mean
-                    numpy.logical_and(self.image > upper_limit, mask, out=mask2)
-                    size2 = mask2.sum()
-                # length of the arc:
-                points = ms.find_pixels(tth[i])
-                seeds = set((i[0], i[1]) for i in points if mask2[i[0], i[1]])
-                # max number of points: 360 points for a full circle
-                azimuthal = chia[points[:, 0].clip(0, shape[0]), points[:, 1].clip(0, shape[1])]
-                nb_deg_azim = numpy.unique(numpy.rad2deg(azimuthal).round()).size
-                keep = int(nb_deg_azim * pts_per_deg)
-                if keep == 0:
-                    continue
-                dist_min = len(seeds) / 2.0 / keep
-                # why 3.0, why not ?
-
-                logger.info("Extracting datapoint for ring %s (2theta = %.2f deg); " +
-                            "searching for %i pts out of %i with I>%.1f, dmin=%.1f",
-                            i, numpy.degrees(tth[i]), keep, size2, upper_limit, dist_min)
-                res = self.massif.peaks_from_area(mask2, Imin=Imin, keep=keep, dmin=dist_min, seed=seeds, ring=i)
-                cp.append(res, i)
-        self.control_points = cp
-        self.geometry_refinement.data = numpy.asarray(cp.getList(), dtype=numpy.float64)
-        return cp
-
-    def get_ai(self):
-        """Create a new azimuthal integrator to be used.
-
-        :return: Azimuthal Integrator instance
-        """
-        config = self.geometry_refinement.get_config()
-        ai = AzimuthalIntegrator()
-        ai.set_config(config)
-        return ai
-
-    def get_wavelength(self):
-        assert self.calibrant.wavelength == self.geometry_refinement.wavelength
-        return self.geometry_refinement.wavelength
-
-    def set_wavelength(self, value):
-        self.calibrant.setWavelength_change2th(value)
-        self.geometry_refinement.set_wavelength(value)
-
-    wavelength = property(get_wavelength, set_wavelength)
-
-
 class GoniometerRefinement(Goniometer):
     """This class allow the translation of a goniometer geometry into a pyFAI
     geometry using a set of parameter to refine.

src/pyFAI/meson.build

@@ -33,6 +33,7 @@ py.install_sources(
     'multi_geometry.py',
     'parallax.py',
     'ring_extraction.py',
+    'single_geometry.py',
     'spline.py',
     'units.py',
     'worker.py'],

src/pyFAI/ring_extraction.py

@@ -50,30 +50,40 @@
 from silx.image import marchingsquares
 
 from .control_points import ControlPoints
-from .goniometer import SingleGeometry
 from .massif import Massif
+from .calibrant import Calibrant
+from .geometryRefinement import GeometryRefinement
 
 logger = logging.getLogger(__name__)
 
 
 class RingExtraction:
     """Class to perform extraction of control points from a calibration image."""
 
-    def __init__(self, single_geometry: SingleGeometry, massif: Optional[Massif] = None):
+    def __init__(
+        self, image: numpy.ndarray, detector, calibrant: Calibrant,
+        geometry_refinement: GeometryRefinement, massif: Optional[Massif] = None,
+    ):
         """
         Parameters
         ----------
-        single_geometry : SingleGeometry
-            instance of pyFAI class
+        image : numpy.ndarray
+            image with Debye-Scherrer rings as 2d numpy array
+        detector : pyFAI.detectors.Detector
+            a pyFAI.detectors.Detector instance or something like that contains the mask to be used
+        calibrant : Calibrant
+            instance of Calibrant class
+        geometry_refinement : GeometryRefinement
+            instance of GeometryRefinement class
         massif : Optional[Massif]
             instance of Massif class, which defines an area around a peak; it is used to find
             neighboring peaks, by default None
         """
 
-        self.single_geometry = single_geometry
-        self.image = self.single_geometry.image
-        self.detector = self.single_geometry.detector
-        self.calibrant = self.single_geometry.calibrant
+        self.image = image
+        self.detector = detector
+        self.calibrant = calibrant
+        self.geometry_refinement = geometry_refinement
 
         if massif:
             self.massif = massif
@@ -84,11 +94,12 @@ def __init__(self, single_geometry: SingleGeometry, massif: Optional[Massif] = N
                 mask = None
             self.massif = Massif(self.image, mask)
 
-        self.two_theta_array = self.single_geometry.geometry_refinement.twoThetaArray()
+        self.two_theta_array = self.geometry_refinement.twoThetaArray()
         self.two_theta_values = self._get_unique_two_theta_values_in_image()
 
     def extract_control_points(
-        self, max_number_of_rings: Optional[int] = None, points_per_degree: float = 1
+        self, max_number_of_rings: Optional[int] = None,
+        points_per_degree: float = 1, min_intensity: float = 0,
     ) -> ControlPoints:
         """
         Primary method of RingExtraction class. Runs extract_control_points_in_one_ring for all
@@ -101,6 +112,8 @@ def extract_control_points(
         points_per_degree : float, optional
             number of control points per azimuthal degree (increase for better precision), by
             default 1.0
+        min_intensity : float
+            minimum of intensity above the background to keep the point
 
         Returns
         -------
@@ -117,7 +130,10 @@ def extract_control_points(
         with ThreadPoolExecutor() as executor:
             for ring_index in range(min(max_number_of_rings, self.two_theta_values.size)):
                 future = executor.submit(
-                    self.extract_list_of_peaks_in_one_ring, ring_index, points_per_degree
+                    self.extract_list_of_peaks_in_one_ring,
+                    ring_index,
+                    min_intensity,
+                    points_per_degree,
                 )
                 tasks[future] = ring_index
 
@@ -130,7 +146,10 @@ def extract_control_points(
         return control_points
 
     def extract_list_of_peaks_in_one_ring(
-        self, ring_index: int, points_per_degree: float = 1.0
+        self,
+        ring_index: int,
+        min_intensity: float,
+        points_per_degree: float = 1.0,
     ) -> Optional[list[tuple[float, float]]]:
         """
         Using massif.peaks_from_area, get all pixel coordinates inside a mask of pixels around a
@@ -141,6 +160,8 @@ def extract_list_of_peaks_in_one_ring(
         ----------
         ring_index : int
             ring number
+        min_intensity : float
+            minimum of intensity above the background to keep the point
         points_per_degree : float, optional
             number of control points per azimuthal degree (increase for better precision), by
             default 1.0
@@ -188,6 +209,7 @@ def extract_list_of_peaks_in_one_ring(
 
                 return self.massif.peaks_from_area(
                     final_mask,
+                    Imin=min_intensity,
                     keep=num_points_to_keep,
                     dmin=min_distance_between_control_points,
                     seed=seeds,
@@ -342,7 +364,9 @@ def _calculate_num_of_points_to_keep(
             Number of points to keep as control points
         """
         image_shape = self.image.shape
-        azimuthal_angles_array = self.single_geometry.geometry_refinement.chiArray(image_shape)
+        azimuthal_angles_array = self.geometry_refinement.chiArray(
+            image_shape
+        )
         azimuthal_degrees_array_in_ring = azimuthal_angles_array[
             pixels_at_two_theta_level[:, 0].clip(0, image_shape[0]),
             pixels_at_two_theta_level[:, 1].clip(0, image_shape[1]),
@@ -391,7 +415,7 @@ def _get_beam_centre_coords(self) -> numpy.ndarray:
         numpy.ndarray
             beam centre coordinates [y,x], to match pyFAI order.
         """
-        fit_2d = self.single_geometry.get_ai().getFit2D()
+        fit_2d = self.geometry_refinement.getFit2D()
         beam_centre_x = fit_2d["centerX"]
         beam_centre_y = fit_2d["centerY"]
         return numpy.array([beam_centre_y, beam_centre_x])