fix: issue with empty screenshots

cryoet_data_portal_neuroglancer/precompute/contrast_limits.py

@@ -9,6 +9,8 @@
 from scipy.signal import find_peaks
 from sklearn.cluster import KMeans
 from sklearn.mixture import GaussianMixture
+from scipy.signal import decimate
+
 
 LOGGER = logging.getLogger(__name__)
 
@@ -423,3 +425,100 @@ def plot_cdf(self, output_filename: Optional[str | Path] = None, real_limits: Op
         else:
             plt.show()
         plt.close(fig)
+
+
+
+class SignalDecimationContrastLimitCalculator(ContrastLimitCalculator):
+
+    def __init__(self, volume: Optional["np.ndarray"] = None):
+        """Initialize the contrast limit calculator.
+
+        Parameters
+        ----------
+            volume: np.ndarray or None, optional.
+                Input volume for calculating contrast limits.
+        """
+        super().__init__(volume)
+        self.cdf = None
+        self.limits = None
+        self.decimation = None
+
+    @compute_with_timer
+    def contrast_limits_from_cdf(self) -> tuple[float, float]:
+        """Calculate the contrast limits using the Cumulative Distribution Function.
+
+        Returns
+        -------
+            tuple[float, float]
+                The calculated contrast limits.
+        """
+        # Calculate the histogram of the volume
+        min_value = np.min(self.volume.flatten())
+        max_value = np.max(self.volume.flatten())
+        hist, bin_edges = np.histogram(self.volume.flatten(), bins=400, range=[min_value, max_value])
+
+        # Calculate the CDF of the histogram
+        cdf = np.cumsum(hist) / np.sum(hist)
+        x = np.linspace(min_value, max_value, 400)
+
+        # Downsampling the CDF
+        downsample_factor = 2
+        y_decimated = decimate(cdf, downsample_factor)
+        x_decimated = np.linspace(np.min(x), np.max(x), len(y_decimated))
+
+        # Change threshold based on the original histogram
+        # change_threshold = 0.009 * (max_value - min_value)  # Or use std of the histogram
+
+        # Calculate the absolute differences between consecutive points in the decimated CDF
+        diff_decimated = np.abs(np.diff(y_decimated))
+        # change_threshold = np.max(diff_decimated[:20]) * 2000
+        change_threshold = np.max(diff_decimated[:20]) * 2
+        lower_change_threshold = change_threshold * 1.5
+        # print(diff_decimated)
+
+        initial_flat = np.mean(cdf[:50])  # Average of first 50 points (assumed flat region)
+        final_flat = np.mean(cdf[-50:])   # Average of last 50 points (assumed flat region)
+        midpoint = (initial_flat + final_flat) / 2
+        change_threshold = 0.01 * midpoint
+        lower_change_threshold = change_threshold * 0.1
+        print(change_threshold)
+
+
+        # Detect start and end of slope
+        start_idx_decimated = np.argmax(diff_decimated > change_threshold)  # First large change
+        end_idx_decimated = np.argmax(diff_decimated[start_idx_decimated + 1:] < lower_change_threshold) + start_idx_decimated
+
+        # Check if we found a start index
+        # if start_idx_decimated < len(diff_decimated):
+        #     if end_idx_decimated >= len(diff_decimated):
+        #         end_idx_decimated = len(diff_decimated) - 1  # Clamp to last index
+        # else:
+        #     end_idx_decimated = -1  # No valid end index found
+
+        # Map back the indices to original values
+        self.cdf = [x, cdf]
+        self.limits = (x_decimated[start_idx_decimated], x_decimated[end_idx_decimated]) if end_idx_decimated != -1 else (None, None)
+
+        return self.limits
+
+    def plot_cdf(self, output_filename: Optional[str | Path] = None, real_limits: Optional[list] = None) -> None:
+        """Plot the CDF and the calculated limits."""
+        fig, ax = plt.subplots()
+
+        ax.plot(self.cdf[0], self.cdf[1])
+        ax.axvline(self.limits[0], color="r")
+        ax.axvline(self.limits[1], color="r")
+
+        if real_limits:
+            ax.axvline(real_limits[0], color="b")
+            ax.axvline(real_limits[1], color="b")
+
+        # ax.plot(*self.decimation, "y")
+        # ax.plot(self.cdf[0], self.first_derivative * 100, "y")
+        # ax.plot(self.cdf[0], self.second_derivative * 100, "g")
+
+        if output_filename:
+            fig.savefig(output_filename)
+        else:
+            plt.show()
+        plt.close(fig)

cryoet_data_portal_neuroglancer/state_generator.py

@@ -183,6 +183,8 @@ def combine_json_layers(
     scale = get_scale(scale)
     if projection_quaternion is not None:
         projection_quaternion = Rotation.from_euler(seq="xyz", angles=(45, 0, 0), degrees=True).as_quat()
+    else:
+        projection_quaternion = []
     combined_json = {
         "dimensions": {dim: [res, units] for dim, res in zip("xyz", scale, strict=False)},
         "crossSectionScale": 1.8,

manual_tests/contrast_limits_from_api.py

@@ -1,6 +1,5 @@
 import argparse
 import json
-import logging
 from pathlib import Path
 
 from cryoet_data_portal import Client, Tomogram
@@ -15,18 +14,20 @@
     ContrastLimitCalculator,
     GMMContrastLimitCalculator,
     KMeansContrastLimitCalculator,
+    # SignalDecimationContrastLimitCalculator,
 )
 from cryoet_data_portal_neuroglancer.state_generator import combine_json_layers, generate_image_layer
 
 # Set up logging - level is info
 # logging.basicConfig(level=logging.INFO, force=True)
 
-OUTPUT_FOLDER = "/media/starfish/LargeSSD/data/cryoET/data/FromAPI"
+# OUTPUT_FOLDER = "/media/starfish/LargeSSD/data/cryoET/data/FromAPI"
+OUTPUT_FOLDER = "./COMPARISON"
 
 id_to_path_map = {
     1000: "1000/16.zarr",
-    706: "706/Position_161.zarr",
-    800: "800/0105.zarr",
+    # 706: "706/Position_161.zarr",
+    # 800: "800/0105.zarr",
     10845: "10845/ay18112021_grid2_lamella3_position7.zarr",
     4279: "4279/dga2018-08-27-600.zarr",
 }
@@ -37,11 +38,11 @@
         "volume": [-0.035, 0.009],
         "gain": -7.6,
     },
-    706: {
-        "slice": [-44499.8, 83143],
-        "volume": [-20221.2, 18767.6],
-        "gain": -7.7,
-    },
+    # 706: {
+    #     "slice": [-44499.8, 83143],
+    #     "volume": [-20221.2, 18767.6],
+    #     "gain": -7.7,
+    # },
     800: {
         "slice": [0.0000748111, 0.00189353],
         "volume": [0.000705811, 0.00152511],
@@ -64,10 +65,10 @@
         "https://files.cryoetdataportal.cziscience.com/10008/16/Tomograms/VoxelSpacing14.080/CanonicalTomogram/16.zarr",
         (1.408e-9, 1.408e-9, 1.408e-9),
     ),
-    706: (
-        "https://files.cryoetdataportal.cziscience.com/10004/Position_161/Tomograms/VoxelSpacing7.560/CanonicalTomogram/Position_161.zarr",
-        (7.56e-10, 7.56e-10, 7.56e-10),
-    ),
+    # 706: (
+    #     "https://files.cryoetdataportal.cziscience.com/10004/Position_161/Tomograms/VoxelSpacing7.560/CanonicalTomogram/Position_161.zarr",
+    #     (7.56e-10, 7.56e-10, 7.56e-10),
+    # ),
     800: (
         "https://files.cryoetdataportal.cziscience.com/10005/0105/Tomograms/VoxelSpacing5.224/CanonicalTomogram/0105.zarr",
         (5.224e-10, 5.224e-10, 5.224e-10),
@@ -131,6 +132,12 @@ def run_all_contrast_limit_calculations(id_, input_data_path, output_path):
     limits_dict["cdf"] = limits
     cdf_calculator.plot_cdf(output_path / "cdf.png", real_limits=volume_limit)
 
+    # # Signal decimation based contrast limits
+    # cdf_calculator = SignalDecimationContrastLimitCalculator(calculator.volume)
+    # limits = cdf_calculator.contrast_limits_from_cdf()
+    # limits_dict["decimation"] = limits
+    # cdf_calculator.plot_cdf(output_path / "decimation.png", real_limits=volume_limit)
+
     # 2D contrast limits
     limits = calculator.contrast_limits_from_percentiles(1.0, 99.0)
     limits_dict["wide_percentile"] = limits
@@ -199,7 +206,7 @@ def create_state(id_, contrast_limit_dict, output_folder):
     return json_state
 
 
-def main(output_folder, take_screenshots=False):
+def main(output_folder, take_screenshots=False, wait_for=60 * 1000):
     url_list = []
     for id_, path in id_to_path_map.items():
         path = Path(output_folder) / path
@@ -222,12 +229,13 @@ def main(output_folder, take_screenshots=False):
     if take_screenshots:
         ids = list(id_to_path_map.keys())
         url_dict = {id_: [url] for id_, url in zip(ids, url_list, strict=False)}
-        run_screenshot_loop(url_dict, Path(output_folder))
+        run_screenshot_loop(url_dict, Path(output_folder), wait_for=wait_for)
 
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
     parser.add_argument("--output_folder", type=str, default=OUTPUT_FOLDER)
     parser.add_argument("--screenshot", action="store_true")
+    parser.add_argument("--wait", default=60, type=int, help="How long to wait before taking the screenshot (in s)")
     args, _ = parser.parse_known_args()
-    main(args.output_folder, args.screenshot)
+    main(args.output_folder, args.screenshot, wait_for=(args.wait * 1000))

manual_tests/screenshot_urls.py

@@ -11,39 +11,40 @@ async def take_screenshot(
     urls: dict[int, list[str]],
     output_dir: "Path",
     ids_to_use: list[int] | None = None,
+    wait_for: int = 60 * 1000,
 ):
     output_screenshot_dir = output_dir / "screenshots"
     output_screenshot_dir.mkdir(parents=True, exist_ok=True)
     if ids_to_use is None:
         ids_to_use = list(urls.keys())
     for id in ids_to_use:
         for i, url in enumerate(urls[id]):
-            print(f"Processing url {i} for id {id}")
+            print(f"Processing url {i} for id {id}, waiting for {wait_for / 1000}s before screenshot")
             browser = await launch()
             page = await browser.newPage()
             await page.setViewport(
                 {
                     "width": 1920,
                     "height": 1080,
-                }
+                },
             )
             await page.goto(
                 url,
-                {
-                    "waitUntil": "networkidle0",
-                },
+                {"waitUntil": "networkidle0", "timeout": 0},
             )
+            await page.waitFor(wait_for)
             await page.screenshot(
                 {"path": str(output_screenshot_dir / f"{id}-screenshot-{i}.png")},
             )
             await browser.close()
 
 
-def run_screenshot_loop(urls: dict[int, list[str]], output_dir: "Path", ids_to_use: list[int] | None = None):
+def run_screenshot_loop(
+    urls: dict[int, list[str]],
+    output_dir: "Path",
+    ids_to_use: list[int] | None = None,
+    wait_for: int = 60 * 1000,
+):
     asyncio.get_event_loop().run_until_complete(
-        take_screenshot(
-            urls=urls,
-            output_dir=output_dir,
-            ids_to_use=ids_to_use,
-        )
+        take_screenshot(urls=urls, output_dir=output_dir, ids_to_use=ids_to_use, wait_for=wait_for),
     )