From 4c07a35028f7fc19e29d806ffbe3a6cbf59e68c1 Mon Sep 17 00:00:00 2001
From: OliDG <35495131+OliDG@users.noreply.github.com>
Date: Fri, 28 Jun 2024 13:28:45 +0200
Subject: [PATCH] Auto masking tool (#46)

* Update center.py

Tool to automatically generate a mask based on the darkest portion of the image. This to skip the manual definition of one and to compensate for possible beam stopper movement (image rotation during acquisition or beam stopper mechanical play).

* formatting, test and documentation

black formatting check, added to the loading by default and updated tutorials section

* improved to avoid hotspot problems, more universal

* Update __init__.py

Black formatting

---------

Co-authored-by: Olivier Donzel-Gargand <olido551@user.uu.se>
---
 docs/tutorials/image.rst | 36 ++++++++++++++++++++++++++++++++++++
 skued/__init__.py        |  1 +
 skued/image/__init__.py  |  2 +-
 skued/image/center.py    | 26 ++++++++++++++++++++++++++
 4 files changed, 64 insertions(+), 1 deletion(-)

diff --git a/docs/tutorials/image.rst b/docs/tutorials/image.rst
index 910146be..387fc73d 100644
--- a/docs/tutorials/image.rst
+++ b/docs/tutorials/image.rst
@@ -27,7 +27,43 @@ Diffraction patterns can come in a variety of exotic file formats. Scikit-ued ha
 * All other file formats supported by `scikit-image`_.
 
 The :func:`diffread` function will transparently distinguish between those formats and dispatch to the right functions. 
+.. _automatic_mask_generation:
 
+Automatic generation of a mask as needed for the autocenter() function
+======================================================================
+
+Automatic mask generation
+-------------------------
+A mask can be used for the :func: `autocenter` function to exclude the beam stopper from the analysis and only keep valuable 
+information. Such a mask can automatically be generated using this function; given a diffraction pattern, the mask will be created to block 
+the darkest regions of the image (default blocks about 10% of the darker values).
+Here is an example: 
+
+.. plot::
+
+	import matplotlib.pyplot as plt
+	from skued import diffread, auto_masking
+
+	image = diffread('data/Cr_1.tif')
+	mask = auto_masking(image, threshold = 0.1)
+
+	# Reduce size of images because of memory usage of ReadTheDocs
+	image = image[::3, ::3]
+	mask = mask[::3, ::3]
+
+	fig , (ax1, ax2) = plt.subplots(1,2, figsize = (9,3))
+	ax1.imshow(image, vmin=0, vmax=150, cmap='inferno')
+	ax2.imshow(mask, vmin=0, vmax=1, cmap='inferno')
+
+	for ax in (ax1, ax2):
+		ax.xaxis.set_visible(False)
+		ax.yaxis.set_visible(False)
+
+	ax1.set_title('Cr_1')
+	ax2.set_title('Mask')
+
+	plt.tight_layout()
+	plt.show()
 .. _alignment:
 
 Automatic center-finding
diff --git a/skued/__init__.py b/skued/__init__.py
index ce5e55af..fdcc97b5 100644
--- a/skued/__init__.py
+++ b/skued/__init__.py
@@ -44,6 +44,7 @@
 from .image import (
     align,
     autocenter,
+    auto_masking,
     azimuthal_average,
     brillouin_zones,
     combine_masks,
diff --git a/skued/image/__init__.py b/skued/image/__init__.py
index 0c57b242..3443bf15 100644
--- a/skued/image/__init__.py
+++ b/skued/image/__init__.py
@@ -4,7 +4,7 @@
 from .alignment import align, ialign, itrack_peak
 from .brillouin import brillouin_zones
 from .calibration import detector_scattvectors, powder_calq
-from .center import autocenter
+from .center import autocenter, auto_masking
 from .indexing import bragg_peaks, bragg_peaks_persistence
 from .metrics import (
     combine_masks,
diff --git a/skued/image/center.py b/skued/image/center.py
index 106a2f5b..0db83606 100644
--- a/skued/image/center.py
+++ b/skued/image/center.py
@@ -124,3 +124,29 @@ def _center_of_intensity(im, mask=None):
     r_ = np.average(rr, weights=weights)
     c_ = np.average(cc, weights=weights)
     return int(r_), int(c_)
+
+
+def auto_masking(im, threshold=0.1):
+    """
+    Generate a mask based on the darkest fraction of an image
+
+    Parameters
+    ----------
+    im : floats, ndarrays of shape (N,M)
+        image used to generate a mask
+    threshold: float, optional
+        fraction of the lowest values to be masked, default = 15%
+
+    Yields
+    ------
+    mask : boolean, ndarrays of shape (N,M)
+        Mask that evaluates to True on valid pixels.
+
+    """
+    # Find the median of the highest intensity value of the image to avoid hot spots
+    max_median = np.median([max(x) for x in np.real(im)])
+    # Set the threshold value
+    lower_limit = threshold * max_median
+    # generate a mask
+    mask = im >= lower_limit
+    return mask