a bug fix in the stitching

islatu/image.py

@@ -31,7 +31,6 @@ class Image:
         data (:py:class:`pandas.DataFrame`, optional): Experimental data about the measurement. Defaults to :py:attr:`None`.
         metadata (:py:attr:`dict`, optional): Metadata regarding the measurement. Defaults to :py:attr:`None`.
         transpose (:py:attr:`bool`, optional): Should the data be rotated by 90 degrees? Defaults to :py:attr:`False`.
-        hot_pixel_threshold (:py:attr:`int`, optional): The number of counts above which a pixel should be assessed to determine if it is hot. Defaults to :py:attr:`200000`.
         pixel_maximum (:py:attr:`int`, optional): The number of counts above which a pixel should be assessed to determine if it is hot. Defaults to :py:attr:`500000`.
         pixel_minimum (:py:attr:`int`, optional): The number of counts above which a pixel should be assessed to determine if it is hot. Defaults to :py:attr:`0`.
     """
@@ -42,12 +41,10 @@ def __init__(
         data=None,
         metadata=None,
         transpose=False,
-        hot_pixel_threshold=200000,
-        pixel_maximum=500000,
         pixel_minimum=0
     ):
         """
-        Initialisation of the :py:class:`islatu.image.Image` class, includes running hot pixel check and assigning uncertainties.
+        Initialisation of the :py:class:`islatu.image.Image` class, includes assigning uncertainties.
         """
         self.file_path = file_path
         self.data = data
@@ -58,8 +55,6 @@ def __init__(
         if transpose:
             array = array.T
         # Remove hot pixels
-        array = _find_hot_pixels(array, threshold=hot_pixel_threshold)
-        array[np.where(array > pixel_maximum)] = 0
         array[np.where(array < pixel_minimum)] = 0
         array_error = np.sqrt(array)
         array_error[np.where(array == 0)] = 1
@@ -181,40 +176,3 @@ def sum(self, axis=None):
         """
         return self.array.sum(axis)
 
-
-def _find_hot_pixels(array, threshold=200000):
-    """
-    Find some dead pixels and assign them with some local average value.
-
-    Args:
-        array (:py:attr:`array_like`): NumPy array describing the image.
-        threshold (:py:attr:`int`, optional): The number of counts above which a pixel should be assessed to determine if it is hot. Defaults to :py:attr:`200000`.
-
-    Returns:
-        :py:attr:`array_like`: NumPy array where hot pixels have been removed.
-    """
-    sorted_array = np.sort(array.flatten())[::-1]
-    for i in sorted_array:
-        if i >= threshold:
-            pos_a, pos_b = np.where(array == i)
-            lower_a = pos_a[0] - 1
-            upper_a = pos_a[0] + 2
-            lower_b = pos_b[0] - 1
-            upper_b = pos_b[0] + 2
-            if pos_a[0] == 0:
-                lower_a = 0
-            elif pos_a[0] == array.shape[0] - 1:
-                upper_a = array.shape[0]
-            if pos_b[0] == 0:
-                lower_b = 0
-            elif pos_b[0] == array.shape[1] - 1:
-                upper_b = array.shape[1]
-            local = np.copy(array[lower_a:upper_a, lower_b:upper_b])
-            local[np.where(local == i)] = 0
-            local_average = np.mean(local)
-            local_std = np.std(local)
-            if i > local_average + 2 * local_std:
-                array[pos_a[0], pos_b[0]] = local_average
-        else:
-            break
-    return array

islatu/refl_data.py

@@ -32,11 +32,15 @@ class Profile:
     """
 
     def __init__(
-        self, file_paths, parser, q_axis_name="qdcd", theta_axis_name="dcdtheta"
+        self, file_paths, parser, q_axis_name="qdcd", theta_axis_name="dcdtheta", pixel_max=10000000, transpose=False
     ):
         self.scans = []
-        for f in file_paths:
-            self.scans.append(Scan(f, parser, q_axis_name, theta_axis_name))
+        try:
+            is_list = len(pixel_max)
+        except TypeError:
+            pixel_max = [pixel_max] * len(file_paths)
+        for i, f in enumerate(file_paths):
+            self.scans.append(Scan(f, parser, q_axis_name, theta_axis_name, pixel_max=pixel_max[i], transpose=transpose))
         self.q_vectors = None
         self.reflected_intensity = None
 
@@ -211,6 +215,8 @@ def __init__(
         q_axis_name="qdcd",
         theta_axis_name="dcdtheta",
         energy=None,
+        pixel_max=1000000,
+        transpose=False,
     ):
         self.file_path = file_path
         self.metadata, self.data = parser(self.file_path)
@@ -237,6 +243,8 @@ def __init__(
         )
         self.R = unp.uarray(np.zeros(self.q.size), np.zeros(self.q.size),)
         self.n_pixels = np.zeros(self.q.size)
+        self.pixel_max = pixel_max
+        self.transpose = transpose
 
     def __str__(self):
         """
@@ -310,21 +318,35 @@ def crop_and_bkg_sub(
             progress (:py:attr:`bool`, optional): Show a progress bar. Requires the :py:mod:`tqdm` package. Defaults to :py:attr:`True`.
         """
         iterator = _get_iterator(unp.nominal_values(self.q), progress)
+        to_remove = []
         for i in iterator:
-            im = image.Image(self.data["file"][i], self.data, self.metadata)
-            if crop_kwargs is None:
-                im.crop(crop_function)
-            else:
-                im.crop(crop_function, **crop_kwargs)
-            if bkg_sub_kwargs is None:
-                im.background_subtraction(bkg_sub_function)
-            else:
-                im.background_subtraction(bkg_sub_function, **bkg_sub_kwargs)
-            self.R[i] = ufloat(im.sum().n, im.sum().s)
-            if im.n_pixel is None:
-                self.n_pixels[i] = None
-            else:
-                self.n_pixels[i] = im.n_pixel.n
+            try:
+                if self.data['roi1_maxval'][i] > self.pixel_max:
+                    to_remove.append(i)
+                    continue
+                im = image.Image(self.data["file"][i], self.data, self.metadata, transpose=self.transpose)
+                if crop_kwargs is None:
+                    im.crop(crop_function)
+                else:
+                    im.crop(crop_function, **crop_kwargs)
+                if bkg_sub_kwargs is None:
+                    im.background_subtraction(bkg_sub_function)
+                else:
+                    im.background_subtraction(bkg_sub_function, **bkg_sub_kwargs)
+                self.R[i] = ufloat(im.sum().n, im.sum().s)
+                if im.n_pixel is None:
+                    self.n_pixels[i] = None
+                else:
+                    self.n_pixels[i] = im.n_pixel.n
+            except (ValueError, RuntimeError) as e:
+                print(e)
+                print('Current file: ', self.data["file"][i])
+                to_remove.append(i)
+                continue
+        self.R = np.delete(self.R, to_remove)
+        self.theta = np.delete(self.theta, to_remove)
+        self.n_pixels = np.delete(self.n_pixels, to_remove)
+        self.q = np.delete(self.q, to_remove)
 
     def footprint_correction(self, beam_width, sample_size):
         """

islatu/stitching.py

@@ -5,9 +5,10 @@
 # Copyright (c) Andrew R. McCluskey
 # Distributed under the terms of the MIT License
 # author: Andrew R. McCluskey
+import warnings
 import numpy as np
 from uncertainties import unumpy as unp
-
+from scipy.stats import linregress
 
 def correct_attentuation(scan_list):
     """
@@ -22,15 +23,22 @@ def correct_attentuation(scan_list):
     for i in range(len(scan_list) - 1):
         overlap_start = scan_list[i + 1].q[0].n
         overlap_end = scan_list[i].q[-1].n
-
-        overlap_start_index = np.argmin(np.abs(scan_list[i].q - overlap_start))
-        overlap_end_index = np.argmin(np.abs(scan_list[i + 1].q - overlap_end))
-
-        target_r = scan_list[i].R[overlap_start_index:]
-        vary_r = scan_list[i + 1].R[: overlap_end_index + 1]
+
+        if overlap_start > overlap_end:
+            warnings.warn('Using extrapolation to correct attenuation between scans {} and {}. Please double check these results.'.format(
+                scan_list[i].file_path, scan_list[i+1].file_path))
+            overlap_start_index = -2
+            overlap_end_index = 1
+        else:
+            overlap_start_index = np.argmin(np.abs(scan_list[i].q - overlap_start))
+            overlap_end_index = np.argmin(np.abs(scan_list[i + 1].q - overlap_end))
+
+        res = linregress(unp.nominal_values(scan_list[i].q[overlap_start_index:]), np.log(unp.nominal_values(scan_list[i].R[overlap_start_index:])))
+        target_r = unp.exp(scan_list[i+1].q[: overlap_end_index + 1] * res.slope + res.intercept)
+        vary_r = scan_list[i+1].R[: overlap_end_index + 1]
 
         ratio = target_r.mean() / vary_r.mean()
-
+        
         scan_list[i + 1].R *= ratio
     return scan_list
 

islatu/tests/test_image.py

@@ -97,57 +97,6 @@ def test_init(self):
         assert_equal((10, 10), test_image.shape)
         assert_almost_equal(expected_image, test_image.n)
 
-    def test_hot_pixel(self):
-        """
-        Test hot pixel in image
-        """
-        b = io.StringIO(EXAMPLE_HOT_PIXEL)
-        buf = io.BytesIO()
-        im = PILIm.fromarray(np.loadtxt(b).astype(np.uint32))
-        im.save(buf, format="png")
-        buf.seek(0)
-        test_image = Image(buf)
-        test_image = image._find_hot_pixels(test_image.n, threshold=6e03)
-        data = io.StringIO(EXAMPLE_HOT_PIXEL)
-        expected_image = np.loadtxt(data)
-        expected_image[2, 2] = 2.333333333333
-        assert_equal((10, 10), test_image.shape)
-        assert_almost_equal(expected_image, test_image)
-
-    def test_hot_pixel_corner_a(self):
-        """
-        Test a hot pixel in the top left corner
-        """
-        b = io.StringIO(EXAMPLES_HOT_PIXEL_CORNERA)
-        buf = io.BytesIO()
-        im = PILIm.fromarray(np.loadtxt(b).astype(np.uint32))
-        im.save(buf, format="png")
-        buf.seek(0)
-        test_image = Image(buf)
-        test_image = image._find_hot_pixels(test_image.n, threshold=6e03)
-        data = io.StringIO(EXAMPLES_HOT_PIXEL_CORNERA)
-        expected_image = np.loadtxt(data)
-        expected_image[0, 0] = 0
-        assert_equal((10, 10), test_image.shape)
-        assert_almost_equal(expected_image, test_image)
-
-    def test_hot_pixel_corner_b(self):
-        """
-        Test a hot pixel in the bottom right corner
-        """
-        b = io.StringIO(EXAMPLES_HOT_PIXEL_CORNERB)
-        buf = io.BytesIO()
-        im = PILIm.fromarray(np.loadtxt(b).astype(np.uint32))
-        im.save(buf, format="png")
-        buf.seek(0)
-        test_image = Image(buf)
-        test_image = image._find_hot_pixels(test_image.n, threshold=6e03)
-        data = io.StringIO(EXAMPLES_HOT_PIXEL_CORNERB)
-        expected_image = np.loadtxt(data)
-        expected_image[9, 9] = 3.25
-        assert_equal((10, 10), test_image.shape)
-        assert_almost_equal(expected_image, test_image)
-
     def test_init_with_transpose(self):
         """
         Test for transposing with reading

islatu/tests/test_stitching.py

@@ -52,7 +52,7 @@ def test_correct_attentuation(self):
         exp_y1 = np.ones(10) * 10
         exp_dy1 = np.ones(10)
         exp_y2 = np.ones(10) * 10
-        exp_dy2 = np.append(np.ones(3), np.ones(7) * 1.2909944)
+        exp_dy2 = np.append(np.ones(3) * 0.8164966, np.ones(7) * 1.1547005)
 
         scan_list = stitching.correct_attentuation(scan_list)
 

setup.py

@@ -15,7 +15,7 @@
 # versioning
 MAJOR = 0
 MINOR = 0
-MICRO = 1
+MICRO = 26 
 ISRELEASED = False
 VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO)