add buggy LM exitation curve correction for test

.gitignore

+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/latest/usage/project/#working-with-version-control
+.pdm.toml
+.pdm-python
+.pdm-build/
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
\ No newline at end of file

ditb/commands/darkwhite.py

@@ -20,7 +20,8 @@ import click
 @click.option('--force', is_flag=True, help="Force to normalize even if the directory is already in the destination folder")
 @click.option('--lower', type=click.FloatRange(0.01,100, clamp=True), help="Set the lower percentile of pixel value set them to zero after normalisation. (The default value is 1)", default=1)
 @click.option('--smooth', type=click.BOOL, default=True, help="Smooth image of white, darks and dark of whites")
-def darkwhite(images, dark, white, darkofwhite, white_z, force, lower, smooth):
+@click.option('--correct_lm_exi', type=click.BOOL, default=False, help="Use the excitation profile of the Lame Matthieu to correct excitation scann of the white")
+def darkwhite(images, dark, white, darkofwhite, white_z, force, lower, smooth, correct_lm_exi):
     """
     🔬 Normalise images with dark (images of the camera noise), 
     white (image stack of the DISCO beam) and "dark of white" images.
@@ -82,7 +83,10 @@ def darkwhite(images, dark, white, darkofwhite, white_z, force, lower, smooth):
     for filein in images:
         fname = Path(filein).name
         if fname not in darkwhite_names or 'normalized' not in fname.lower() or 'dark' not in fname.lower or 'white' not in fname.lower:
-            output_file = Path(filein).parent / Path(filein).name.replace('.zarr', '_normalized.zarr') 
+            out_suffix = '_normalized'
+            if correct_lm_exi:
+                out_suffix = '_normalized_LME'
+            output_file = Path(filein).parent / Path(filein).name.replace('.zarr', f'{out_suffix}.zarr') 
             print(f'\n▷ normalize {filein} to {output_file}')
         
             if force:
@@ -91,7 +95,7 @@ def darkwhite(images, dark, white, darkofwhite, white_z, force, lower, smooth):
                 if output_file.is_dir():
                     shutil.rmtree(str(output_file))
                 try:
-                    correct_dark_white(filein, dark, white, darkofwhite, white_z, lower_percentile=lower)
+                    correct_dark_white(filein, dark, white, darkofwhite, white_z, lower_percentile=lower, correct_LM_excitation=correct_lm_exi)
                 except Exception as e:
                     print(f'Error for {filein}')
                     print(e)
@@ -99,7 +103,7 @@ def darkwhite(images, dark, white, darkofwhite, white_z, force, lower, smooth):
             else:
                 if not output_file.is_dir():
                     try:
-                        correct_dark_white(filein, dark, white, darkofwhite, white_z, lower_percentile=lower)
+                        correct_dark_white(filein, dark, white, darkofwhite, white_z, lower_percentile=lower, correct_LM_excitation=correct_lm_exi)
                     except Exception as e:
                         print(f'Error for {filein}')
                         print(e)

ditb/exposure/darkwhite_correction.py

@@ -16,7 +16,7 @@ import matplotlib.pyplot as plt
 def correct_dark_white(images: str, dark: str = '', white: str = '' ,
                        darkofwhite: str = '', white_z: [int, str] = 'center',
                        lower_percentile: float = 1, smoothsize: int = 2,
-                       remove_cosmics = False):
+                       remove_cosmics = False, correct_LM_excitation=False):
     """
     Do the dark white correction of images.
     
@@ -31,6 +31,16 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
         the path to the dark images of the stack
     - white: str,
         the path to the white
+    - darkofwhite: str,
+        the path to the dark of white images
+    - white_z: int or 'center':
+        The plane in the z-stack of the white to use in the correction. 
+        This option is only used for regular stack correction
+    - remove_cosmics: bool:
+        Use and opening/closing filter to remove cosmics on the image
+    - correct_LM_excitation: bool,
+        Use the excitation spectra of the Lame Matthieu to correct the excitation 
+        scanning profile.
     """
     
     whitein = Path(white)
@@ -44,7 +54,11 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
             raise ValueError(f'the file {f} is not a zarr file')
     
     # Prepare output name
-    out_path = imagesin.parent / imagesin.name.replace('.zarr', '_normalized.zarr')
+    out_suffix = '_normalized'
+    if correct_LM_excitation:
+        out_suffix = '_normalized_LME'
+
+    out_path = imagesin.parent / imagesin.name.replace('.zarr', f'{out_suffix}.zarr')
     
     white = mmzarr_read(str(whitein))
     meta_white = read_mmzarr_metadata(str(whitein))
@@ -66,9 +80,20 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
             dow_stack = dow_stack.map_blocks(remove_cosmic, 2, dtype=dow_stack.dtype)
         
         diffwhite = white_stack - dow_stack
+        # If we correct by the white excitation, we need to 
+        # take the maximum of the white for each excitation
+        # wavelength, to only correct for spatial excitation 
+        # inhomogenity
+        if correct_LM_excitation:
+            print('process White for LME')
+            max_white = diffwhite.max(axis=(-2,-1))
+            min_ratio = white_stack.min() / max_white
+            normed_white = diffwhite / max_white[:, None, None]
+            # Cut at the maximum of the minimum ratio
+            normed_white[normed_white < min_ratio.max()] = min_ratio.max()
+        else:
             max_white = diffwhite.max()
             min_ratio = white_stack.min() / max_white
-
             normed_white = diffwhite / max_white
             normed_white[normed_white < min_ratio] = min_ratio
 
@@ -161,6 +186,7 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
         # Save an image of the normalised white profile, and the percentile
         # value as a function of the excitation wavelength
         w = [float(get_mono_position(meta_white, i, excitation_as_z=True)) for i in range(n_white.shape[2])]
+        """
         plt.figure(figsize=(15,4))
         plt.subplot(131)
         plt.plot(w, n_white.max(axis=(-2,-1)).squeeze())
@@ -183,12 +209,45 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
         plt.suptitle(str(imagesin.name))
         plt.tight_layout()
         plt.savefig(str(imagesin.parent / imagesin.name.replace('.zarr', '_normalized.pdf')))
+        """
+
+        # Compute the transfert function
+        if correct_LM_excitation:
+            # The excitation of the white is not flat
+            print('Use the transfer function from the excitation spectra of the LM')
+            xlm, ylm = load_excitation_LM()
+
+            # Interpolate the excitation spectra of the white
+            # 07/05/24 seems thats their is an offset of 7nm in wavelength on 
+            # TELEMOS mono compared to Fluoromax
+            ylmi = np.interp(w, xlm+7, ylm)
+
+            # Normalise this profile by it's maximum
+            ylmi_norm = ylmi / ylmi.max()
+
+            # Compute the normed excitation profile of the white
+            max_white = diffwhite.max()
+            min_ratio = white_stack.min() / max_white
+            normed_white = diffwhite / max_white
+            normed_white[normed_white < min_ratio] = min_ratio
+            normed_white = normed_white.compute()
+
+            # Correct from spatial difference in the white
+            # tf = 1/n_white[None,...]
+            tf = 1
+
+            # the transfert function is the excitation of the white divided by the normed white
+            tf = tf*(ylmi_norm[:, None, None] / normed_white)[None, ...]
+        else:
+            # The transfer function is just the inverse of the white 
+            # (i.e. the excitation spectra of the white is flat)
+            tf = 1 / n_white[None, ...]
 
         # Mask pixel with value lower than the min percentile
         masked_images = dask.array.ma.masked_less_equal(diff_im_darks, np.array(pmin)[None, None, None, :, None, None])
 
         # normalise image stack
-        normed_images = masked_images / n_white[None,...]
+        normed_images = masked_images * tf
 
         # Need to clip images max to 1 !
         normed_images[normed_images > 1.0] = 1
@@ -282,7 +341,26 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
         diff_im_darks = images - darks
         diff_im_darks[diff_im_darks < 0 ] = 0
 
-        normed_images = diff_im_darks / n_white
+        # Compute the transfert function
+        if correct_LM_excitation:
+            # The excitation of the white is not flat
+            print('Use the transfer function from the excitation spectra of the LM')
+            xlm, ylm = load_excitation_LM()
+
+            # Interpolate the excitation spectra of the white
+            ylmi = np.interp(all_white_ex_wavelength, xlm, ylm)
+
+            # Normalise this profile by it's maximum
+            ylmi_norm = ylmi / ylmi.max()
+
+            # the transfert function is the excitation of the white divided by the normed white
+            tf = ylmi_norm / n_white
+        else:
+            # The transfer function is just the inverse of the white 
+            # (i.e. the excitation spectra of the white is flat)
+            tf = 1 / n_white
+
+        normed_images = diff_im_darks * tf
 
         # Need to clip images max to 1 !
         normed_images[normed_images > 1.0] = 1
@@ -295,3 +373,42 @@ def correct_dark_white(images: str, dark: str = '', white: str = '' ,
     # The classical one excitation correction
     else:
         normalizeTELEMOS(imagesin, darkin, whitein, darkofwhitein, white_z, remove_cosmics)
+
+
+def load_excitation_LM():
+    """
+    Measurements of the Lame Matthieu's excitation profile
+
+    Measurement done the 29 septembre 2023 on the Fluoromax4 @ SOLEIL
+    """
+
+    # The wavelength
+    xlm = np.array([240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252,
+                    253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265,
+                    266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278,
+                    279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291,
+                    292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304,
+                    305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317,
+                    318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330,
+                    331, 332, 333, 334, 335, 336, 337, 338, 339, 340])
+
+    # The relative fluorescence measured (emission 530nm, slit 19 nm)
+    ylm = np.array([1.00e+08, 1.01e+08, 1.01e+08, 1.00e+08, 1.00e+08, 9.97e+07,
+                    9.94e+07, 9.88e+07, 9.83e+07, 9.75e+07, 9.71e+07, 9.66e+07,
+                    9.57e+07, 9.51e+07, 9.46e+07, 9.38e+07, 9.30e+07, 9.23e+07,
+                    9.18e+07, 9.13e+07, 9.05e+07, 9.00e+07, 8.95e+07, 8.86e+07,
+                    8.82e+07, 8.74e+07, 8.67e+07, 8.62e+07, 8.54e+07, 8.44e+07,
+                    8.35e+07, 8.22e+07, 8.11e+07, 7.97e+07, 7.81e+07, 7.65e+07,
+                    7.50e+07, 7.32e+07, 7.15e+07, 6.97e+07, 6.81e+07, 6.63e+07,
+                    6.48e+07, 6.31e+07, 6.17e+07, 6.02e+07, 5.87e+07, 5.72e+07,
+                    5.57e+07, 5.37e+07, 5.18e+07, 5.09e+07, 4.95e+07, 4.81e+07,
+                    4.68e+07, 4.55e+07, 4.42e+07, 4.30e+07, 4.19e+07, 4.09e+07,
+                    3.99e+07, 3.90e+07, 3.81e+07, 3.73e+07, 3.66e+07, 3.58e+07,
+                    3.52e+07, 3.45e+07, 3.38e+07, 3.33e+07, 3.27e+07, 3.20e+07,
+                    3.14e+07, 3.09e+07, 3.06e+07, 3.03e+07, 2.99e+07, 2.95e+07,
+                    2.91e+07, 2.86e+07, 2.82e+07, 2.76e+07, 2.71e+07, 2.69e+07,
+                    2.68e+07, 2.66e+07, 2.63e+07, 2.61e+07, 2.59e+07, 2.57e+07,
+                    2.56e+07, 2.55e+07, 2.54e+07, 2.54e+07, 2.53e+07, 2.53e+07,
+                    2.53e+07, 2.53e+07, 2.53e+07, 2.52e+07, 2.51e+07])
+
+    return xlm, ylm
\ No newline at end of file

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "DITB"
-version = "20240311rc13"
+version = "20240519rc1"
 readme = "README.md"
 requires-python = ">=3.11"
 authors = [{ name = "Hugo Chauvet", email = "hugo.chauvet@synchrotron-soleil.fr" }]