Tif support

docs/source/Data_processing/Overview.rst

@@ -11,8 +11,10 @@ The package is an image processing pipeline, that can be divided intro three mai
 Use the jupyter notebook ``snazzy-processing-pipeline.ipynb`` to run the pipeline.
 
 The data processing requires a ``.tif`` file.
-There is built in support for formatting ``.nd2`` files to ``.tif``.
-If your raw data is in another format, you must first convert if to ``.tif``.
+There is built in support for converting ``.nd2`` files to ``.tif``.
+This means that you can feed either ``.tif`` or ``.nd2`` files into the pipeline.
+If your raw data is in another format, you must first convert it to ``.tif``.
+ImageJ for example provides several plugins to convert files to tif, including the excellent `BioFormats extension <https://imagej.net/formats/bio-formats>`__.
 
 Before actually running the pipeline, which is the last cell of the jupyter notebook, we must determine from where to crop each movie in the raw data.
 

docs/source/Data_processing/Process_raw_data.rst

@@ -7,6 +7,7 @@ There is a considerable amount of background pixels that can be ignored in the r
 This already saves considerable ROM memory but most importantly, it means we can easily load individual movies in the RAM of a regular computer (8~16 GB RAM), without needing to use memory mapped files.
 
 The algorithm to process the raw image can be resumed as:
+
 1. Get the maximum projection of each pixel for the first 10 frames
 2. Automatic threshold (Triangle method)
 3. Binarize the image 

docs/source/Data_processing/ROIs_and_signal_intensity.rst

@@ -6,6 +6,7 @@ By default, a single ROI is calculated for groups of 10 frames to speed up the p
 This is a good approximation and the speed up justifies the eventual errors in readings caused by movement (see `activity.ipynb` for details about the error in activity caused by downsampling).
 
 The ROI algorithm can be resumed as:
+
 1. Average the group of frames into a single 2D matrix
 2. Automatic threshold (Otsu's method)
 3. Binarize the image 
@@ -14,4 +15,18 @@ The ROI algorithm can be resumed as:
 6. Return a mask that matches the largest label
 
 To calculate the signal intensity, we apply the mask to the embryo and calculate the mean pixel value.
-The active and structural channel measurements are exported as a `.csv` file and further processed using the code from ``snazzy_analysis``.
+The active and structural channel measurements are exported as a `.csv` file and further processed using the code from ``snazzy_analysis``.
+
+Visualizing calculated ROIs
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The ROIs can be inspected visually by running the ``plot_countours.py`` script.
+The script displays a matplotlib animation with an overlayed ROI contour.
+To display it, ``cd`` into the ``snazzy_processing`` directory, and run the file:
+
+.. code:: bash
+
+    python3 scripts/plot_contours.py
+
+It will look for any experiment directories you have inside the ``./data`` directory and present the available options in the terminal.
+Animations can be paused by pressing any key.

environment.yml

@@ -1,4 +1,4 @@
-name: snazzy-env
+name: snazzy-test-env
 channels:
   - conda-forge
   - defaults
@@ -17,18 +17,15 @@ dependencies:
   - pydantic[version='<2']
   - scikit-learn=1.6.1
   - scikit-image=0.25.0
+  - coverage==7.10.2
+  - ipympl==0.9.4
+  - pyqt==6.9.1
+  - pyqtgraph==0.13.7   
+  - pytest==8.3.5
+  - pytest-qt==4.5.0
+  - statannotations==0.7.1
+  - sphinx-book-theme==1.1.2
+  - sphinx==7.3.7
   - pip:
-    - coverage==7.10.2
-    - ipympl==0.9.4
-    - ipython-genutils==0.2.0
-    - pyqt6==6.8.0
-    - pyqt6-qt6==6.8.1
-    - pyqt6-sip==13.9.1
-    - pyqtgraph==0.13.7   
-    - pytest==8.3.5
-    - pytest-qt==4.5.0
-    - statannotations==0.7.1
-    - sphinx-book-theme==1.1.2
-    - sphinx==7.3.7
     - -e ./snazzy_analysis
     - -e ./snazzy_processing

snazzy_analysis/snazzy_analysis/gui/compare_plot_window.py

@@ -187,7 +187,7 @@ def sna_duration(self, save=False, save_dir=None):
             self._save_plot(save_dir, "sna_duration.png")
 
     def plot_AUC(self, save=False, save_dir=None):
-        """Binned area under the curve."""
+        """Area under the curve binned by developmental time."""
         self.clear_axes()
         data = {"group": [], "auc": [], "bin": []}
 
@@ -208,13 +208,15 @@ def plot_AUC(self, save=False, save_dir=None):
                 if not isinstance(bin_idxs, np.ndarray):
                     continue
 
-                data["group"].append(group.name)
+                data["group"].extend([group.name] * len(bin_idxs))
                 data["auc"].extend(trace.peak_aucs)
                 data["bin"].extend(bin_idxs)
 
         bins.append(first_bin + bin_width * n_bins)
         x_labels = [f"{s:.1f}~{e:.1f}" for (s, e) in zip(bins[:-1], bins[1:])]
-        ax = sns.pointplot(data=data, x="bin", y="auc", linestyle="None", ax=self.ax)
+        ax = sns.pointplot(
+            data=data, x="bin", y="auc", hue="group", linestyle="None", ax=self.ax
+        )
         ax.set_xticks(ticks=list(range(n_bins)), labels=x_labels)
         ax.set_title(f"Binned AUC")
         ax.set_ylabel("AUC [activity*t]")

snazzy_analysis/snazzy_analysis/gui/gui.py

@@ -99,6 +99,8 @@ def _show_experiment_dialog(
             **exp_params,
             "dff_strategy": dff_strategy,
         }
+        del dialog_params["acquisition_period"]
+
         self.exp_params_dialog = ExperimentParamsDialog(dialog_params, parent=self)
 
         self.exp_params_dialog.accepted.connect(

snazzy_analysis/snazzy_analysis/gui/image_window.py

@@ -15,7 +15,7 @@
 )
 import pyqtgraph as pg
 
-from snazzy_analysis import Embryo
+from snazzy_analysis import Embryo, utils
 
 
 class ImageWindow(QWidget):
@@ -78,7 +78,10 @@ def init_file_selector(self):
         self.selector_label = QLabel("Select a file:")
         self.combo_box = QComboBox()
         embs_path = self.directory.joinpath("embs")
-        file_names = [str(f) for f in embs_path.iterdir() if "ch1.tif" in f.name]
+        file_names = sorted(
+            [str(f) for f in embs_path.iterdir() if "ch1.tif" in f.name],
+            key=utils.emb_id_from_filename,
+        )
         self.combo_box.addItems(file_names)
         self.open_button = QPushButton("Open Viewer")
 

snazzy_analysis/snazzy_analysis/utils.py

@@ -8,9 +8,22 @@ def split_in_bins(arr: np.ndarray, bins: int):
 
 
 def emb_id(emb: Path | str) -> int:
-    """Retuns the number that identifies a given embryos.
+    """Retun the number that identifies a given embryos.
 
     Assumes that embryos are always named as emb + id, e.g: `emb21`."""
     if isinstance(emb, Path):
         emb = emb.stem
     return int(emb[3:])
+
+
+def emb_id_from_filename(emb_path: str) -> int:
+    """Return the embryo id based on filename.
+
+    Assumes that files are named as embXX-chY.
+
+    Parameters:
+        emb_path (str):
+            Full path to embryo file, as a string.
+    """
+    emb_name = Path(emb_path).stem
+    return int(emb_name.split("-")[0][3:])

snazzy_processing/notebooks/snazzy-processing-pipeline.ipynb

@@ -34,7 +34,7 @@
     "res_dir.mkdir(parents=True, exist_ok=True)\n",
     "\n",
     "# Provide an absolute path for the raw tif image\n",
-    "# If you are starting with another image format, ignore this step and convert the\n",
+    "# If you are starting with another image format, ignore this line and convert the\n",
     "# image in the next cell\n",
     "img_path = ''"
    ]
@@ -43,7 +43,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If the raw data is in nd2 format, it must first be converted to a tif file, and that tif file should be used in the other cells of this jupyter notebook."
+    "If the raw data is in nd2 format, it must first be converted to a tif file, and that tif file should be used in the other cells of this jupyter notebook.\n",
+    "\n",
+    "Adjust the image paths to where the nd2 file is, and where the tif file will be saved."
    ]
   },
   {
@@ -56,14 +58,29 @@
     "# Path where the raw image is located:\n",
     "nd2_path = base_path.joinpath(\"Documents\", \"raw_data\", f\"{experiment_name}.nd2\")\n",
     "# Path where the new tiff file will be saved:\n",
-    "img_path = base_path.joinpath(\"Documents\", \"raw_data\", f\"{experiment_name}.tiff\")\n",
-    "\n",
+    "img_path = base_path.joinpath(\"Documents\", \"raw_data\", f\"{experiment_name}.tif\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Extract the first frames from the image. These frames are used to determine the embryo positions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "first_frames = f\"first_frames.tif\"\n",
     "first_frames_path = root_dir.joinpath(\"results\", experiment_name, first_frames)\n",
     "if first_frames_path.exists():\n",
     "    print(f\"{first_frames_path.stem} already exists.\\nPlease choose another path.\")\n",
     "else:\n",
-    "    slice_img.save_first_frames_as_tiff(nd2_path, first_frames_path, 10)"
+    "    file_path = nd2_path if nd2_path else img_path\n",
+    "    slice_img.save_first_frames_as_tiff(file_path, first_frames_path, 10)"
    ]
   },
   {
@@ -83,7 +100,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "img = slice_img.get_first_image_from_mmap(first_frames_path)\n",
+    "img = slice_img.get_first_image(first_frames_path)\n",
     "\n",
     "coords = slice_img.calculate_slice_coordinates(\n",
     "    first_frames_path, n_cols=4, thres_adjust=-5\n",
@@ -127,6 +144,7 @@
    "source": [
     "If the image above looks good, select which embryos you want to analyze, by changing the values in the `embryos` list in the next cell.\n",
     "Embryos will be saved under the `data` directory, for the corresponding experiment.\n",
+    "To process all embryos you can just pass an empty `embryos` list.\n",
     "\n",
     "The length and activity data will be saved in the `results` directory."
    ]
@@ -141,7 +159,7 @@
     "clean_up_data = False\n",
     "# List of the ids of the embryos that should be processed\n",
     "# The ids are the bbox numbers from the previous cell output\n",
-    "embryos = list(range(1, 4))\n",
+    "embryos = [1, 2, 3, 4, 5]\n",
     "# Interval (number of frames) used to calculate VNC length\n",
     "vnc_length_interval = 10\n",
     "# Window (number of frames) to calculate VNC ROI (which is then used to calculate activity)\n",
@@ -151,9 +169,10 @@
     "embs_dest = root_dir.joinpath(\"data\", experiment_name, \"embs\")\n",
     "embs_dest.mkdir(parents=True, exist_ok=True)\n",
     "\n",
-    "# nd2 to tiff\n",
-    "print(\"Converting from nd2 to tif\")\n",
-    "slice_img.save_as_tiff(nd2_path, img_path)\n",
+    "# nd2 to tif\n",
+    "if nd2_path:\n",
+    "    print(\"Converting from nd2 to tif\")\n",
+    "    slice_img.save_as_tiff(nd2_path, img_path)\n",
     "\n",
     "# tif to individual movies\n",
     "print(\"Cropping individual movies\")\n",

snazzy_processing/snazzy_processing/centerline.py

@@ -128,6 +128,7 @@ def centerline_mask(img_shape: tuple, predictor: RANSACRegressor.predict) -> np.
     # the RANSAC estimations might fall out of the image range
     # make sure that the estimate is within the image dimensions
     rr = np.clip(rr, 0, rows - 1)
+    cc = np.clip(cc, 0, cols - 1)
 
     mask = np.zeros(img_shape, dtype=np.bool_)
     mask[rr, cc] = True
@@ -209,10 +210,11 @@ def view_centerline_dist(binary_image: np.ndarray, ax: Axes, thres_rel=0.6, min_
 
     rr, cc = line(y_start, x_start, y_end, x_end)
 
+    rr = np.clip(rr, 0, rows - 1)
+    cc = np.clip(cc, 0, cols - 1)
     mask = binary_image[rr, cc] == 1
 
     rr_m, cc_m = rr[mask], cc[mask]
-    rr = np.clip(rr, 0, rows - 1)
 
     inliers = estimator.inlier_mask_
     outliers = np.logical_not(inliers)

snazzy_processing/snazzy_processing/slice_img.py

@@ -6,7 +6,7 @@
 from skimage.filters import threshold_triangle
 from skimage.morphology import binary_closing, octagon
 from skimage.exposure import equalize_hist
-from tifffile import imwrite, TiffFile
+from tifffile import imread, imwrite, TiffFile
 import numpy as np
 
 from snazzy_processing import utils
@@ -43,26 +43,30 @@ def save_as_tiff(file: Path, dest_path: Path):
 
 
 def save_first_frames_as_tiff(file: Path, dest_path: Path, n: int):
-    """Save the first frames of an nd2 file as tif.
+    """Save the first frames as tif.
 
     Does not overwrite the file if `dest_path` exists.
 
     Parameters:
         file (Path):
-            Path to nd2 file.
+            Path to tif or nd2 file.
         dest_path (Path):
-            Path to save tiff file.
+            Path to save tif file.
         n (int):
             Number of frames to save.
     """
     dest = Path(dest_path)
     if dest.exists():
         print(f"File '{dest.name}' already exists.")
         return
-    with ND2File(file) as f:
-        darray = f.to_dask()
-        initial_frames = darray[:n].compute()
+    if file.suffix == ".tif" or file.suffix == ".tiff":
+        initial_frames = imread(file, key=slice(0, n))
         imwrite(dest_path, initial_frames)
+    else:
+        with ND2File(file) as f:
+            darray = f.to_dask()
+            initial_frames = darray[:n].compute()
+            imwrite(dest_path, initial_frames)
 
 
 def get_threshold(img: np.ndarray, thres_adjust=0) -> float:
@@ -333,6 +337,7 @@ def cut_movies(
             Directory where the movies will be saved.
         embryos (list[int]):
             List of embryo numbers. Used to select a subgroup of embryos.
+            If not provided, all embryos will be processed.
         active_ch (1 | 2):
             Indicates the image active channel.
             Defaults to 1.
@@ -416,13 +421,11 @@ def get_initial_frames_from_mmap(img_path: Path, n=10):
     return read_mmap(img_path, num_frames=n)
 
 
-def get_first_image_from_mmap(img_path: Path):
-    """Returns the first image from a mmap file, for plotting.
+def get_first_image(img_path: Path):
+    """Returns the first image from for plotting.
 
-    The image is the average of the first 10 slices for channel 2.
-    It is also equalized, since this method is supposed to be used for
-    displaying the image."""
-    img = get_initial_frames_from_mmap(img_path, n=10)
+    The channel 2 frames are averaged and equalized, for better visualization."""
+    img = imread(img_path)
     first_frame = np.average(img[:, 1, :, :], axis=0)
     return equalize_hist(first_frame)