Threshold filter

Author: MaximLippeveld

scip_tiff.yml

@@ -15,6 +15,7 @@ masking:
   bbox_channel_index: 0
   export: false
   kwargs:
+    smooth: 1
 normalization:
   lower: 0
   upper: 1

scip_zarr.yml

@@ -15,6 +15,7 @@ masking:
   export: false
   kwargs:
     noisy_channels: [0]
+filter:
 normalization:
 export:
   format: parquet

src/scip/filters/__init__.py

@@ -0,0 +1,27 @@
+from typing import Mapping
+
+from scip.utils.util import copy_without
+import dask.bag
+import scipy.stats
+
+
+def feature_partition(part):
+    for i in range(len(part)):
+        if "pixels" in part[i]:
+            part[i]["filter_sum"] = part[i]["pixels"][0].sum()
+    return part
+
+
+def item(bag: dask.bag.Bag) -> Mapping[str, dask.bag.Item]:
+    bag = bag.filter(lambda a: "filter_sum" in a).map(lambda a: a["filter_sum"])
+    mu = bag.mean()
+    std = bag.std()
+    return dict(mu=mu, std=std)
+
+
+def predicate(x, *, mu, std):
+    q5 = scipy.stats.norm.ppf(0.05, loc=mu, scale=std)
+    if ("filter_sum" in x) and (x["filter_sum"] > q5):
+        return x
+    else:
+        return copy_without(x, without=["mask", "pixels"])

src/scip/filters/threshold.py

@@ -23,18 +23,46 @@
 from pathlib import Path
 import re
 from typing import Tuple
+import copy
+
+
+def reload_image_partition(
+    partition,
+    channels,
+    clip: int,
+    regex: str,
+    limit: int = -1
+):
+    z = zarr.open(partition[0]["path"])
+    indices = [p["zarr_idx"] for p in partition]
+    data = z.get_coordinate_selection(indices)
+    shapes = numpy.array(z.attrs["shape"])[indices]
+
+    newpartition = copy.deepcopy(partition)
+    for i in range(len(partition)):
+        if "mask" not in partition[i]:
+            continue
+
+        if clip is not None:
+            newpartition[i]["pixels"] = numpy.clip(data[i].reshape(shapes[i])[channels], 0, clip)
+        else:
+            newpartition[i]["pixels"] = data[i].reshape(shapes[i])[channels]
+        newpartition[i]["pixels"] = newpartition[i]["pixels"].astype(numpy.float32)
+    return newpartition
 
 
 def load_image_partition(partition, z, channels, clip):
+
     start, end = partition[0]["zarr_idx"], partition[-1]["zarr_idx"]
     data = z[start:end + 1]
     shapes = z.attrs["shape"][start:end + 1]
-    for i, event in enumerate(partition):
+
+    for i in range(len(partition)):
         if clip is not None:
-            event["pixels"] = numpy.clip(data[i].reshape(shapes[i])[channels], 0, clip)
+            partition[i]["pixels"] = numpy.clip(data[i].reshape(shapes[i])[channels], 0, clip)
         else:
-            event["pixels"] = data[i].reshape(shapes[i])[channels]
-    event["pixels"] = event["pixels"].astype(numpy.float32)
+            partition[i]["pixels"] = data[i].reshape(shapes[i])[channels]
+        partition[i]["pixels"] = partition[i]["pixels"].astype(numpy.float32)
     return partition
 
 

src/scip/loading/zarr.py

@@ -47,10 +47,10 @@ def get_images_bag(
     channels: List[int],
     config: dict,
     partition_size: int,
-    gpu_accelerated: bool
+    gpu_accelerated: bool,
+    loader_module
 ) -> Tuple[dask.bag.Bag, int, dict]:
 
-    loader_module = import_module('scip.loading.%s' % config["loading"]["format"])
     loader = partial(
         loader_module.bag_from_directory,
         channels=channels,
@@ -109,7 +109,7 @@ def get_schema(event):
 
 
 def remove_pixels(event):
-    newevent = copy_without(event, ["pixels", "shape", "mask"])
+    newevent = copy_without(event, ["pixels", "mask"])
     newevent["shape"] = list(event["mask"].shape)
     newevent["mask"] = event["mask"].ravel()
     return newevent
@@ -203,13 +203,15 @@ def main(
 
         logger.debug("loading images in to bags")
 
+        loader_module = import_module('scip.loading.%s' % config["loading"]["format"])
         with dask.config.set(**{'array.slicing.split_large_chunks': False}):
             images, maximum_pixel_value, loader_meta = get_images_bag(
                 paths=paths,
                 channels=channels,
                 config=config,
                 partition_size=partition_size,
-                gpu_accelerated=gpu > 0
+                gpu_accelerated=gpu > 0,
+                loader_module=loader_module
             )
 
         futures = []
@@ -309,6 +311,23 @@ def main(
                     output=output
                 ))
 
+        if config["filter"] is not None:
+            filter_module = import_module('scip.loading.%s' % config["filter"]["name"])
+
+            images = images.map_partitions(filter_module.feature_partition)
+            images = images.map(copy_without, without=["pixels"]).persist()
+            filter_items = filter_module.grouped_item(images, key="group")
+
+            filter_items = filter_items["mu"].compute()
+
+            images = images.map(filter_module.predicate, **filter_items)
+
+            images = images.map_partitions(
+                loader_module.reload_image_partition,
+                channels=channels,
+                **(config["loading"]["loader_kwargs"] or dict())
+            )
+
         quantiles = None
         if config["normalization"] is not None:
             logger.debug("performing normalization")