diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9487c5ea6663c6aab78c46aa49dc77c4fe275264..5059b9ce9c584546f3d51568560ad48c79104b8e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -38,7 +38,15 @@ Local Tests:
stage: tests
script:
# - poetry install -vv
- - poetry run pytest ./tests --ignore ./tests/aliby/network --ignore ./tests/aliby/pipeline
+ - poetry run coverage run -m --branch pytest ./tests --ignore ./tests/aliby/network --ignore ./tests/aliby/pipeline
+ - poetry run coverage report -m
+ - poetry run coverage xml
+ coverage: '/(?i)total.*? (100(?:\.0+)?\%|[1-9]?\d(?:\.\d+)?\%)$/'
+ artifacts:
+ reports:
+ coverage_report:
+ coverage_format: cobertura
+ path: coverage.xml
Network Tools Tests:
stage: tests
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 5b051ca2fd056d45496cce9dcb83c6438b1e2af7..b62a8f465d37cbb7b581b901db68b3732da0fa76 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -1,9 +1,9 @@
# Contributing
-We focus our work on python 3.7 due to the current neural network being developed on tensorflow 1. In the near future we will migrate the networ to pytorch to support more recent versions of all packages.
+We focus our work on python 3.8 due to the current neural network being developed on tensorflow 1. In the near future we will migrate the network to pytorch to support more recent versions of all packages.
## Issues
-All issues are managed within the gitlab [ repository ](https://git.ecdf.ed.ac.uk/swain-lab/aliby/aliby/-/issues), if you don't have an account on the University of Edinburgh's gitlab instance and would like to submit issues please get in touch with [Prof. Peter Swain](mailto:peter.swain@ed.ac.uk ).
+All issues are managed within the gitlab [ repository ](https://gitlab.com/aliby/aliby/-/issues), if you don't have an account on the University of Edinburgh's gitlab instance and would like to submit issues please get in touch with [Prof. Peter Swain](mailto:peter.swain@ed.ac.uk ).
## Data aggregation
diff --git a/README.md b/README.md
index 5b5557465d4372bfaacb245645a6aff2543deb96..0837e85a75506fcb2252dfb11a28d13d780ba764 100644
--- a/README.md
+++ b/README.md
@@ -2,47 +2,38 @@
[](https://aliby.readthedocs.io/en/latest)
[](https://badge.fury.io/py/aliby)
-[](https://git.ecdf.ed.ac.uk/swain-lab/aliby/aliby/-/pipelines)
-[](https://git.ecdf.ed.ac.uk/swain-lab/aliby/aliby/-/commits/dev)
+[](https://gitlab.com/aliby/aliby/-/pipelines)
+[](https://gitlab.com/aliby/aliby/-/commits/dev)
+[](https://gitlab.com/aliby/aliby/-/commits/dev)
End-to-end processing of cell microscopy time-lapses. ALIBY automates segmentation, tracking, lineage predictions, post-processing and report production. It leverages the existing Python ecosystem and open-source scientific software available to produce seamless and standardised pipelines.
## Quickstart Documentation
+Installation of [VS Studio](https://visualstudio.microsoft.com/downloads/#microsoft-visual-c-redistributable-for-visual-studio-2022) Native MacOS support for is under work, but you can use containers (e.g., Docker, Podman) in the meantime.
-We use (and recommend) [OMERO](https://www.openmicroscopy.org/omero/) to manage our microscopy database, but ALIBY can process both locally-stored experiments and remote ones hosted on a server.
+For analysing local data
+ ```bash
+pip install aliby # aliby[network] if you want to access an OMERO server
+ ```
-### Setting up a server
-For testing and development, the easiest way to set up an OMERO server is by
-using Docker images.
-[The software carpentry](https://software-carpentry.org/) and the [Open
- Microscopy Environment](https://www.openmicroscopy.org), have provided
-[instructions](https://ome.github.io/training-docker/) to do this.
+See our [installation instructions]( https://aliby.readthedocs.io/en/latest/INSTALL.html ) for more details.
-The `docker-compose.yml` file can be used to create an OMERO server with an
-accompanying PostgreSQL database, and an OMERO web server.
-It is described in detail
-[here](https://ome.github.io/training-docker/12-dockercompose/).
+### CLI
-Our version of the `docker-compose.yml` has been adapted from the above to
-use version 5.6 of OMERO.
+ ```bash
+aliby-run --expt_id EXPT_PATH --distributed 4 --tps None
+ ```
-To start these containers (in background):
-```shell script
-cd pipeline-core
-docker-compose up -d
-```
-Omit the `-d` to run in foreground.
+And to run Omero servers, the basic arguments are shown:
+ ```bash
+ aliby-run --expt_id XXX --host SERVER.ADDRESS --user USER --password PASSWORD
+ ```
-To stop them, in the same directory, run:
-```shell script
-docker-compose stop
-```
+The output is a folder with the original logfiles and a set of hdf5 files, one with the results of each multidimensional inside.
-### Installation
-
-See our [installation instructions]( https://aliby.readthedocs.io/en/latest/INSTALL.html ) for more details.
+## Using specific components
-### Raw data access
+### Access raw data
ALIBY's tooling can also be used as an interface to OMERO servers, taking care of fetching data when needed.
```python
diff --git a/docs/requirements.txt b/docs/requirements.txt
index c44e65a11d24a30475d5a910ed3166660a08ca5f..431210e514a8601bc7af22afbb0b0f8621a6b6e2 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,4 +1,11 @@
numpydoc>=1.3.1
aliby[network]>=0.1.43
+sphinx-autodoc-typehints==1.19.2
+sphinx-rtd-theme==1.0.0
+sphinxcontrib-applehelp==1.0.2
+sphinxcontrib-devhelp==1.0.2
+sphinxcontrib-htmlhelp==2.0.0
+sphinxcontrib-jsmath==1.0.1
+sphinxcontrib-qthelp==1.0.3
+sphinxcontrib-serializinghtml==1.1.5
myst-parser
-sphinx-autodoc-typehints
diff --git a/docs/source/INSTALL.md b/docs/source/INSTALL.md
index 24ff46c5be432271bac78f16a9eb78fcfa911339..94abc570f87bc636041e8d9203c7247a782ccd0b 100644
--- a/docs/source/INSTALL.md
+++ b/docs/source/INSTALL.md
@@ -40,26 +40,71 @@ Or using [pyenv](https://github.com/pyenv/pyenv) with pyenv-virtualenv:
### Pip version
Once you have created and activated your virtual environment, run:
-If you are analysing data locally:
+If you are not using an OMERO server setup:
$ pip install aliby
-If you are contacting an OMERO server:
+Otherwise, if you are contacting an OMERO server:
$ pip install aliby[network]
NOTE: Support for OMERO servers in GNU/Linux computers requires building ZeroC-Ice, thus it requires build tools. The versions for Windows and MacOS are provided as Python wheels and thus installation is faster.
+### FAQ
+- Installation fails during zeroc-ice compilation (Windows and MacOS).
+
+
+For Windows, the simplest way to install it is using conda (or mamba). You can install the (OMERO) network components separately:
+
+ $ conda create -n aliby -c conda-forge python=3.8 omero-py
+ $ conda activate aliby
+ $ cd c:/Users/Public/Repos/aliby
+ $ \PATH\TO\POETRY\LOCATION\poetry install
+
+ - MacOS
+ Under work (See issue https://github.com/ome/omero-py/issues/317)
+
### Git version
-We use [ poetry ](https://python-poetry.org/docs/#installation) for dependency management.
+Install [ poetry ](https://python-poetry.org/docs/#installation) for dependency management.
In case you want to have local version:
- $ git clone git@git.ecdf.ed.ac.uk:swain-lab/aliby/aliby.git
+ $ git clone git@gitlab.com/aliby/aliby.git
$ cd aliby && poetry install --all-extras
-This will automatically install the [ BABY ](https://git.ecdf.ed.ac.uk/swain-lab/aliby/baby) segmentation software. Support for additional segmentation and tracking algorithms is under development.
+This will automatically install the [ BABY ](https://gitlab.com/aliby/baby) segmentation software. Support for additional segmentation and tracking algorithms is under development.
+
+## Omero Server
+
+We use (and recommend) [OMERO](https://www.openmicroscopy.org/omero/) to manage our microscopy database, but ALIBY can process both locally-stored experiments and remote ones hosted on a server.
+
+### Setting up a server
+For testing and development, the easiest way to set up an OMERO server is by
+using Docker images.
+[The software carpentry](https://software-carpentry.org/) and the [Open
+ Microscopy Environment](https://www.openmicroscopy.org), have provided
+[instructions](https://ome.github.io/training-docker/) to do this.
+
+The `docker-compose.yml` file can be used to create an OMERO server with an
+accompanying PostgreSQL database, and an OMERO web server.
+It is described in detail
+[here](https://ome.github.io/training-docker/12-dockercompose/).
+
+Our version of the `docker-compose.yml` has been adapted from the above to
+use version 5.6 of OMERO.
+
+To start these containers (in background):
+```shell script
+cd pipeline-core
+docker-compose up -d
+```
+Omit the `-d` to run in foreground.
+
+To stop them, in the same directory, run:
+```shell script
+docker-compose stop
+```
### Troubleshooting
diff --git a/docs/source/specifications/metadata.org b/docs/source/specifications/metadata.org
index 6618730c2dfb78e133a64cb62a07a9fd8b90c93a..52d514a3465236a9079a8bba3ee721b013808a28 100644
--- a/docs/source/specifications/metadata.org
+++ b/docs/source/specifications/metadata.org
@@ -4,7 +4,7 @@ Draft for recommended metadata for images to provide a standard interface for al
* Essential data
- DimensionOrder: str
- Order of dimensions (e.g., CTZYX for Channel, Time, Z,Y,X)
+ Order of dimensions (e.g., TCZYX for Time, Channel, Z, Y, X)
- PixelSize: float
Size of pixel, useful for segmentation.
- Channels: List[str]
diff --git a/examples/extraction/pos_example.py b/examples/extraction/pos_example.py
index 9372fa1dd1147d7c6584fba695f20238e3dd1167..da44f464b54d81dbba3f5902125d1f7b1fe4abc9 100644
--- a/examples/extraction/pos_example.py
+++ b/examples/extraction/pos_example.py
@@ -14,5 +14,4 @@ params = Parameters(
ext = Extractor(params, omero_id=19310)
-# ext.extract_exp(tile_size=117)
d = ext.extract_tp(tp=1, tile_size=117)
diff --git a/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003.tif b/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003.tif
new file mode 100755
index 0000000000000000000000000000000000000000..569f33a72ab18f65c79e28e38869a4f17cd30cca
Binary files /dev/null and b/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003.tif differ
diff --git a/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003_square.tif b/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003_square.tif
new file mode 100755
index 0000000000000000000000000000000000000000..f813c15a9c8aef8a6e629ff17e460aef4acdb630
Binary files /dev/null and b/examples/tiler/pypipeline_unit_test_00_000001_Brightfield_003_square.tif differ
diff --git a/pyproject.toml b/pyproject.toml
index cb5706e46b26bfb104ace8ad8e34d47d056e59c9..239cb892d0517c669efd9a353bde906ab80e8fe5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
[tool.poetry]
name = "aliby"
-version = "0.1.51"
+version = "0.1.53"
description = "Process and analyse live-cell imaging data"
authors = ["Alan Munoz <alan.munoz@ed.ac.uk>"]
packages = [
@@ -12,6 +12,13 @@ packages = [
]
readme = "README.md"
+[tool.poetry.scripts]
+aliby-run = "aliby.bin.run:run"
+
+[build-system]
+requires = ["setuptools", "poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
+
[tool.poetry.dependencies]
python = ">=3.8, <3.11"
PyYAML = "^6.0"
@@ -26,7 +33,7 @@ py-find-1st = "^1.1.5" # Fast indexing
scikit-learn = ">=1.0.2" # Used for an extraction metric
scipy = ">=1.7.3"
-# [tool.poetry.group.pipeline.dependencies]
+# Pipeline + I/O
aliby-baby = "^0.1.15"
dask = "^2021.12.0"
imageio = "2.8.0" # For image-visualisation utilities
@@ -35,13 +42,13 @@ scikit-image = ">=0.18.1"
tqdm = "^4.62.3" # progress bars
xmltodict = "^0.13.0" # read ome-tiff metadata
-# [tool.poetry.group.postprocessor.dependencies]
+# Postprocessing
leidenalg = "^0.8.8"
more-itertools = "^8.12.0"
pathos = "^0.2.8" # Lambda-friendly multithreading
pycatch22 = "^0.4.2"
-# [tool.poetry.group.network.dependencies]
+# Networking
omero-py = { version = ">=5.6.2", optional = true } # contact omero server
zeroc-ice = { version="3.6.5", optional = true } # networking interface, slow to build
GitPython = "^3.1.27"
@@ -86,10 +93,6 @@ optional = true
pytest = "^6.2.5"
-[build-system]
-requires = ["setuptools", "poetry-core>=1.0.0"]
-build-backend = "poetry.core.masonry.api"
-
[tool.black]
line-length = 79
target-version = ['py38']
diff --git a/src/agora/abc.py b/src/agora/abc.py
index 9c7578092e7ab869e0600ba5bc79abbf3a673f44..c396b4b1503c7558dc51eaf9032b7cb14485bc38 100644
--- a/src/agora/abc.py
+++ b/src/agora/abc.py
@@ -1,13 +1,17 @@
+import logging
import typing as t
from abc import ABC, abstractmethod
from collections.abc import Iterable
from copy import copy
from pathlib import Path, PosixPath
+from time import perf_counter
from typing import Union
from flatten_dict import flatten
from yaml import dump, safe_load
+from agora.logging import timer
+
atomic = t.Union[int, float, str, bool]
@@ -198,6 +202,11 @@ class ProcessABC(ABC):
def run(self):
pass
+ def _log(self, message: str, level: str = "warn"):
+ # Log messages in the corresponding level
+ logger = logging.getLogger("aliby")
+ getattr(logger, level)(f"{self.__class__.__name__}: {message}")
+
def check_type_recursive(val1, val2):
same_types = True
@@ -217,3 +226,28 @@ def check_type_recursive(val1, val2):
for k in val2.keys():
same_types = same_types and check_type_recursive(val1[k], val2[k])
return same_types
+
+
+class StepABC(ProcessABC):
+ """
+ Base class that expands on ProcessABC to include tools used by Aliby steps.
+ It adds a setup step, logging and benchmarking for time benchmarks.
+ """
+
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+
+ @abstractmethod
+ def _run_tp(self):
+ pass
+
+ @timer
+ def run_tp(self, tp: int, **kwargs):
+ """
+ Time and log the timing of a step.
+ """
+ return self._run_tp(tp, **kwargs)
+
+ def run(self):
+ # Replace run withn run_tp
+ raise Warning("Steps use run_tp instead of run")
diff --git a/src/agora/io/bridge.py b/src/agora/io/bridge.py
index 478408f37d1f0724a232876f0f5ce0da5b02b585..788db5e9d7cc24b1a19c261d80fdf9cc74589427 100644
--- a/src/agora/io/bridge.py
+++ b/src/agora/io/bridge.py
@@ -2,9 +2,10 @@
Tools to interact with h5 files and handle data consistently.
"""
import collections
+import logging
+import typing as t
from itertools import chain, groupby, product
from typing import Union
-import typing as t
import h5py
import numpy as np
@@ -25,6 +26,11 @@ class BridgeH5:
self._hdf = h5py.File(filename, flag)
self._filecheck
+ def _log(self, message: str, level: str = "warn"):
+ # Log messages in the corresponding level
+ logger = logging.getLogger("aliby")
+ getattr(logger, level)(f"{self.__class__.__name__}: {message}")
+
def _filecheck(self):
assert "cell_info" in self._hdf, "Invalid file. No 'cell_info' found."
@@ -151,12 +157,6 @@ def attrs_from_h5(fpath: str):
return dict(f.attrs)
-def parameters_from_h5(fpath: str):
- """Return parameters from an h5 file."""
- attrs = attrs_from_h5(fpath)
- return yaml.safe_load(attrs["parameters"])
-
-
def image_creds_from_h5(fpath: str):
"""Return image id and server credentials from an h5."""
attrs = attrs_from_h5(fpath)
diff --git a/src/agora/io/cells.py b/src/agora/io/cells.py
index d110869142a071d5d3f054b6a1e172b535e297b0..35dfbc5bf208793fc702d89431cd631e5248858d 100644
--- a/src/agora/io/cells.py
+++ b/src/agora/io/cells.py
@@ -36,6 +36,11 @@ class Cells:
def from_source(cls, source: t.Union[PosixPath, str]):
return cls(Path(source))
+ def _log(self, message: str, level: str = "warn"):
+ # Log messages in the corresponding level
+ logger = logging.getLogger("aliby")
+ getattr(logger, level)(f"{self.__class__.__name__}: {message}")
+
@staticmethod
def _asdense(array: np.ndarray):
if not isdense(array):
@@ -300,7 +305,7 @@ class Cells:
)
else:
mothers_daughters = np.array([])
- # print("Warning:Cells: No mother-daughters assigned")
+ self._log("No mother-daughters assigned")
return mothers_daughters
diff --git a/src/agora/io/metadata.py b/src/agora/io/metadata.py
index 7fde2b10d03ea2dfe2bfeace981c9a6345a6f700..c058e9f430672a9fc4cb6a7e25d4e328368a2044 100644
--- a/src/agora/io/metadata.py
+++ b/src/agora/io/metadata.py
@@ -85,7 +85,12 @@ def datetime_to_timestamp(time, locale="Europe/London"):
def find_file(root_dir, regex):
- file = glob.glob(os.path.join(str(root_dir), regex))
+ file = [
+ f
+ for f in glob.glob(os.path.join(str(root_dir), regex))
+ if Path(f).name != "aliby.log" # Skip filename reserved for aliby
+ ]
+
if len(file) > 1:
print(
"Warning:Metadata: More than one logfile found. Defaulting to first option."
@@ -163,8 +168,9 @@ def get_meta_swainlab(parsed_metadata: dict):
def get_meta_from_legacy(parsed_metadata: dict):
- channels = parsed_metadata["channels/channel"]
- return {"channels": channels}
+ result = parsed_metadata
+ result["channels"] = result["channels/channel"]
+ return result
def parse_swainlab_metadata(filedir: t.Union[str, PosixPath]):
diff --git a/src/agora/io/signal.py b/src/agora/io/signal.py
index 84d6f049418cb0f10dd97874cc22d18a0e1c219e..c0c3e840f171d8c092ffec03ee4efb0d841f5a92 100644
--- a/src/agora/io/signal.py
+++ b/src/agora/io/signal.py
@@ -42,21 +42,15 @@ class Signal(BridgeH5):
"Cy5",
"pHluorin405",
)
- # Alan: why "equivalences"? this variable is unused.
- equivalences = {
- "m5m": ("extraction/GFP/max/max5px", "extraction/GFP/max/median")
- }
def __getitem__(self, dsets: t.Union[str, t.Collection]):
"""Get and potentially pre-process data from h5 file and return as a dataframe."""
- if isinstance(dsets, str):
- # no pre-processing
+ if isinstance(dsets, str): # no pre-processing
df = self.get_raw(dsets)
return self.add_name(df, dsets)
- elif isinstance(dsets, list):
- # pre-processing
+ elif isinstance(dsets, list): # pre-processing
is_bgd = [dset.endswith("imBackground") for dset in dsets]
- # Alan: what does this error message mean?
+ # Check we are not comaring tile-indexed and cell-indexed data
assert sum(is_bgd) == 0 or sum(is_bgd) == len(
dsets
), "Tile data and cell data can't be mixed"
@@ -72,18 +66,12 @@ class Signal(BridgeH5):
df.name = name
return df
- # def cols_in_mins_old(self, df: pd.DataFrame):
- # """Convert numerical columns in a dataframe to minutes."""
- # try:
- # df.columns = (df.columns * self.tinterval // 60).astype(int)
- # except Exception as e:
- # print(f"Warning:Signal: Unable to convert columns to minutes: {e}")
- # return df
-
def cols_in_mins(self, df: pd.DataFrame):
- """Convert numerical columns in a dataframe to minutes."""
- if type(self.tinterval) == int:
- df.columns *= self.tinterval
+ # Convert numerical columns in a dataframe to minutes
+ try:
+ df.columns = (df.columns * self.tinterval // 60).astype(int)
+ except Exception as e:
+ self._log(f"Unable to convert columns to minutes: {e}", "debug")
return df
@cached_property
@@ -239,7 +227,8 @@ class Signal(BridgeH5):
with h5py.File(self.filename, "r") as f:
f.visititems(self.store_signal_path)
except Exception as e:
- print("Error visiting h5: {}".format(e))
+ self._log("Exception when visiting h5: {}".format(e), "exception")
+
return self._available
def get_merged(self, dataset):
@@ -267,10 +256,10 @@ class Signal(BridgeH5):
def get_raw(
self,
- dataset: str,
+ dataset: str or t.List[str],
in_minutes: bool = True,
lineage: bool = False,
- ) -> pd.DataFrame:
+ ) -> pd.DataFrame or t.List[pd.DataFrame]:
"""
Load data from a h5 file and return as a dataframe.
@@ -289,10 +278,11 @@ class Signal(BridgeH5):
if in_minutes:
df = self.cols_in_mins(df)
elif isinstance(dataset, list):
- # Alan: no mother_labels in this case?
- return [self.get_raw(dset) for dset in dataset]
- if lineage:
- # assumes that df is sorted
+ return [
+ self.get_raw(dset, in_minutes=in_minutes, lineage=lineage)
+ for dset in dataset
+ ]
+ if lineage: # assume that df is sorted
mother_label = np.zeros(len(df), dtype=int)
lineage = self.lineage()
a, b = validate_association(
@@ -304,7 +294,7 @@ class Signal(BridgeH5):
df = add_index_levels(df, {"mother_label": mother_label})
return df
except Exception as e:
- print(f"Could not fetch data set {dataset}")
+ self._log(f"Could not fetch dataset {dataset}: {e}", "error")
raise e
def get_merges(self):
@@ -355,7 +345,10 @@ class Signal(BridgeH5):
fullname: str,
node: t.Union[h5py.Dataset, h5py.Group],
):
- """Store the name of a signal if it is a leaf node (a group with no more groups inside) and if it starts with extraction."""
+ """
+ Store the name of a signal if it is a leaf node
+ (a group with no more groups inside) and if it starts with extraction.
+ """
if isinstance(node, h5py.Group) and np.all(
[isinstance(x, h5py.Dataset) for x in node.values()]
):
diff --git a/src/agora/io/utils.py b/src/agora/io/utils.py
index 0acca82cb57fbab596d990a9e9d554f0c8b80344..b32b9314c3291c66f8b4c01336c25db236308e00 100644
--- a/src/agora/io/utils.py
+++ b/src/agora/io/utils.py
@@ -35,29 +35,6 @@ def imread(path):
return cv2.imread(str(path), -1)
-class ImageCache:
- """HDF5-based image cache for faster loading of the images once they've
- been read.
- """
-
- def __init__(self, file, name, shape, remote_fn):
- self.store = h5py.File(file, "a")
- # Create a dataset
- self.dataset = self.store.create_dataset(
- name, shape, dtype=np.float, fill_value=np.nan
- )
- self.remote_fn = remote_fn
-
- def __getitem__(self, item):
- cached = self.dataset[item]
- if np.any(np.isnan(cached)):
- full = self.remote_fn(item)
- self.dataset[item] = full
- return full
- else:
- return cached
-
-
class Cache:
"""
Fixed-length mapping to use as a cache.
diff --git a/src/agora/io/writer.py b/src/agora/io/writer.py
index f11827fa9adb6818e354c0574f91a419770f8f9a..264478482875c6992241ed9ad00613f0ef96fbf0 100644
--- a/src/agora/io/writer.py
+++ b/src/agora/io/writer.py
@@ -1,4 +1,3 @@
-import itertools
import logging
from collections.abc import Iterable
from pathlib import Path
@@ -55,6 +54,11 @@ class DynamicWriter:
if Path(file).exists():
self.metadata = load_attributes(file)
+ def _log(self, message: str, level: str = "warn"):
+ # Log messages in the corresponding level
+ logger = logging.getLogger("aliby")
+ getattr(logger, level)(f"{self.__class__.__name__}: {message}")
+
def _append(self, data, key, hgroup):
"""
Append data to existing dataset in the h5 file otherwise create a new one.
@@ -135,10 +139,6 @@ class DynamicWriter:
# write all data, signified by the empty tuple
hgroup[key][()] = data
- # def _check_key(self, key):
- # if key not in self.datatypes:
- # raise KeyError(f"No defined data type for key {key}")
-
def write(self, data: dict, overwrite: list, meta: dict = {}):
"""
Write data and metadata to h5 file.
@@ -174,7 +174,9 @@ class DynamicWriter:
self._append(value, key, hgroup)
except Exception as e:
print(key, value)
- raise (e)
+ self._log(
+ f"{key}:{value} could not be written: {e}", "error"
+ )
# write metadata
for key, value in meta.items():
hgroup.attrs[key] = value
@@ -222,220 +224,6 @@ class TilerWriter(DynamicWriter):
super().write(data=data, overwrite=overwrite, meta=meta)
-# Alan: we use complex numbers because...
-@timed()
-def save_complex(array, dataset):
- # append array, an 1D array of complex numbers, onto dataset, a 2D array of real numbers
- n = len(array)
- if n > 0:
- dataset.resize(dataset.shape[0] + n, axis=0)
- dataset[-n:, 0] = array.real
- dataset[-n:, 1] = array.imag
-
-
-@timed()
-def load_complex(dataset):
- # convert 2D dataset into a 1D array of complex numbers
- array = dataset[:, 0] + 1j * dataset[:, 1]
- return array
-
-
-class BabyWriter(DynamicWriter):
- """
- Write data stored in a Baby instance to h5 files.
-
- Assumes the edgemasks are of form ((max_ncells, max_tps, tile_size, tile_size), bool).
- """
-
- compression = "gzip"
- max_ncells = 2e5 # Alan: Could just make this None
- max_tps = 1e3 # Could just make this None
- # the number of cells in a chunk for edge masks
- chunk_cells = 25
- default_tile_size = 117
- datatypes = {
- "centres": ((None, 2), np.uint16),
- "position": ((None,), np.uint16),
- "angles": ((None,), h5py.vlen_dtype(np.float32)),
- "radii": ((None,), h5py.vlen_dtype(np.float32)),
- "edgemasks": (
- (max_ncells, max_tps, default_tile_size, default_tile_size),
- bool,
- ),
- "ellipse_dims": ((None, 2), np.float32),
- "cell_label": ((None,), np.uint16),
- "trap": ((None,), np.uint16),
- "timepoint": ((None,), np.uint16),
- # "mother_assign": ((None,), h5py.vlen_dtype(np.uint16)),
- "mother_assign_dynamic": ((None,), np.uint16),
- "volumes": ((None,), np.float32),
- }
- group = "cell_info"
-
- def __init__(self, *args, **kwargs):
- super().__init__(*args, **kwargs)
- self._traps_initialised = False
-
- def __init_trap_info(self):
- # requires traps to have been initialised
- trap_metadata = load_attributes(self.file, "trap_info")
- tile_size = trap_metadata.get("tile_size", self.default_tile_size)
- max_tps = self.metadata["time_settings/ntimepoints"][0]
- self.datatypes["edgemasks"] = (
- (self.max_ncells, max_tps, tile_size, tile_size),
- bool,
- )
- self._traps_initialised = True
-
- def __init_edgemasks(self, hgroup, edgemasks, current_indices, n_cells):
- # create the values dataset in the h5 file
- # holds the edge masks and has shape (n_tps, n_cells, tile_size, tile_size)
- key = "edgemasks"
- max_shape, dtype = self.datatypes[key]
- shape = (n_cells, 1) + max_shape[2:]
- chunks = (self.chunk_cells, 1) + max_shape[2:]
- val_dset = hgroup.create_dataset(
- "values",
- shape=shape,
- maxshape=max_shape,
- dtype=dtype,
- chunks=chunks,
- compression=self.compression,
- )
- val_dset[:, 0] = edgemasks
- # create index dataset in the h5 file:
- # the (trap, cell_id) description used to index into the values and has shape (n_cells, 2)
- ix_max_shape = (max_shape[0], 2)
- ix_shape = (0, 2)
- ix_dtype = np.uint16
- ix_dset = hgroup.create_dataset(
- "indices",
- shape=ix_shape,
- maxshape=ix_max_shape,
- dtype=ix_dtype,
- compression=self.compression,
- )
- save_complex(current_indices, ix_dset)
-
- def __append_edgemasks(self, hgroup, edgemasks, current_indices):
- val_dset = hgroup["values"]
- ix_dset = hgroup["indices"]
- existing_indices = load_complex(ix_dset)
- # check if there are any new labels
- available = np.in1d(current_indices, existing_indices)
- missing = current_indices[~available]
- all_indices = np.concatenate([existing_indices, missing])
- # resizing
- debug_t = perf_counter() # for timing code for debugging
- n_tps = val_dset.shape[1] + 1
- n_add_cells = len(missing)
- # resize dataset for Time and Cells
- new_shape = (val_dset.shape[0] + n_add_cells, n_tps) + val_dset.shape[
- 2:
- ]
- val_dset.resize(new_shape)
- logging.debug(f"Timing:resizing:{perf_counter() - debug_t}")
- # write data
- cell_indices = np.where(np.in1d(all_indices, current_indices))[0]
- for ix, mask in zip(cell_indices, edgemasks):
- try:
- val_dset[ix, n_tps - 1] = mask
- except Exception as e:
- logging.debug(
- "Exception: {}:{}, {}, {}".format(
- e, ix, n_tps, val_dset.shape
- )
- )
- # save the index values
- save_complex(missing, ix_dset)
-
- def write_edgemasks(self, data, keys, hgroup):
- """
- Write edgemasks to h5 file.
-
- Parameters
- ----------
- data: list of arrays
- Data to be written, in the form (trap_ids, cell_labels, edgemasks)
- keys: list of str
- Names corresponding to the elements of data.
- For example: ["trap", "cell_label", "edgemasks"]
- hgroup: group object
- Group to write to in h5 file.
- """
- if not self._traps_initialised:
- self.__init_trap_info()
- key = "edgemasks"
- val_key = "values"
- traps, cell_labels, edgemasks = data
- n_cells = len(cell_labels)
- hgroup = hgroup.require_group(key)
- # create complex indices with traps as real part and cell_labels as imaginary part
- current_indices = np.array(traps) + 1j * np.array(cell_labels)
- if val_key not in hgroup:
- self.__init_edgemasks(hgroup, edgemasks, current_indices, n_cells)
- else:
- self.__append_edgemasks(hgroup, edgemasks, current_indices)
-
- def write(
- self, data: dict, overwrite: list, tp: int = None, meta: dict = {}
- ):
- """
- Write data from a Baby instance, including edgemasks.
-
- Parameters
- ----------
- data: dict
- A dict of datasets and data
- overwrite: list of str
- A list of datasets to overwrite
- tp: int
- The time point of interest
- meta: dict, optional
- Metadata to be written as attributes of the h5 file
- """
- with h5py.File(self.file, "a") as store:
- hgroup = store.require_group(self.group)
- # write data
- for key, value in data.items():
- if key not in self.datatypes:
- raise KeyError(
- f"BabyWriter: No defined data type for key {key}"
- )
- else:
- try:
- if key.startswith("attrs/"):
- # metadata
- key = key.split("/")[1]
- hgroup.attrs[key] = value
- elif key in overwrite:
- # delete and replace existing dataset
- self._overwrite(value, key, hgroup)
- elif key == "edgemasks":
- keys = ["trap", "cell_label", "edgemasks"]
- value = [data[x] for x in keys]
- edgemask_dset = hgroup.get(key + "/values", None)
- if (
- edgemask_dset
- and tp < edgemask_dset[()].shape[1]
- ):
- # data already exists
- print(
- f"BabyWriter: Skipping edgemasks in tp {tp}"
- )
- else:
- self.write_edgemasks(value, keys, hgroup)
- else:
- # append or create new dataset
- self._append(value, key, hgroup)
- except Exception as e:
- print(key, value)
- raise (e)
- # write metadata
- for key, value in meta.items():
- hgroup.attrs[key] = value
-
-
class LinearBabyWriter(DynamicWriter):
"""
Write data stored in a Baby instance to h5 files.
@@ -456,7 +244,6 @@ class LinearBabyWriter(DynamicWriter):
"cell_label": ((None,), np.uint16),
"trap": ((None,), np.uint16),
"timepoint": ((None,), np.uint16),
- # "mother_assign": ((None,), h5py.vlen_dtype(np.uint16)),
"mother_assign_dynamic": ((None,), np.uint16),
"volumes": ((None,), np.float32),
}
@@ -667,14 +454,6 @@ class Writer(BridgeH5):
if overwrite == "overwrite": # TODO refactor overwriting
if path in f:
del f[path]
- # elif overwrite == "accumulate": # Add a number if needed
- # if path in f:
- # parent, name = path.rsplit("/", maxsplit=1)
- # n = sum([x.startswith(name) for x in f[path]])
- # path = path + str(n).zfill(3)
- # elif overwrite == "skip":
- # if path in f:
- # logging.debug("Skipping dataset {}".format(path))
logging.debug(
"{} {} to {} and {} metadata fields".format(
overwrite, type(data), path, len(meta)
@@ -791,7 +570,7 @@ class Writer(BridgeH5):
dset[()] = df.index.get_level_values(level=name).tolist()
# create dataset and write columns
if (
- df.columns.dtype == np.int
+ df.columns.dtype == int
or df.columns.dtype == np.dtype("uint")
or df.columns.name == "timepoint"
):
@@ -828,9 +607,7 @@ class Writer(BridgeH5):
)
# split indices in existing and additional
new = df.index.tolist()
- if (
- df.index.nlevels == 1
- ):
+ if df.index.nlevels == 1:
# cover cases with a single index
new = [(x,) for x in df.index.tolist()]
(
@@ -931,8 +708,6 @@ class Writer(BridgeH5):
return existing_cells, new_cells
-
-
def locate_indices(existing, new):
if new.any():
if new.shape[1] > 1:
diff --git a/src/agora/logging.py b/src/agora/logging.py
new file mode 100644
index 0000000000000000000000000000000000000000..a004b0182363d495c32ffcb2edf9d51cd8875303
--- /dev/null
+++ b/src/agora/logging.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env jupyter
+"""
+Add general logging functions and decorators
+"""
+
+import logging
+from time import perf_counter
+
+
+def timer(func):
+ # Log duration of a function into aliby logfile
+ def wrap_func(*args, **kwargs):
+ t1 = perf_counter()
+ result = func(*args, **kwargs)
+ logging.getLogger("aliby").debug(
+ f"{func.__qualname__} took {(perf_counter()-t1):.4f}s"
+ )
+ return result
+
+ return wrap_func
diff --git a/src/agora/track_abc.py b/src/agora/track_abc.py
new file mode 100644
index 0000000000000000000000000000000000000000..0a11799eb8531a0ea249455b5a92042d8caa5c99
--- /dev/null
+++ b/src/agora/track_abc.py
@@ -0,0 +1,273 @@
+#!/usr/bin/env jupyter
+
+import typing as t
+from abc import ABC
+from os.path import join
+
+import numpy as np
+from skimage.measure import regionprops_table
+
+from aliby.track.utils import calc_barycentre, pick_baryfun
+
+
+class FeatureCalculator(ABC):
+ """
+ Base class for making use of regionprops-based features.
+ If no features are offered it uses most of them.
+
+ This class is not to be used directly
+ """
+
+ a_ind = None
+ ma_ind = None
+ x_ind = None
+ y_ind = None
+
+ def __init__(
+ self,
+ feats2use: t.Collection[str],
+ trapfeats: t.Optional[t.Collection[str]] = None,
+ extrafeats: t.Optional[t.Collection[str]] = None,
+ aweights: t.Optional[bool] = None,
+ pixel_size: t.Optional[float] = None,
+ ) -> None:
+
+ self.feats2use = feats2use
+
+ if trapfeats is None:
+ trapfeats = ()
+ self.trapfeats = tuple(trapfeats)
+
+ if extrafeats is None:
+ extrafeats = ()
+ self.extrafeats = tuple(extrafeats)
+
+ if aweights is None:
+ aweights = None
+ self.aweights = aweights
+
+ if pixel_size is None:
+ pixel_size = 0.182
+ self.pixel_size = pixel_size
+
+ self.outfeats = self.get_outfeats()
+
+ self.set_named_ids()
+
+ self.tfeats = self.outfeats + self.tmp_outfeats + self.trapfeats
+ self.ntfeats = len(self.tfeats)
+
+ def get_outfeats(
+ self, feats2use: t.Optional[t.Collection[str]] = None
+ ) -> tuple:
+ if feats2use is None:
+ feats2use = self.feats2use
+ outfeats = tuple(
+ regionprops_table(np.diag((1, 0)), properties=feats2use).keys()
+ )
+ return outfeats
+
+ def set_named_ids(self):
+ # Manage calling feature outputs by name
+ d = {"centroid-0": "xind", "centroid-1": "yind", "area": "aind"}
+ tmp_d = {
+ "barydist": ["centroid", "area"],
+ "baryangle": ["centroid", "area"],
+ "distance": ["centroid"],
+ }
+
+ nonbase_feats = self.trapfeats + self.extrafeats
+
+ tmp_infeats = np.unique([j for x in nonbase_feats for j in tmp_d[x]])
+ self.tmp_infeats = tuple(
+ [f for f in tmp_infeats if f not in self.feats2use]
+ )
+
+ # feats that are only used to calculate others
+ tmp_outfeats = (
+ self.get_outfeats(feats2use=tmp_infeats)
+ if len(tmp_infeats)
+ else []
+ )
+
+ self.tmp_outfeats = []
+ for feat in tmp_outfeats:
+ if feat in self.outfeats:
+ setattr(self, d[feat], self.outfeats.index(feat))
+ else: # Only add them if not in normal outfeats
+ self.tmp_outfeats.append(feat)
+ setattr(
+ self,
+ d[feat],
+ len(self.outfeats) + self.tmp_outfeats.index(feat),
+ )
+
+ self.tmp_outfeats = tuple(self.tmp_outfeats)
+ self.out_merged = self.outfeats + self.tmp_outfeats
+
+ def load_model(self, path, fname):
+ model_file = join(path, fname)
+ with open(model_file, "rb") as file_to_load:
+ model = pickle.load(file_to_load)
+
+ return model
+
+ def calc_feats_from_mask(
+ self,
+ masks: np.ndarray,
+ feats2use: t.Optional[t.Tuple[str]] = None,
+ trapfeats: t.Optional[t.Tuple[str]] = None,
+ scale: t.Optional[bool] = True,
+ pixel_size: t.Optional[float] = None,
+ ):
+ """
+ Calculate feature ndarray from ndarray of cell masks
+ ---
+ input
+
+ :masks: ndarray (ncells, x_size, y_size), typically dtype bool
+ :feats2use: list of strings with the feature properties to extract.
+ If it is None it uses the ones set in self.feats2use.
+ :trapfeats: List of str with additional features to use
+ calculated immediately after basic features.
+ :scale: bool, if True scales mask to a defined pixel_size.
+ :pixel_size: float, used to rescale the object features.
+
+ returns
+
+ (ncells, nfeats) ndarray of features for input masks
+ """
+ if pixel_size is None:
+ pixel_size = self.pixel_size
+
+ if feats2use is None:
+ feats2use = self.feats2use + self.tmp_infeats
+
+ if trapfeats is None:
+ trapfeats = self.trapfeats
+
+ ncells = masks.shape[0] if masks.ndim == 3 else masks.max()
+ feats = np.empty((ncells, self.ntfeats)) # ncells * nfeats
+ if masks.any():
+ if masks.ndim == 3: # Individual cells in dim 0
+ assert masks.sum(
+ axis=(1, 2)
+ ).all(), "Dimension with at least one empty outline slice"
+
+ cell_feats = np.array(
+ [
+ [
+ x[0]
+ for x in regionprops_table(
+ mask.astype(int), properties=feats2use
+ ).values()
+ ]
+ for mask in masks
+ ]
+ )
+
+ elif masks.ndim == 2: # No overlap between cells
+ cell_feats = np.array(
+ [
+ x
+ for x in regionprops_table(
+ masks.astype(int), properties=feats2use
+ ).values()
+ ]
+ ).T
+ else:
+ raise Exception(
+ "TrackerException: masks do not have the appropiate dimensions"
+ )
+
+ if scale:
+ cell_feats = self.scale_feats(cell_feats, pixel_size)
+
+ # Fill first sector, with directly extracted features
+ feats[:, : len(self.out_merged)] = cell_feats
+ if trapfeats: # Add additional features
+
+ tfeats = self.calc_trapfeats(feats)
+ feats[:, len(self.out_merged) :] = tfeats
+ else:
+
+ feats = np.zeros((0, self.ntfeats))
+
+ return feats
+
+ def calc_trapfeats(self, basefeats):
+ """
+ Calculate trap-based features
+ using basic ones.
+ :basefeats: (n basic outfeats) 1-D array with features outputed by
+ skimage.measure.regionprops_table
+
+
+ requires
+ self.aind
+ self.aweights
+ self.xind
+ self.yind
+ self.trapfeats
+
+ returns
+ (ntrapfeats) 1-D array with
+ """
+ if self.aweights is not None:
+ weights = basefeats[:, self.aind]
+ else:
+ weights = None
+
+ barycentre = calc_barycentre(
+ basefeats[:, [self.xind, self.yind]], weights=weights
+ )
+
+ trapfeat_nd = np.empty((basefeats.shape[0], len(self.trapfeats)))
+ for i, trapfeat in enumerate(self.trapfeats):
+ trapfeat_nd[:, i] = pick_baryfun(trapfeat)(
+ basefeats[:, [self.xind, self.yind]], barycentre
+ )
+
+ return trapfeat_nd
+
+ def scale_feats(self, feats: np.ndarray, pixel_size: float):
+ """
+ input
+
+ :feats: np.ndarray (ncells * nfeatures)
+ :pixel_size: float Value used to normalise the images.
+
+ returns
+ Rescaled list of feature values
+ """
+ area = pixel_size**2
+
+ scaling = {None: 1, "linear": pixel_size, "square": area}
+ degrees_feats = {
+ None: ["eccentricity", "extent", "orientation", "solidity"],
+ "linear": [
+ "centroid-0",
+ "centroid-1",
+ "minor_axis_length",
+ "major_axis_length",
+ "perimeter",
+ "perimeter_crofton",
+ "equivalent_diameter",
+ ],
+ "square": [
+ "area",
+ "convex_area",
+ "bbox_area",
+ "equivalent_diameter",
+ ],
+ }
+
+ scaler = {
+ feat: scaling[k]
+ for k, feats in degrees_feats.items()
+ for feat in feats
+ }
+ # for k in feats.keys():
+ # feats[k] /= scaler[k]
+
+ return feats * [scaler[feat] for feat in self.out_merged]
diff --git a/src/agora/utils/example.py b/src/agora/utils/example.py
deleted file mode 100644
index e3ff571acc0cb780e2feb92007111dbe041a91fb..0000000000000000000000000000000000000000
--- a/src/agora/utils/example.py
+++ /dev/null
@@ -1,53 +0,0 @@
-"""This is an example module to show the structure."""
-from typing import Union
-
-
-class ExampleClass:
- """This is an example class to show the structure."""
-
- def __init__(self, parameter: int):
- """This class takes one parameter and is used to add one to that
- parameter.
-
- :param parameter: The parameter for this class
- """
- self.parameter = parameter
-
- def add_one(self):
- """Takes the parameter and adds one.
-
- >>> x = ExampleClass(1)
- >>> x.add_one()
- 2
-
- :return: the parameter + 1
- """
- return self.parameter + 1
-
- def add_n(self, n: int):
- """Adds n to the class instance's parameter.
-
- For instance
- >>> x = ExampleClass(1)
- >>> x.add_n(10)
- 11
-
- :param n: The number to add
- :return: the parameter + n
- """
- return self.parameter + n
-
-
-def example_function(parameter: Union[int, str]):
- """This is a factory function for an ExampleClass.
-
- :param parameter: the parameter to give to the example class
- :return: An example class
- """
- try:
- return ExampleClass(int(parameter))
- except ValueError as e:
- raise ValueError(
- f"The parameter {parameter} could not be turned "
- f"into an integer."
- ) from e
diff --git a/src/agora/utils/lineage.py b/src/agora/utils/lineage.py
index 5b6686863f0262e515a6164db29c17dbecd80920..52fb552bd530e872ec013b119040cc5d4bba8764 100644
--- a/src/agora/utils/lineage.py
+++ b/src/agora/utils/lineage.py
@@ -20,60 +20,3 @@ def mb_array_to_dict(mb_array: np.ndarray):
for mo, daughters in groupsort(mo_da).items()
}
-
-def mb_array_to_indices(mb_array: np.ndarray):
- """
- Convert a lineage ndarray (trap, mother_id, daughter_id)
- into a dictionary of lists ( mother_id ->[daughters_ids] )
- """
- return pd.MultiIndex.from_arrays(mb_array[:, :2].T).union(
- pd.MultiIndex.from_arrays(mb_array[:, [0, 2]].T)
- )
-
-
-def group_matrix(
- matrix: np.ndarray,
- n_keys: int = 2,
-) -> t.Dict[t.Tuple[int], t.List[int]]:
- """Group a matrix of integers by grouping the first two columns
- and setting the third one in a list.
-
-
- Parameters
- ----------
- matrix : np.ndarray
- id_matrix, generally its columns are three integers indicating trap,
- mother and daughter.
- n_keys : int
- number of keys to use to determine groups.
-
- Returns
- -------
- t.Dict[t.Tuple[int], t.Collection[int, ...]]
- The column(s) not used for generaeting keys are grouped as values.
-
- Examples
- --------
- FIXME: Add docs.
-
- """
- lineage_dict = {}
- if len(matrix):
-
- daughter = matrix[:, n_keys]
- mother_global_id = matrix[:, :n_keys]
-
- iterator = groupby(
- zip(mother_global_id, daughter), lambda x: str(x[0])
- )
- lineage_dict = {key: [x[1] for x in group] for key, group in iterator}
-
- def str_to_tuple(k: str) -> t.Tuple[int, ...]:
- return tuple([int(x) for x in re.findall("[0-9]+", k)])
-
- # Convert keys from str to tuple
- lineage_dict = {
- str_to_tuple(k): sorted(v) for k, v in lineage_dict.items()
- }
-
- return lineage_dict
diff --git a/src/aliby/baby_client.py b/src/aliby/baby_client.py
index f13c5b027924c0f05a733edbc58a7d2f7ecb7280..3f0c174560d4763e41b0c379eb514887c8a2d379 100644
--- a/src/aliby/baby_client.py
+++ b/src/aliby/baby_client.py
@@ -13,7 +13,7 @@ import baby.errors
import h5py
import numpy as np
import requests
-from agora.abc import ParametersABC
+from agora.abc import ParametersABC, StepABC
from baby import modelsets
from baby.brain import BabyBrain
from baby.crawler import BabyCrawler
@@ -128,7 +128,7 @@ class BabyParameters(ParametersABC):
self.update("model_config", weights_flattener)
-class BabyRunner:
+class BabyRunner(StepABC):
"""A BabyRunner object for cell segmentation.
Does segmentation one time point at a time."""
@@ -157,92 +157,20 @@ class BabyRunner:
.swapaxes(1, 2)
)
- def run_tp(self, tp, with_edgemasks=True, assign_mothers=True, **kwargs):
+ def _run_tp(self, tp, with_edgemasks=True, assign_mothers=True, **kwargs):
"""Simulating processing time with sleep"""
# Access the image
t = perf_counter()
img = self.get_data(tp)
- logging.debug(f"Timing:BF_fetch:{perf_counter()-t}s")
- t = perf_counter()
segmentation = self.crawler.step(
img,
with_edgemasks=with_edgemasks,
assign_mothers=assign_mothers,
**kwargs,
)
- logging.debug(f"Timing:crawler_step:{perf_counter()-t}s")
return format_segmentation(segmentation, tp)
-class BabyClient:
- """A dummy BabyClient object for Dask Demo.
-
-
- Does segmentation one time point at a time.
- Should work better with the parallelisation.
- """
-
- bf_channel = 0
- model_name = "prime95b_brightfield_60x_5z"
- url = "http://localhost:5101"
- max_tries = 50
- sleep_time = 0.1
-
- def __init__(self, tiler):
- self.tiler = tiler
- self._session = None
-
- @property
- def session(self):
- if self._session is None:
- r_session = requests.get(self.url + f"/session/{self.model_name}")
- r_session.raise_for_status()
- self._session = r_session.json()["sessionid"]
- return self._session
-
- def get_data(self, tp):
- return self.tiler.get_tp_data(tp, self.bf_channel).swapaxes(1, 3)
-
- # def queue_image(self, img, **kwargs):
- # bit_depth = img.dtype.itemsize * 8 # bit depth = byte_size * 8
- # data = create_request(img.shape, bit_depth, img, **kwargs)
- # status = requests.post(
- # self.url + f"/segment?sessionid={self.session}",
- # data=data,
- # headers={"Content-Type": data.content_type},
- # )
- # status.raise_for_status()
- # return status
-
- def get_segmentation(self):
- try:
- seg_response = requests.get(
- self.url + f"/segment?sessionid={self.session}", timeout=120
- )
- seg_response.raise_for_status()
- result = seg_response.json()
- except Timeout as e:
- raise e
- except HTTPError as e:
- raise e
- return result
-
- def run_tp(self, tp, **kwargs):
- # Get data
- img = self.get_data(tp)
- # Queue image
- _ = self.queue_image(img, **kwargs)
- # Get segmentation
- for _ in range(self.max_tries):
- try:
- seg = self.get_segmentation()
- break
- except (Timeout, HTTPError):
- time.sleep(self.sleep_time)
- continue
- return format_segmentation(seg, tp)
-
-
def choose_model_from_params(
modelset_filter=None,
camera="prime95b",
diff --git a/src/aliby/bin/__init__.py b/src/aliby/bin/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..a051bb08ab05961eccb8250a04ca9bf4f0e382d9
--- /dev/null
+++ b/src/aliby/bin/__init__.py
@@ -0,0 +1 @@
+#!/usr/bin/env jupyter
diff --git a/src/aliby/bin/run.py b/src/aliby/bin/run.py
new file mode 100644
index 0000000000000000000000000000000000000000..e8495580d456783a35135a89641b15321c4fb02e
--- /dev/null
+++ b/src/aliby/bin/run.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env jupyter
+
+
+def run():
+ import argparse
+
+ from aliby.pipeline import Pipeline, PipelineParameters
+
+ parser = argparse.ArgumentParser(
+ prog="aliby-run",
+ description="Run a default microscopy analysis pipeline",
+ )
+
+ param_values = {
+ "expt_id": None,
+ "distributed": 2,
+ "tps": 2,
+ "directory": "./data",
+ "filter": 0,
+ "host": None,
+ "username": None,
+ "password": None,
+ }
+
+ def _cast_str(x: str or None):
+ """
+ Cast string as int if possible. If Nonetype return None.
+ """
+ if x:
+ try:
+ return int(x)
+ except:
+ return x
+
+ for k in param_values:
+ parser.add_argument(f"--{k}", action="store")
+
+ args = parser.parse_args()
+
+ for k in param_values:
+ if passed_value := _cast_str(getattr(args, k)):
+
+ param_values[k] = passed_value
+
+ params = PipelineParameters.default(general=param_values)
+ p = Pipeline(params)
+
+ p.run()
diff --git a/src/aliby/haystack.py b/src/aliby/haystack.py
index 340584c5bf95e7d57c9417fc4731c80fb6fe02be..d1368ffd7715a4c6285cf6a6b2e2e79edf5d6319 100644
--- a/src/aliby/haystack.py
+++ b/src/aliby/haystack.py
@@ -37,8 +37,6 @@ def timer(func, *args, **kwargs):
################## CUSTOM OBJECTS ##################################
-
-
class ModelPredictor:
"""Generic object that takes a NN and returns the prediction.
@@ -77,32 +75,3 @@ class ModelPredictorWriter(DynamicWriter):
"timepoint": ((None,), np.uint16),
}
self.group = f"{self.name}_info"
-
-
-class Saver:
- channel_names = {0: "BrightField", 1: "GFP"}
-
- def __init__(self, tiler, save_directory, pos_name):
- """This class straight up saves the trap data for use with neural networks in the future."""
- self.tiler = tiler
- self.name = pos_name
- self.save_dir = Path(save_directory)
-
- def channel_dir(self, index):
- ch_dir = self.save_dir / self.channel_names[index]
- if not ch_dir.exists():
- ch_dir.mkdir()
- return ch_dir
-
- def get_data(self, tp, ch):
- return self.tiler.get_tp_data(tp, ch).swapaxes(1, 3).swapaxes(1, 2)
-
- def cache(self, tp):
- # Get a given time point
- # split into channels
- for ch in self.channel_names:
- ch_dir = self.channel_dir(ch)
- data = self.get_data(tp, ch)
- for tid, trap in enumerate(data):
- np.save(ch_dir / f"{self.name}_{tid}_{tp}.npy", trap)
- return
diff --git a/src/aliby/io/dataset.py b/src/aliby/io/dataset.py
index b09567ed933cdafcf4b411d5a1c2b6d4f217318d..927a35d4c2b5ab5d0b482d7f07ebe5c2aa85bad5 100644
--- a/src/aliby/io/dataset.py
+++ b/src/aliby/io/dataset.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python3
"""
Dataset is a group of classes to manage multiple types of experiments:
- - Remote experiments on an OMERO server
+ - Remote experiments on an OMERO server (located in src/aliby/io/omero.py)
- Local experiments in a multidimensional OME-TIFF image containing the metadata
- Local experiments in a directory containing multiple positions in independent images with or without metadata
"""
@@ -11,7 +11,6 @@ import time
import typing as t
from abc import ABC, abstractproperty, abstractmethod
from pathlib import Path, PosixPath
-from typing import Union
try:
import omero
@@ -19,7 +18,6 @@ except ModuleNotFoundError:
print("Warning: Cannot import omero.")
from agora.io.bridge import BridgeH5
from aliby.io.image import ImageLocalOME
-from aliby.io.omero import BridgeOmero
def dispatch_dataset(expt_id: int or str, **kwargs):
@@ -35,7 +33,11 @@ def dispatch_dataset(expt_id: int or str, **kwargs):
Callable Dataset instance, either network-dependent or local.
"""
if isinstance(expt_id, int): # Is an experiment online
+
+ from aliby.io.omero import Dataset
+
return Dataset(expt_id, **kwargs)
+
elif isinstance(expt_id, str): # Files or Dir
expt_path = Path(expt_id)
if expt_path.is_dir():
@@ -54,7 +56,7 @@ class DatasetLocalABC(ABC):
_valid_suffixes = ("tiff", "png")
_valid_meta_suffixes = ("txt", "log")
- def __init__(self, dpath: Union[str, PosixPath], *args, **kwargs):
+ def __init__(self, dpath: t.Union[str, PosixPath], *args, **kwargs):
self.path = Path(dpath)
def __enter__(self):
@@ -110,7 +112,7 @@ class DatasetLocalDir(DatasetLocalABC):
It relies on ImageDir to manage images.
"""
- def __init__(self, dpath: Union[str, PosixPath], *args, **kwargs):
+ def __init__(self, dpath: t.Union[str, PosixPath], *args, **kwargs):
super().__init__(dpath)
@property
@@ -141,7 +143,7 @@ class DatasetLocalOME(DatasetLocalABC):
It uses the standard OME-TIFF file format.
"""
- def __init__(self, dpath: Union[str, PosixPath], *args, **kwargs):
+ def __init__(self, dpath: t.Union[str, PosixPath], *args, **kwargs):
super().__init__(dpath)
assert len(self.get_images()), "No .tiff files found"
@@ -157,102 +159,3 @@ class DatasetLocalOME(DatasetLocalABC):
for suffix in self._valid_suffixes
for f in self.path.glob(f"*.{suffix}")
}
-
-
-class Dataset(BridgeOmero):
- def __init__(self, expt_id, **server_info):
- self.ome_id = expt_id
-
- super().__init__(**server_info)
-
- @property
- def name(self):
- return self.ome_class.getName()
-
- @property
- def date(self):
- return self.ome_class.getDate()
-
- @property
- def unique_name(self):
- return "_".join(
- (
- str(self.ome_id),
- self.date.strftime("%Y_%m_%d").replace("/", "_"),
- self.name,
- )
- )
-
- def get_images(self):
- return {
- im.getName(): im.getId() for im in self.ome_class.listChildren()
- }
-
- @property
- def files(self):
- if not hasattr(self, "_files"):
- self._files = {
- x.getFileName(): x
- for x in self.ome_class.listAnnotations()
- if isinstance(x, omero.gateway.FileAnnotationWrapper)
- }
- if not len(self._files):
- raise Exception(
- "exception:metadata: experiment has no annotation files."
- )
- elif len(self.file_annotations) != len(self._files):
- raise Exception("Number of files and annotations do not match")
-
- return self._files
-
- @property
- def tags(self):
- if self._tags is None:
- self._tags = {
- x.getname(): x
- for x in self.ome_class.listAnnotations()
- if isinstance(x, omero.gateway.TagAnnotationWrapper)
- }
- return self._tags
-
- def cache_logs(self, root_dir):
- valid_suffixes = ("txt", "log")
- for name, annotation in self.files.items():
- filepath = root_dir / annotation.getFileName().replace("/", "_")
- if (
- any([str(filepath).endswith(suff) for suff in valid_suffixes])
- and not filepath.exists()
- ):
- # save only the text files
- with open(str(filepath), "wb") as fd:
- for chunk in annotation.getFileInChunks():
- fd.write(chunk)
- return True
-
- @classmethod
- def from_h5(
- cls,
- filepath: t.Union[str, PosixPath],
- ):
- """Instatiate Dataset from a hdf5 file.
-
- Parameters
- ----------
- cls : Image
- Image class
- filepath : t.Union[str, PosixPath]
- Location of hdf5 file.
-
- Examples
- --------
- FIXME: Add docs.
-
- """
- # metadata = load_attributes(filepath)
- bridge = BridgeH5(filepath)
- dataset_keys = ("omero_id", "omero_id,", "dataset_id")
- for k in dataset_keys:
- if k in bridge.meta_h5:
- return cls(
- bridge.meta_h5[k], **cls.server_info_from_h5(filepath)
- )
diff --git a/src/aliby/io/image.py b/src/aliby/io/image.py
index fe6c0759f9d746fe91ad78d93a19db8f23a9e0f3..7e021fcdfa6e582443161addcb5d5eb9e4d8bc1d 100644
--- a/src/aliby/io/image.py
+++ b/src/aliby/io/image.py
@@ -1,50 +1,35 @@
#!/usr/bin/env python3
+"""
+Image: Loads images and registers them.
+
+Image instances loads images from a specified directory into an object that
+also contains image properties such as name and metadata. Pixels from images
+are stored in dask arrays; the standard way is to store them in 5-dimensional
+arrays: T(ime point), C(channel), Z(-stack), Y, X.
+
+This module consists of a base Image class (BaseLocalImage). ImageLocalOME
+handles local OMERO images. ImageDir handles cases in which images are split
+into directories, with each time point and channel having its own image file.
+ImageDummy is a dummy class for silent failure testing.
+"""
import typing as t
-from abc import ABC, abstractproperty
+from abc import ABC, abstractmethod, abstractproperty
from datetime import datetime
+from importlib_resources import files
from pathlib import Path, PosixPath
import dask.array as da
-import numpy as np
import xmltodict
-from dask import delayed
from dask.array.image import imread
-
-try:
- from omero.model import enums as omero_enums
-except ModuleNotFoundError:
- print("Warning: Cannot import omero_enums.")
from tifffile import TiffFile
-from yaml import safe_load
-from agora.io.bridge import BridgeH5
from agora.io.metadata import dir_to_meta
-try:
- from aliby.io.omero import BridgeOmero
-except ModuleNotFoundError:
-
- class BridgeOmero:
- # dummy class necessary when the omero module is not installed
- pass
-
- print("Warning: Cannot import BridgeOmero.")
-# convert OMERO definitions into numpy types
-try:
- PIXEL_TYPES = {
- omero_enums.PixelsTypeint8: np.int8,
- omero_enums.PixelsTypeuint8: np.uint8,
- omero_enums.PixelsTypeint16: np.int16,
- omero_enums.PixelsTypeuint16: np.uint16,
- omero_enums.PixelsTypeint32: np.int32,
- omero_enums.PixelsTypeuint32: np.uint32,
- omero_enums.PixelsTypefloat: np.float32,
- omero_enums.PixelsTypedouble: np.float64,
- }
-except NameError:
- print("Warning: PIXEL_TYPES are not defined.")
+def get_examples_dir():
+ """Get examples directory which stores dummy image for tiler"""
+ return files("aliby").parent.parent / "examples" / "tiler"
def get_image_class(source: t.Union[str, int, t.Dict[str, str], PosixPath]):
@@ -52,6 +37,8 @@ def get_image_class(source: t.Union[str, int, t.Dict[str, str], PosixPath]):
Wrapper to pick the appropiate Image class depending on the source of data.
"""
if isinstance(source, int):
+ from aliby.io.omero import Image
+
instatiator = Image
elif isinstance(source, dict) or (
isinstance(source, (str, PosixPath)) and Path(source).is_dir()
@@ -67,10 +54,10 @@ def get_image_class(source: t.Union[str, int, t.Dict[str, str], PosixPath]):
class BaseLocalImage(ABC):
"""
- Base class to set path and provide context management method.
+ Base Image class to set path and provide context management method.
"""
- _default_dimorder = "tczxy"
+ _default_dimorder = "tczyx"
def __init__(self, path: t.Union[str, PosixPath]):
# If directory, assume contents are naturally sorted
@@ -79,6 +66,12 @@ class BaseLocalImage(ABC):
def __enter__(self):
return self
+ def __exit__(self, *exc):
+ for e in exc:
+ if e is not None:
+ print(e)
+ return False
+
def rechunk_data(self, img):
# Format image using x and y size from metadata.
@@ -88,12 +81,16 @@ class BaseLocalImage(ABC):
1,
1,
1,
- self._meta["size_x"],
self._meta["size_y"],
+ self._meta["size_x"],
),
)
return self._rechunked_img
+ @abstractmethod
+ def get_data_lazy(self) -> da.Array:
+ pass
+
@abstractproperty
def name(self):
pass
@@ -106,23 +103,130 @@ class BaseLocalImage(ABC):
def data(self):
return self.get_data_lazy()
- def __enter__(self):
- return self
-
- def __exit__(self, *exc):
- for e in exc:
- if e is not None:
- print(e)
- return False
-
@property
def metadata(self):
return self._meta
+class ImageDummy(BaseLocalImage):
+ """
+ Dummy Image class.
+
+ ImageDummy mimics the other Image classes in such a way that it is accepted
+ by Tiler. The purpose of this class is for testing, in particular,
+ identifying silent failures. If something goes wrong, we should be able to
+ know whether it is because of bad parameters or bad input data.
+
+ For the purposes of testing parameters, ImageDummy assumes that we already
+ know the tiler parameters before Image instances are instantiated. This is
+ true for a typical pipeline run.
+ """
+
+ def __init__(self, tiler_parameters: dict):
+ """Builds image instance
+
+ Parameters
+ ----------
+ tiler_parameters : dict
+ Tiler parameters, in dict form. Following
+ aliby.tile.tiler.TilerParameters, the keys are: "tile_size" (size of
+ tile), "ref_channel" (reference channel for tiling), and "ref_z"
+ (reference z-stack, 0 to choose a default).
+ """
+ self.ref_channel = tiler_parameters["ref_channel"]
+ self.ref_z = tiler_parameters["ref_z"]
+
+ # Goal: make Tiler happy.
+ @staticmethod
+ def pad_array(
+ image_array: da.Array,
+ dim: int,
+ n_empty_slices: int,
+ image_position: int = 0,
+ ):
+ """Extends a dimension in a dask array and pads with zeros
+
+ Extends a dimension in a dask array that has existing content, then pads
+ with zeros.
+
+ Parameters
+ ----------
+ image_array : da.Array
+ Input dask array
+ dim : int
+ Dimension in which to extend the dask array.
+ n_empty_slices : int
+ Number of empty slices to extend the dask array by, in the specified
+ dimension/axis.
+ image_position : int
+ Position within the new dimension to place the input arary, default 0
+ (the beginning).
+
+ Examples
+ --------
+ ```
+ extended_array = pad_array(
+ my_da_array, dim = 2, n_empty_slices = 4, image_position = 1)
+ ```
+ Extends a dask array called `my_da_array` in the 3rd dimension
+ (dimensions start from 0) by 4 slices, filled with zeros. And puts the
+ original content in slice 1 of the 3rd dimension
+ """
+ # Concats zero arrays with same dimensions as image_array, and puts
+ # image_array as first element in list of arrays to be concatenated
+ zeros_array = da.zeros_like(image_array)
+ return da.concatenate(
+ [
+ *([zeros_array] * image_position),
+ image_array,
+ *([zeros_array] * (n_empty_slices - image_position)),
+ ],
+ axis=dim,
+ )
+
+ # Logic: We want to return a image instance
+ def get_data_lazy(self) -> da.Array:
+ """Return 5D dask array. For lazy-loading multidimensional tiff files. Dummy image."""
+ examples_dir = get_examples_dir()
+ # TODO: Make this robust to having multiple TIFF images, one for each z-section,
+ # all falling under the same "pypipeline_unit_test_00_000001_Brightfield_*.tif"
+ # naming scheme. The aim is to create a multidimensional dask array that stores
+ # the z-stacks.
+ img_filename = "pypipeline_unit_test_00_000001_Brightfield_003.tif"
+ img_path = examples_dir / img_filename
+ # img is a dask array has three dimensions: z, x, y
+ # TODO: Write a test to confirm this: If everything worked well,
+ # z = 1, x = 1200, y = 1200
+ img = imread(str(img_path))
+ # Adds t & c dimensions
+ img = da.reshape(
+ img, (1, 1, img.shape[-3], img.shape[-2], img.shape[-1])
+ )
+ # Pads t, c, and z dimensions
+ img = self.pad_array(
+ img, dim=0, n_empty_slices=199
+ ) # 200 timepoints total
+ img = self.pad_array(img, dim=1, n_empty_slices=2) # 3 channels
+ img = self.pad_array(
+ img, dim=2, n_empty_slices=4, image_position=self.ref_z
+ ) # 5 z-stacks
+ return img
+
+ @property
+ def name(self):
+ pass
+
+ @property
+ def dimorder(self):
+ pass
+
+
class ImageLocalOME(BaseLocalImage):
"""
- Fetch image from OMEXML data format, in which a multidimensional tiff image contains the metadata.
+ Local OMERO Image class.
+
+ This is a derivative Image class. It fetches an image from OMEXML data format,
+ in which a multidimensional tiff image contains the metadata.
"""
def __init__(self, path: str, dimorder=None):
@@ -222,7 +326,7 @@ class ImageDir(BaseLocalImage):
"""
Image class for the case in which all images are split in one or
multiple folders with time-points and channels as independent files.
- It inherits from Imagelocal so we only override methods that are critical.
+ It inherits from BaseLocalImage so we only override methods that are critical.
Assumptions:
- One folders per position.
@@ -278,139 +382,3 @@ class ImageDir(BaseLocalImage):
return [
k.split("_")[-1] for k in self._meta.keys() if k.startswith("size")
]
-
-
-class Image(BridgeOmero):
- """
- Load images from OMERO and gives access to the data and metadata.
- """
-
- def __init__(self, image_id: int, **server_info):
- """
- Establish the connection to the OMERO server via the Argo
- base class.
-
- Parameters
- ----------
- image_id: integer
- server_info: dictionary
- Specifies the host, username, and password as strings
- """
- self.ome_id = image_id
- super().__init__(**server_info)
-
- def init_interface(self, ome_id: int):
- self.set_id(ome_id)
- self.ome_class = self.conn.getObject("Image", ome_id)
-
- @classmethod
- def from_h5(
- cls,
- filepath: t.Union[str, PosixPath],
- ):
- """
- Instantiate Image from a hdf5 file.
-
- Parameters
- ----------
- cls : Image
- Image class
- filepath : t.Union[str, PosixPath]
- Location of hdf5 file.
-
- Examples
- --------
- FIXME: Add docs.
-
- """
- # metadata = load_attributes(filepath)
- bridge = BridgeH5(filepath)
- image_id = bridge.meta_h5["image_id"]
- # server_info = safe_load(bridge.meta_h5["parameters"])["general"][
- # "server_info"
- # ]
- return cls(image_id, **cls.server_info_from_h5(filepath))
-
- @property
- def name(self):
- return self.ome_class.getName()
-
- @property
- def data(self):
- return get_data_lazy(self.ome_class)
-
- @property
- def metadata(self):
- """
- Store metadata saved in OMERO: image size, number of time points,
- labels of channels, and image name.
- """
- meta = dict()
- meta["size_x"] = self.ome_class.getSizeX()
- meta["size_y"] = self.ome_class.getSizeY()
- meta["size_z"] = self.ome_class.getSizeZ()
- meta["size_c"] = self.ome_class.getSizeC()
- meta["size_t"] = self.ome_class.getSizeT()
- meta["channels"] = self.ome_class.getChannelLabels()
- meta["name"] = self.ome_class.getName()
- return meta
-
-
-class UnsafeImage(Image):
- """
- Loads images from OMERO and gives access to the data and metadata.
- This class is a temporary solution while we find a way to use
- context managers inside napari. It risks resulting in zombie connections
- and producing freezes in an OMERO server.
-
- """
-
- def __init__(self, image_id, **server_info):
- """
- Establishes the connection to the OMERO server via the Argo
- base class.
-
- Parameters
- ----------
- image_id: integer
- server_info: dictionary
- Specifies the host, username, and password as strings
- """
- super().__init__(image_id, **server_info)
- self.create_gate()
- self.init_wrapper()
-
- @property
- def data(self):
- try:
- return get_data_lazy(self.ome_class)
- except Exception as e:
- print(f"ERROR: Failed fetching image from server: {e}")
- self.conn.connect(False)
-
-
-def get_data_lazy(image) -> da.Array:
- """
- Get 5D dask array, with delayed reading from OMERO image.
- """
- nt, nc, nz, ny, nx = [getattr(image, f"getSize{x}")() for x in "TCZYX"]
- pixels = image.getPrimaryPixels()
- dtype = PIXEL_TYPES.get(pixels.getPixelsType().value, None)
- # using dask
- get_plane = delayed(lambda idx: pixels.getPlane(*idx))
-
- def get_lazy_plane(zct):
- return da.from_delayed(get_plane(zct), shape=(ny, nx), dtype=dtype)
-
- # 5D stack: TCZXY
- t_stacks = []
- for t in range(nt):
- c_stacks = []
- for c in range(nc):
- z_stack = []
- for z in range(nz):
- z_stack.append(get_lazy_plane((z, c, t)))
- c_stacks.append(da.stack(z_stack))
- t_stacks.append(da.stack(c_stacks))
-
- return da.stack(t_stacks)
diff --git a/src/aliby/io/omero.py b/src/aliby/io/omero.py
index 9ae5ffba80b42accfdc0be12c623e2926763bad7..62475a0ccce8e7e4f123531f92e93ee0e89806e7 100644
--- a/src/aliby/io/omero.py
+++ b/src/aliby/io/omero.py
@@ -1,15 +1,33 @@
+"""
+Tools to manage I/O using a remote OMERO server.
+"""
+
+import re
import typing as t
from abc import abstractmethod
from pathlib import PosixPath
-import re
+
+import dask.array as da
+import numpy as np
+import omero
+from dask import delayed
+from omero.gateway import BlitzGateway
+from omero.model import enums as omero_enums
+from yaml import safe_load
from agora.io.bridge import BridgeH5
-try:
- from omero.gateway import BlitzGateway
-except ModuleNotFoundError:
- print("Warning: Cannot import BlitzGateway.")
-from yaml import safe_load
+# convert OMERO definitions into numpy types
+PIXEL_TYPES = {
+ omero_enums.PixelsTypeint8: np.int8,
+ omero_enums.PixelsTypeuint8: np.uint8,
+ omero_enums.PixelsTypeint16: np.int16,
+ omero_enums.PixelsTypeuint16: np.uint16,
+ omero_enums.PixelsTypeint32: np.int32,
+ omero_enums.PixelsTypeuint32: np.uint32,
+ omero_enums.PixelsTypefloat: np.float32,
+ omero_enums.PixelsTypedouble: np.float64,
+}
class BridgeOmero:
@@ -21,9 +39,10 @@ class BridgeOmero:
def __init__(
self,
- host="islay.bio.ed.ac.uk",
- username="upload",
- password="***REMOVED***",
+ host: str = None,
+ username: str = None,
+ password: str = None,
+ ome_id: int = None,
):
"""
Parameters
@@ -32,22 +51,34 @@ class BridgeOmero:
web address of OMERO host
username: string
password : string
+ ome_id: Optional int
+ Unique identifier on Omero database. Used to fetch specific objects.
"""
+ # assert all((host, username, password)), str(f"Invalid credentials host:{host}, user:{username}, pass:{pass}")
+ assert all(
+ (host, username, password)
+ ), f"Invalid credentials. host: {host}, user: {username}, pwd: {password}"
+
self.conn = None
self.host = host
self.username = username
self.password = password
+ self.ome_id = ome_id
# standard method required for Python's with statement
def __enter__(self):
self.create_gate()
- self.init_wrapper()
return self
- def init_wrapper(self):
+ @property
+ def ome_class(self):
# Initialise Omero Object Wrapper for instances when applicable.
- if hasattr(self, "ome_id"):
+ if not hasattr(self, "_ome_class"):
+ assert (
+ self.conn.isConnected() and self.ome_id is not None
+ ), "No Blitz connection or valid omero id"
+
ome_type = [
valid_name
for valid_name in ("Dataset", "Image")
@@ -57,7 +88,11 @@ class BridgeOmero:
re.IGNORECASE,
)
][0]
- self.ome_class = self.conn.getObject(ome_type, self.ome_id)
+ self._ome_class = self.conn.getObject(ome_type, self.ome_id)
+
+ assert self._ome_class, f"{ome_type} {self.ome_id} not found."
+
+ return self._ome_class
def create_gate(self) -> bool:
self.conn = BlitzGateway(
@@ -106,10 +141,6 @@ class BridgeOmero:
def set_id(self, ome_id: int):
self.ome_id = ome_id
- @abstractmethod
- def init_interface(self):
- ...
-
@property
def file_annotations(self):
valid_annotations = [
@@ -137,3 +168,226 @@ class BridgeOmero:
**kwargs,
)
self.ome_class.linkAnnotation(file_annotation)
+
+
+class Dataset(BridgeOmero):
+ def __init__(self, expt_id: str or int, **server_info):
+ super().__init__(ome_id=expt_id, **server_info)
+
+ @property
+ def name(self):
+ return self.ome_class.getName()
+
+ @property
+ def date(self):
+ return self.ome_class.getDate()
+
+ @property
+ def unique_name(self):
+ return "_".join(
+ (
+ str(self.ome_id),
+ self.date.strftime("%Y_%m_%d").replace("/", "_"),
+ self.name,
+ )
+ )
+
+ def get_images(self):
+ return {
+ im.getName(): im.getId() for im in self.ome_class.listChildren()
+ }
+
+ @property
+ def files(self):
+ if not hasattr(self, "_files"):
+ self._files = {
+ x.getFileName(): x
+ for x in self.ome_class.listAnnotations()
+ if isinstance(x, omero.gateway.FileAnnotationWrapper)
+ }
+ if not len(self._files):
+ raise Exception(
+ "exception:metadata: experiment has no annotation files."
+ )
+ elif len(self.file_annotations) != len(self._files):
+ raise Exception("Number of files and annotations do not match")
+
+ return self._files
+
+ @property
+ def tags(self):
+ if self._tags is None:
+ self._tags = {
+ x.getname(): x
+ for x in self.ome_class.listAnnotations()
+ if isinstance(x, omero.gateway.TagAnnotationWrapper)
+ }
+ return self._tags
+
+ def cache_logs(self, root_dir):
+ valid_suffixes = ("txt", "log")
+ for _, annotation in self.files.items():
+ filepath = root_dir / annotation.getFileName().replace("/", "_")
+ if (
+ any([str(filepath).endswith(suff) for suff in valid_suffixes])
+ and not filepath.exists()
+ ):
+ # save only the text files
+ with open(str(filepath), "wb") as fd:
+ for chunk in annotation.getFileInChunks():
+ fd.write(chunk)
+ return True
+
+ @classmethod
+ def from_h5(
+ cls,
+ filepath: t.Union[str, PosixPath],
+ ):
+ """Instatiate Dataset from a hdf5 file.
+
+ Parameters
+ ----------
+ cls : Image
+ Image class
+ filepath : t.Union[str, PosixPath]
+ Location of hdf5 file.
+
+ Examples
+ --------
+ FIXME: Add docs.
+
+ """
+ # metadata = load_attributes(filepath)
+ bridge = BridgeH5(filepath)
+ dataset_keys = ("omero_id", "omero_id,", "dataset_id")
+ for k in dataset_keys:
+ if k in bridge.meta_h5:
+ return cls(
+ bridge.meta_h5[k], **cls.server_info_from_h5(filepath)
+ )
+
+
+class Image(BridgeOmero):
+ """
+ Loads images from OMERO and gives access to the data and metadata.
+ """
+
+ def __init__(self, image_id: int, **server_info):
+ """
+ Establishes the connection to the OMERO server via the Argo
+ base class.
+
+ Parameters
+ ----------
+ image_id: integer
+ server_info: dictionary
+ Specifies the host, username, and password as strings
+ """
+ super().__init__(ome_id=image_id, **server_info)
+
+ @classmethod
+ def from_h5(
+ cls,
+ filepath: t.Union[str, PosixPath],
+ ):
+ """Instatiate Image from a hdf5 file.
+
+ Parameters
+ ----------
+ cls : Image
+ Image class
+ filepath : t.Union[str, PosixPath]
+ Location of hdf5 file.
+
+ Examples
+ --------
+ FIXME: Add docs.
+
+ """
+ # metadata = load_attributes(filepath)
+ bridge = BridgeH5(filepath)
+ image_id = bridge.meta_h5["image_id"]
+ return cls(image_id, **cls.server_info_from_h5(filepath))
+
+ @property
+ def name(self):
+ return self.ome_class.getName()
+
+ @property
+ def data(self):
+ return get_data_lazy(self.ome_class)
+
+ @property
+ def metadata(self):
+ """
+ Store metadata saved in OMERO: image size, number of time points,
+ labels of channels, and image name.
+ """
+ meta = dict()
+ meta["size_x"] = self.ome_class.getSizeX()
+ meta["size_y"] = self.ome_class.getSizeY()
+ meta["size_z"] = self.ome_class.getSizeZ()
+ meta["size_c"] = self.ome_class.getSizeC()
+ meta["size_t"] = self.ome_class.getSizeT()
+ meta["channels"] = self.ome_class.getChannelLabels()
+ meta["name"] = self.ome_class.getName()
+ return meta
+
+
+class UnsafeImage(Image):
+ """
+ Loads images from OMERO and gives access to the data and metadata.
+ This class is a temporary solution while we find a way to use
+ context managers inside napari. It risks resulting in zombie connections
+ and producing freezes in an OMERO server.
+
+ """
+
+ def __init__(self, image_id, **server_info):
+ """
+ Establishes the connection to the OMERO server via the Argo
+ base class.
+
+ Parameters
+ ----------
+ image_id: integer
+ server_info: dictionary
+ Specifies the host, username, and password as strings
+ """
+ super().__init__(image_id, **server_info)
+ self.create_gate()
+
+ @property
+ def data(self):
+ try:
+ return get_data_lazy(self.ome_class)
+ except Exception as e:
+ print(f"ERROR: Failed fetching image from server: {e}")
+ self.conn.connect(False)
+
+
+def get_data_lazy(image) -> da.Array:
+ """
+ Get 5D dask array, with delayed reading from OMERO image.
+ """
+ nt, nc, nz, ny, nx = [getattr(image, f"getSize{x}")() for x in "TCZYX"]
+ pixels = image.getPrimaryPixels()
+ dtype = PIXEL_TYPES.get(pixels.getPixelsType().value, None)
+ # using dask
+ get_plane = delayed(lambda idx: pixels.getPlane(*idx))
+
+ def get_lazy_plane(zct):
+ return da.from_delayed(get_plane(zct), shape=(ny, nx), dtype=dtype)
+
+ # 5D stack: TCZXY
+ t_stacks = []
+ for t in range(nt):
+ c_stacks = []
+ for c in range(nc):
+ z_stack = []
+ for z in range(nz):
+ z_stack.append(get_lazy_plane((z, c, t)))
+ c_stacks.append(da.stack(z_stack))
+ t_stacks.append(da.stack(c_stacks))
+
+ return da.stack(t_stacks)
diff --git a/src/aliby/io/utils.py b/src/aliby/io/utils.py
deleted file mode 100644
index d59a2f65de12ece9fc0ed63ff37e7adfa2f1e801..0000000000000000000000000000000000000000
--- a/src/aliby/io/utils.py
+++ /dev/null
@@ -1,50 +0,0 @@
-import re
-import struct
-
-
-def clean_ascii(text):
- return re.sub(r"[^\x20-\x7F]", ".", text)
-
-
-def xxd(x, start=0, stop=None):
- if stop is None:
- stop = len(x)
- for i in range(start, stop, 8):
- # Row number
- print("%04d" % i, end=" ")
- # Hexadecimal bytes
- for r in range(i, i + 8):
- print("%02x" % x[r], end="")
- if (r + 1) % 4 == 0:
- print(" ", end="")
- # ASCII
- print(
- " ",
- clean_ascii(x[i : i + 8].decode("utf-8", errors="ignore")),
- " ",
- end="",
- )
- # Int32
- print(
- "{:>10} {:>10}".format(*struct.unpack("II", x[i : i + 8])),
- end=" ",
- )
- print("") # Newline
- return
-
-
-# Buffer reading functions
-def read_int(buffer, n=1):
- res = struct.unpack("I" * n, buffer.read(4 * n))
- if n == 1:
- res = res[0]
- return res
-
-
-def read_string(buffer):
- return "".join([x.decode() for x in iter(lambda: buffer.read(1), b"\x00")])
-
-
-def read_delim(buffer, n):
- delim = read_int(buffer, n)
- assert all([x == 0 for x in delim]), "Unknown nonzero value in delimiter"
diff --git a/src/aliby/lineage/__init__.py b/src/aliby/lineage/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..a051bb08ab05961eccb8250a04ca9bf4f0e382d9
--- /dev/null
+++ b/src/aliby/lineage/__init__.py
@@ -0,0 +1 @@
+#!/usr/bin/env jupyter
diff --git a/src/aliby/lineage/bud_tracker.py b/src/aliby/lineage/bud_tracker.py
new file mode 100644
index 0000000000000000000000000000000000000000..ef6a2e46820e15908dc92dc31a3eec9d398b1a14
--- /dev/null
+++ b/src/aliby/lineage/bud_tracker.py
@@ -0,0 +1,191 @@
+#!/usr/bin/env jupyter
+
+import pickle
+from os.path import join
+
+import numpy as np
+from skimage.draw import polygon
+
+from agora.track_abc import FeatureCalculator
+
+models_path = join(dirname(__file__), "./models")
+
+
+class BudTracker(FeatureCalculator):
+ def __init__(self, model=None, feats2use=None, **kwargs):
+
+ if model is None:
+ model_file = join(models_path, "mb_model_20201022.pkl")
+ with open(model_file, "rb") as file_to_load:
+ model = pickle.load(file_to_load)
+ self.model = model
+
+ if feats2use is None:
+ feats2use = ["centroid", "area", "minor_axis_length"]
+ super().__init__(feats2use, **kwargs)
+
+ self.a_ind = self.outfeats.index("area")
+ self.ma_ind = self.outfeats.index("minor_axis_length")
+
+ ### Assign mother-
+ def calc_mother_bud_stats(self, p_budneck, p_bud, masks, feats=None):
+ """
+ ---
+
+ input
+
+ :p_budneck: 2d ndarray (size_x, size_y) giving the probability that a
+ pixel corresponds to a bud neck
+ :p_bud: 2d ndarray (size_x, size_y) giving the probability that a pixel
+ corresponds to a bud
+ :masks: 3d ndarray (ncells, size_x, size_y)
+ :feats: ndarray (ncells, nfeats)
+
+ NB: ASSUMES FEATS HAVE ALREADY BEEN NORMALISED!
+
+ returns
+
+ :n2darray: 2d ndarray (ncells x ncells, n_feats) specifying,
+ for each pair of cells in the masks array, the features used for
+ mother-bud pair prediction (as per 'feats2use')
+ """
+
+ if feats is None:
+ feats = self.calc_feats_from_mask(masks)
+ elif len(feats) != len(masks):
+ raise Exception("number of features must match number of masks")
+
+ ncells = len(masks)
+
+ # Entries will be NaN unless validly specified below
+ p_bud_mat = np.nan * np.ones((ncells, ncells))
+ p_budneck_mat = np.nan * np.ones((ncells, ncells))
+ budneck_ratio_mat = np.nan * np.ones((ncells, ncells))
+ size_ratio_mat = np.nan * np.ones((ncells, ncells))
+ adjacency_mat = np.nan * np.ones((ncells, ncells))
+
+ for m in range(ncells):
+ for d in range(ncells):
+ if m == d:
+ # Mother-bud pairs can only be between different cells
+ continue
+
+ p_bud_mat[m, d] = np.mean(p_bud[masks[d].astype("bool")])
+
+ a_i = self.a_ind
+ size_ratio_mat[m, d] = feats[m, a_i] / feats[d, a_i]
+
+ # Draw rectangle
+ r_points = self.get_rpoints(feats, d, m)
+ if r_points is None:
+ continue
+ rr, cc = polygon(
+ r_points[0, :], r_points[1, :], p_budneck.shape
+ )
+ if len(rr) == 0:
+ # Rectangles with zero size are not informative
+ continue
+
+ r_im = np.zeros(p_budneck.shape, dtype="bool")
+ r_im[rr, cc] = True
+
+ # Calculate the mean of bud neck probabilities greater than some threshold
+ pbn = p_budneck[r_im].flatten()
+
+ pbn = pbn[pbn > 0.2]
+ p_budneck_mat[m, d] = np.mean(pbn) if len(pbn) > 0 else 0
+
+ # Normalise number of bud-neck positive pixels by the scale of
+ # the bud (a value proportional to circumference):
+ raw_circumf_est = np.sqrt(feats[d, a_i]) / self.pixel_size
+ budneck_ratio_mat[m, d] = pbn.sum() / raw_circumf_est
+
+ # Adjacency is the proportion of the joining rectangle that overlaps the mother daughter union
+ md_union = masks[m] | masks[d]
+ adjacency_mat[m, d] = np.sum(md_union & r_im) / np.sum(r_im)
+
+ return np.hstack(
+ [
+ s.flatten()[:, np.newaxis]
+ for s in (
+ p_bud_mat,
+ size_ratio_mat,
+ p_budneck_mat,
+ budneck_ratio_mat,
+ adjacency_mat,
+ )
+ ]
+ )
+
+ def predict_mother_bud(self, p_budneck, p_bud, masks, feats=None):
+ """
+ ---
+
+ input
+
+ :p_budneck: 2d ndarray (size_x, size_y) giving the probability that a
+ pixel corresponds to a bud neck
+ :p_bud: 2d ndarray (size_x, size_y) giving the probability that a pixel
+ corresponds to a bud
+ :masks: 3d ndarray (ncells, size_x, size_y)
+ :feats: ndarray (ncells, nfeats)
+
+ returns
+
+ :n2darray: 2d ndarray (ncells, ncells) giving probability that a cell
+ (row) is a mother to another cell (column)
+ """
+
+ ncells = len(masks)
+
+ mb_stats = self.calc_mother_bud_stats(p_budneck, p_bud, masks, feats)
+
+ good_stats = ~np.isnan(mb_stats).any(axis=1)
+ # Assume probability of bud assignment for any rows that are NaN will
+ # be zero
+ ba_probs = np.zeros(ncells**2)
+ if good_stats.any():
+ ba_probs[good_stats] = self.model.predict_proba(
+ mb_stats[good_stats, :]
+ )[:, 1]
+ ba_probs = ba_probs.reshape((ncells,) * 2)
+
+ return ba_probs
+
+ def get_rpoints(self, feats, d, m):
+ """
+ Draw a rectangle in the budneck of cells
+ ---
+
+ NB: ASSUMES FEATS HAVE ALREADY BEEN NORMALISED!
+
+ input
+
+ feats: 2d ndarray (ncells, nfeats)
+
+ returns
+
+ r_points: 2d ndarray (2,4) with the coordinates of the rectangle corner
+
+ """
+
+ # Get un-scaled features for m-d pair
+ descaled_feats = feats / self.pixel_size
+ m_centre = descaled_feats[m, :2]
+ d_centre = descaled_feats[d, :2]
+ r_width = np.max((2, descaled_feats[d, self.ma_ind] * 0.25))
+
+ # Draw connecting rectangle
+ r_hvec = d_centre - m_centre
+ r_wvec = np.matmul(np.array([[0, -1], [1, 0]]), r_hvec)
+ r_wvec_len = np.linalg.norm(r_wvec)
+ if r_wvec_len == 0:
+ return None
+ r_wvec = r_width * r_wvec / r_wvec_len
+ r_points = np.zeros((2, 4))
+ r_points[:, 0] = m_centre - 0.5 * r_wvec
+ r_points[:, 1] = r_points[:, 0] + r_hvec
+ r_points[:, 2] = r_points[:, 1] + r_wvec
+ r_points[:, 3] = r_points[:, 2] - r_hvec
+
+ return r_points
diff --git a/src/aliby/pipeline.py b/src/aliby/pipeline.py
index ea56674601afe7d9c60637a761b7a14baf55f5cd..2977cada71b37a131f30e25103d6d668853878d0 100644
--- a/src/aliby/pipeline.py
+++ b/src/aliby/pipeline.py
@@ -9,8 +9,6 @@ import typing as t
from copy import copy
from importlib.metadata import version
from pathlib import Path, PosixPath
-from time import perf_counter
-from typing import Union
import h5py
import numpy as np
@@ -22,7 +20,7 @@ from agora.abc import ParametersABC, ProcessABC
from agora.io.metadata import MetaData, parse_logfiles
from agora.io.reader import StateReader
from agora.io.signal import Signal
-from agora.io.writer import ( # BabyWriter,
+from agora.io.writer import (
LinearBabyWriter,
StateWriter,
TilerWriter,
@@ -36,17 +34,20 @@ from extraction.core.extractor import Extractor, ExtractorParameters
from extraction.core.functions.defaults import exparams_from_meta
from postprocessor.core.processor import PostProcessor, PostProcessorParameters
-# from postprocessor.compiler import ExperimentCompiler, PageOrganiser
-
-logging.basicConfig(
- filename="aliby.log",
- filemode="w",
- format="%(name)s - %(levelname)s - %(message)s",
- level=logging.DEBUG,
-)
-
class PipelineParameters(ParametersABC):
+ """
+ Parameters that host what is run and how. It takes a list of dictionaries, one for
+ general in collection:
+ pass dictionary for each step
+ --------------------
+ expt_id: int or str Experiment id (if integer) or local path (if string).
+ directory: str Directory into which results are dumped. Default is "../data"
+
+ Provides default parameters for the entire pipeline. This downloads the logfiles and sets the default
+ timepoints and extraction parameters from there.
+ """
+
_pool_index = None
def __init__(
@@ -67,16 +68,7 @@ class PipelineParameters(ParametersABC):
baby={},
extraction={},
postprocessing={},
- # reporting={},
):
- """
- Load unit test experiment
- :expt_id: Experiment id
- :directory: Output directory
-
- Provides default parameters for the entire pipeline. This downloads the logfiles and sets the default
- timepoints and extraction parameters from there.
- """
expt_id = general.get("expt_id", 19993)
if isinstance(expt_id, PosixPath):
expt_id = str(expt_id)
@@ -85,7 +77,8 @@ class PipelineParameters(ParametersABC):
directory = Path(general.get("directory", "../data"))
with dispatch_dataset(
- expt_id, **general.get("server_info", {})
+ expt_id,
+ **{k: general.get(k) for k in ("host", "username", "password")},
) as conn:
directory = directory / conn.unique_name
if not directory.exists():
@@ -95,11 +88,13 @@ class PipelineParameters(ParametersABC):
try:
meta_d = MetaData(directory, None).load_logs()
except Exception as e:
+ logging.getLogger("aliby").warn(
+ f"WARNING:Metadata: error when loading: {e}"
+ )
minimal_default_meta = {
"channels": ["Brightfield"],
"ntps": [2000],
}
- print(f"WARNING:Metadata: error when loading: {e}")
# Set minimal metadata
meta_d = minimal_default_meta
@@ -118,6 +113,8 @@ class PipelineParameters(ParametersABC):
thresh_trap_area=0.9,
ntps_to_eval=5,
),
+ logfile_level="INFO",
+ use_explog=True,
)
}
@@ -187,6 +184,29 @@ class Pipeline(ProcessABC):
store = Path(store)
self.store = store
+ @staticmethod
+ def setLogger(
+ folder, file_level: str = "INFO", stream_level: str = "WARNING"
+ ):
+
+ logger = logging.getLogger("aliby")
+ logger.setLevel(getattr(logging, file_level))
+ formatter = logging.Formatter(
+ "%(asctime)s - %(levelname)s:%(message)s",
+ datefmt="%Y-%m-%dT%H:%M:%S%z",
+ )
+
+ ch = logging.StreamHandler()
+ ch.setLevel(getattr(logging, stream_level))
+ ch.setFormatter(formatter)
+ logger.addHandler(ch)
+
+ # create file handler which logs even debug messages
+ fh = logging.FileHandler(Path(folder) / "aliby.log", "w+")
+ fh.setLevel(getattr(logging, file_level))
+ fh.setFormatter(formatter)
+ logger.addHandler(fh)
+
@classmethod
def from_yaml(cls, fpath):
# This is just a convenience function, think before implementing
@@ -257,23 +277,34 @@ class Pipeline(ProcessABC):
return cls(pipeline_parameters, store=directory)
+ @property
+ def _logger(self):
+ return logging.getLogger("aliby")
+
def run(self):
- # Config holds the general information, use in main
- # Steps holds the description of tasks with their parameters
- # Steps: all holds general tasks
- # steps: strain_name holds task for a given strain
+ """
+ Config holds the general information, use in main
+ Steps: all holds general tasks
+ steps: strain_name holds task for a given strain
+ """
+
config = self.parameters.to_dict()
expt_id = config["general"]["id"]
distributed = config["general"]["distributed"]
pos_filter = config["general"]["filter"]
root_dir = Path(config["general"]["directory"])
+ self.server_info = {
+ k: config["general"].get(k)
+ for k in ("host", "username", "password")
+ }
- print("Searching OMERO")
+ dispatcher = dispatch_dataset(expt_id, **self.server_info)
+ logging.getLogger("aliby").info(
+ f"Fetching data using {dispatcher.__class__.__name__}"
+ )
# Do all all initialisations
- with dispatch_dataset(
- expt_id, **self.general.get("server_info", {})
- ) as conn:
+ with dispatcher as conn:
image_ids = conn.get_images()
directory = self.store or root_dir / conn.unique_name
@@ -288,6 +319,8 @@ class Pipeline(ProcessABC):
self.parameters.general["directory"] = str(directory)
config["general"]["directory"] = directory
+ self.setLogger(directory)
+
# Filter TODO integrate filter onto class and add regex
def filt_int(d: dict, filt: int):
return {k: v for i, (k, v) in enumerate(d.items()) if i == filt}
@@ -295,7 +328,7 @@ class Pipeline(ProcessABC):
def filt_str(image_ids: dict, filt: str):
return {k: v for k, v in image_ids.items() if re.search(filt, k)}
- def pick_filter(image_ids: dict, filt: Union[int, str]):
+ def pick_filter(image_ids: dict, filt: int or str):
if isinstance(filt, str):
image_ids = filt_str(image_ids, filt)
elif isinstance(filt, int):
@@ -332,7 +365,7 @@ class Pipeline(ProcessABC):
def create_pipeline(
self,
- image_id: t.Tuple[str, t.Union[str, PosixPath, int]],
+ image_id: t.Tuple[str, str or PosixPath or int],
index: t.Optional[int] = None,
):
""" """
@@ -370,7 +403,7 @@ class Pipeline(ProcessABC):
min_process_from = min(process_from.values())
with get_image_class(image_id)(
- image_id, **self.general.get("server_info", {})
+ image_id, **self.server_info
) as image:
# Initialise Steps
@@ -432,15 +465,10 @@ class Pipeline(ProcessABC):
for step in self.iterative_steps:
if i >= process_from[step]:
- t = perf_counter()
result = steps[step].run_tp(
i, **run_kwargs.get(step, {})
)
- logging.debug(
- f"Timing:{step}:{perf_counter() - t}s"
- )
if step in loaded_writers:
- t = perf_counter()
loaded_writers[step].write(
data=result,
overwrite=writer_ow_kwargs.get(
@@ -449,16 +477,13 @@ class Pipeline(ProcessABC):
tp=i,
meta={"last_processed": i},
)
- logging.debug(
- f"Timing:Writing-{step}:{perf_counter() - t}s"
- )
# Step-specific actions
if (
step == "tiler"
and i == min_process_from
):
- print(
+ logging.getLogger("aliby").info(
f"Found {steps['tiler'].n_traps} traps in {image.name}"
)
elif (
@@ -482,18 +507,15 @@ class Pipeline(ProcessABC):
frac_clogged_traps = self.check_earlystop(
filename, earlystop, steps["tiler"].tile_size
)
- logging.debug(
- f"Quality:Clogged_traps:{frac_clogged_traps}"
+ self._log(
+ f"{name}:Clogged_traps:{frac_clogged_traps}"
)
frac = np.round(frac_clogged_traps * 100)
pbar.set_postfix_str(f"{frac} Clogged")
else: # Stop if more than X% traps are clogged
- logging.debug(
- f"EarlyStop:{earlystop['thresh_pos_clogged']*100}% traps clogged at time point {i}"
- )
- print(
- f"Stopping analysis at time {i} with {frac_clogged_traps} clogged traps"
+ self._log(
+ f"{name}:Analysis stopped early at time {i} with {frac_clogged_traps} clogged traps"
)
meta.add_fields({"end_status": "Clogged"})
break
@@ -507,14 +529,14 @@ class Pipeline(ProcessABC):
)
PostProcessor(filename, post_proc_params).run()
+ self._log("Analysis finished successfully.", "info")
return 1
- except Exception as e: # bug during setup or runtime
+ except Exception as e: # Catch bugs during setup or runtime
logging.exception(
- f"Caught exception in worker thread (x = {name}):",
+ f"{name}: Exception caught.",
exc_info=True,
)
- print(f"Caught exception in worker thread (x = {name}):")
# This prints the type, value, and stack trace of the
# current exception being handled.
traceback.print_exc()
@@ -522,15 +544,6 @@ class Pipeline(ProcessABC):
finally:
_close_session(session)
- # try:
- # compiler = ExperimentCompiler(None, filepath)
- # tmp = compiler.run()
- # po = PageOrganiser(tmp, grid_spec=(3, 2))
- # po.plot()
- # po.save(fullpath / f"{directory}report.pdf")
- # except Exception as e:
- # print("Report failed: {}".format(e))
-
@staticmethod
def check_earlystop(filename: str, es_parameters: dict, tile_size: int):
s = Signal(filename)
@@ -643,9 +656,7 @@ class Pipeline(ProcessABC):
directory = general_config["directory"]
trackers_state: t.List[np.ndarray] = []
- with get_image_class(image_id)(
- image_id, **self.general.get("server_info", {})
- ) as image:
+ with get_image_class(image_id)(image_id, **self.server_info) as image:
filename = Path(f"{directory}/{image.name}.h5")
meta = MetaData(directory, filename)
@@ -664,6 +675,7 @@ class Pipeline(ProcessABC):
# If no previous segmentation and keep tiler
if filename.exists():
+ self._log("Result file exists.", "info")
if not ow["tiler"]:
steps["tiler"] = Tiler.from_hdf5(image, filename)
try:
@@ -685,30 +697,9 @@ class Pipeline(ProcessABC):
except Exception:
pass
- # Delete datasets to overwrite and update pipeline data
- # Use existing parameters
- # with h5py.File(filename, "a") as f:
- # pparams = PipelineParameters.from_yaml(
- # f.attrs["parameters"]
- # ).to_dict()
-
- # for k, v in ow.items():
- # if v:
- # for gname in self.writer_groups[k]:
- # if gname in f:
- # del f[gname]
-
- # pparams[k] = config[k]
- # meta.add_fields(
- # {
- # "parameters": PipelineParameters.from_dict(
- # pparams
- # ).to_yaml()
- # },
- # overwrite=True,
- # )
-
- meta.run()
+ if config["general"]["use_explog"]:
+ meta.run()
+
meta.add_fields( # Add non-logfile metadata
{
"aliby_version": version("aliby"),
diff --git a/src/aliby/tile/tiler.py b/src/aliby/tile/tiler.py
index 905cecf59a52b77a1c5d22012becb8929424e70f..34aa89d1c2af76976f134a0fd38d58919cedc9bb 100644
--- a/src/aliby/tile/tiler.py
+++ b/src/aliby/tile/tiler.py
@@ -15,7 +15,7 @@ One key method is Tiler.run.
The image-processing is performed by traps/segment_traps.
-The experiment is stored as an array with a standard indexing order of (Time, Channels, Z-stack, Y, X).
+The experiment is stored as an array with a standard indexing order of (Time, Channels, Z-stack, X, Y).
"""
import re
import typing as t
@@ -28,9 +28,9 @@ import h5py
import numpy as np
from skimage.registration import phase_cross_correlation
-from agora.abc import ParametersABC, ProcessABC
+from agora.abc import ParametersABC, StepABC
from agora.io.writer import BridgeH5
-from aliby.io.image import Image, ImageLocalOME, ImageDir
+from aliby.io.image import ImageLocalOME, ImageDir, ImageDummy
from aliby.tile.traps import segment_traps
@@ -200,7 +200,7 @@ class TilerParameters(ParametersABC):
_defaults = {"tile_size": 117, "ref_channel": "Brightfield", "ref_z": 0}
-class Tiler(ProcessABC):
+class Tiler(StepABC):
"""
Remote Timelapse Tiler.
@@ -243,12 +243,49 @@ class Tiler(ProcessABC):
for ch, zsect in zip(self.channels, metadata["zsections"])
}
except Exception as e:
- print(f"Warning:Tiler: No z_perchannel data: {e}")
+ self._log(f"No z_perchannel data: {e}")
self.tile_size = self.tile_size or min(self.image.shape[-2:])
@classmethod
- def from_image(cls, image: Image, parameters: TilerParameters):
+ def dummy(cls, parameters: dict):
+ """
+ Instantiate dummy Tiler from dummy image
+
+ If image.dimorder exists dimensions are saved in that order.
+ Otherwise default to "tczyx".
+
+ Parameters
+ ----------
+ parameters: dictionary output of an instance of TilerParameters
+ """
+ imgdmy_obj = ImageDummy(parameters)
+ dummy_image = imgdmy_obj.get_data_lazy()
+ # Default to "tczyx" if image.dimorder is None
+ dummy_omero_metadata = {
+ f"size_{dim}": dim_size
+ for dim, dim_size in zip(
+ imgdmy_obj.dimorder or "tczyx", dummy_image.shape
+ )
+ }
+ dummy_omero_metadata.update(
+ {
+ "channels": [
+ parameters["ref_channel"],
+ *(["nil"] * (dummy_omero_metadata["size_c"] - 1)),
+ ],
+ "name": "",
+ }
+ )
+
+ return cls(
+ imgdmy_obj.data,
+ dummy_omero_metadata,
+ TilerParameters.from_dict(parameters),
+ )
+
+ @classmethod
+ def from_image(cls, image, parameters: TilerParameters):
"""
Instantiate Tiler from an Image instance
@@ -262,7 +299,9 @@ class Tiler(ProcessABC):
@classmethod
def from_h5(
cls,
- image: t.Union[Image, ImageLocalOME, ImageDir],
+ image: t.Union[
+ ImageLocalOME, ImageDir
+ ], # TODO provide baseclass instead
filepath: t.Union[str, PosixPath],
parameters: TilerParameters = None,
):
@@ -279,8 +318,6 @@ class Tiler(ProcessABC):
trap_locs = TrapLocations.read_hdf5(filepath)
metadata = BridgeH5(filepath).meta_h5
metadata["channels"] = image.metadata["channels"]
- # metadata["zsectioning/nsections"] = image.metadata["zsectioning/nsections"]
- # metadata["channels/zsect"] = image.metadata["channels/zsect"]
if parameters is None:
parameters = TilerParameters.default()
tiler = cls(
@@ -321,15 +358,10 @@ class Tiler(ProcessABC):
@property
def shape(self):
"""
- Returns properties of the time-lapse experiment
- no of channels
- no of time points
- no of z stacks
- no of pixels in y direction
- no of pixels in z direction
+ Returns properties of the time-lapse as shown by self.image.shape
+
"""
- c, t, z, y, x = self.image.shape
- return (c, t, x, y, z)
+ return self.image.shape
@property
def n_processed(self):
@@ -453,7 +485,7 @@ class Tiler(ProcessABC):
ndtrap = self.ifoob_pad(full, trap.as_range(tp))
return ndtrap
- def run_tp(self, tp):
+ def _run_tp(self, tp):
"""
Find traps if they have not yet been found.
Determine any translational drift of the current image from the
@@ -491,9 +523,8 @@ class Tiler(ProcessABC):
return None
def get_traps_timepoint(self, *args, **kwargs):
- #
- print(
- DeprecationWarning("Deprecated:Use get_tiles_timepoint instead.")
+ self._log(
+ "get_trap_timepoints is deprecated; get_tiles_timepoint instead."
)
return self.get_tiles_timepoint(*args, **kwargs)
@@ -622,7 +653,7 @@ def find_channel_index(image_channels: t.List[str], channel: str):
found = re.match(channel, ch, re.IGNORECASE)
if found:
if len(found.string) - (found.endpos - found.start()):
- print(f"WARNING: channel {channel} matched {ch} using regex")
+ self._log(f"Channel {channel} matched {ch} using regex")
return i
diff --git a/src/aliby/tile/traps.py b/src/aliby/tile/traps.py
index 8d1d776c8bbc17a5c843fa3d145c324a20baba5e..4eddeb7e45a0f39ea0de28c865b79b685061da5b 100644
--- a/src/aliby/tile/traps.py
+++ b/src/aliby/tile/traps.py
@@ -140,9 +140,6 @@ def segment_traps(
return traps_retry
-###
-
-
def identify_trap_locations(
image, trap_template, optimize_scale=True, downscale=0.35, trap_size=None
):
@@ -240,103 +237,10 @@ def identify_trap_locations(
def stretch_image(image):
+ # FIXME Used in aliby.utils.imageViewer
image = ((image - image.min()) / (image.max() - image.min())) * 255
minval = np.percentile(image, 2)
maxval = np.percentile(image, 98)
image = np.clip(image, minval, maxval)
image = (image - minval) / (maxval - minval)
- return image
-
-
-def get_tile_shapes(x, tile_size):
- half_size = tile_size // 2
- xmin = int(x[0] - half_size)
- ymin = max(0, int(x[1] - half_size))
-
- return xmin, xmin + tile_size, ymin, ymin + tile_size
-
-
-def in_image(img, xmin, xmax, ymin, ymax, xidx=2, yidx=3):
- if xmin >= 0 and ymin >= 0:
- if xmax < img.shape[xidx] and ymax < img.shape[yidx]:
- return True
- else:
- return False
-
-
-def get_xy_tile(img, xmin, xmax, ymin, ymax, xidx=2, yidx=3, pad_val=None):
- if pad_val is None:
- pad_val = np.median(img)
- # Get the tile from the image
- idx = [slice(None)] * len(img.shape)
- idx[xidx] = slice(max(0, xmin), min(xmax, img.shape[xidx]))
- idx[yidx] = slice(max(0, ymin), min(ymax, img.shape[yidx]))
- tile = img[tuple(idx)]
- # Check if the tile is in the image
- if in_image(img, xmin, xmax, ymin, ymax, xidx, yidx):
- return tile
- else:
- # Add padding
- pad_shape = [(0, 0)] * len(img.shape)
- pad_shape[xidx] = (max(-xmin, 0), max(xmax - img.shape[xidx], 0))
- pad_shape[yidx] = (max(-ymin, 0), max(ymax - img.shape[yidx], 0))
- tile = np.pad(tile, pad_shape, constant_values=pad_val)
- return tile
-
-
-def tile_where(centre, x, y, MAX_X, MAX_Y):
- # Find the position of the tile
- xmin = int(centre[1] - x // 2)
- ymin = int(centre[0] - y // 2)
- xmax = xmin + x
- ymax = ymin + y
- # What do we actually have available?
- r_xmin = max(0, xmin)
- r_xmax = min(MAX_X, xmax)
- r_ymin = max(0, ymin)
- r_ymax = min(MAX_Y, ymax)
- return xmin, ymin, xmax, ymax, r_xmin, r_ymin, r_xmax, r_ymax
-
-
-def get_tile(shape, center, raw_expt, ch, t, z):
- """Returns a tile from the raw experiment with a given shape.
-
- :param shape: The shape of the tile in (C, T, Z, Y, X) order.
- :param center: The x,y position of the centre of the tile
- :param
- """
- _, _, x, y, _ = shape
- _, _, MAX_X, MAX_Y, _ = raw_expt.shape
- tile = np.full(shape, np.nan)
-
- # Find the position of the tile
- xmin = int(center[1] - x // 2)
- ymin = int(center[0] - y // 2)
- xmax = xmin + x
- ymax = ymin + y
- # What do we actually have available?
- r_xmin = max(0, xmin)
- r_xmax = min(MAX_X, xmax)
- r_ymin = max(0, ymin)
- r_ymax = min(MAX_Y, ymax)
-
- # Fill values
- tile[
- :,
- :,
- (r_xmin - xmin) : (r_xmax - xmin),
- (r_ymin - ymin) : (r_ymax - ymin),
- :,
- ] = raw_expt[ch, t, r_xmin:r_xmax, r_ymin:r_ymax, z]
- # fill_val = np.nanmedian(tile)
- # np.nan_to_num(tile, nan=fill_val, copy=False)
- return tile
-
-
-def centre(img, percentage=0.3):
- y, x = img.shape
- cropx = int(np.ceil(x * percentage))
- cropy = int(np.ceil(y * percentage))
- startx = int(x // 2 - (cropx // 2))
- starty = int(y // 2 - (cropy // 2))
- return img[starty : starty + cropy, startx : startx + cropx]
+ return image
\ No newline at end of file
diff --git a/src/aliby/track/__init__.py b/src/aliby/track/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..643e3ab57855d4ad5d1742a491806f2a9d45bb0c
--- /dev/null
+++ b/src/aliby/track/__init__.py
@@ -0,0 +1,3 @@
+"""
+Classes that link outlines within and between time-points.
+"""
diff --git a/src/aliby/track/benchmark.py b/src/aliby/track/benchmark.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c671d033d349c7ff413fc28210972ce13c9ed55
--- /dev/null
+++ b/src/aliby/track/benchmark.py
@@ -0,0 +1,382 @@
+import numpy as np
+import pandas as pd
+from scipy.ndimage import binary_fill_holes
+
+from baby.io import load_tiled_image
+from baby.tracker.core import CellTracker
+
+
+class CellBenchmarker: # TODO Simplify this by inheritance
+ """
+ Takes a metadata dataframe and a model and estimates the prediction in a trap-wise manner.
+
+ This class can also produce confusion matrices for a given Tracker and validation dataset.
+ """
+
+ def __init__(self, meta, model, bak_model, nstepsback=None):
+ self.indices = ["experimentID", "position", "trap", "tp"]
+ self.cindices = self.indices + ["cellLabels"]
+ self.meta = meta.copy()
+ self.meta["cont_list_index"] = (*range(len(self.meta)),)
+
+ self.tracker = CellTracker(model=model, bak_model=bak_model)
+ if nstepsback is None:
+ self.nstepsback = self.tracker.nstepsback
+ self.traps_loc
+
+ @property
+ def traps_loc(self):
+ """
+ Generates a list of trap locations using the metadata.
+ """
+ if not hasattr(self, "_traps_loc"):
+
+ traps = np.unique(
+ [
+ ind[: self.indices.index("trap") + 1]
+ for ind in self.meta.index
+ ],
+ axis=0,
+ )
+ # str->int conversion
+ traps = [(ind[0], *map(int, ind[1:])) for ind in traps]
+ self._traps_loc = (*map(tuple, traps),)
+
+ return self._traps_loc
+
+ @property
+ def masks(self):
+ if not hasattr(self, "_masks"):
+ self._masks = [
+ load_tiled_image(fname)[0] for fname in self.meta["filename"]
+ ]
+ for i, mask in enumerate(self._masks):
+ for j in range(mask.shape[2]):
+ self._masks[i][..., j] = binary_fill_holes(
+ self._masks[i][..., j]
+ )
+
+ self._masks = [np.moveaxis(mask, 2, 0) for mask in self._masks]
+
+ return self._masks
+
+ def predict_lbls_from_tpimgs(self, tp_img_tuple):
+ max_lbl = 0
+ prev_feats = []
+ cell_lbls = []
+ for tp, masks in tp_img_tuple:
+ lastn_lbls = cell_lbls[-self.nstepsback :]
+ lastn_feats = prev_feats[-self.nstepsback :]
+ new_lbl, feats, max_lbl = self.tracker.get_new_lbls(
+ masks, lastn_lbls, lastn_feats, max_lbl
+ )
+ cell_lbls = cell_lbls + [new_lbl]
+ prev_feats = prev_feats + [feats]
+
+ return (tp, cell_lbls)
+
+ def df_get_imglist(self, exp, pos, trap, tp=None):
+ df = self.meta.loc[(exp, pos, trap), ["cont_list_index", "cellLabels"]]
+ return zip(df.index, [self.masks[i] for i in df["cont_list_index"]])
+
+ def predict_set(self, exp, pos, trap, tp=None):
+ """
+ Predict labels using tp1-tp2 accuracy of prediction
+ """
+ # print("Processing trap {}".format(exp, pos, trap))
+ tp_img_tuple = (*self.df_get_imglist(exp, pos, trap),)
+ tp, lbl_list = self.predict_lbls_from_tpimgs(tp_img_tuple)
+ # print("loc {}, {}, {}, labels: {}".format(exp, pos, trap, lbl_list))
+ return lbl_list
+
+ def compare_traps(self, exp, pos, trap):
+ """
+ Error calculator for testing model and assignment heuristics.
+
+ Uses the trap id to compare the amount of cells correctly predicted.
+ This uses local indices, not whole timepoints. It returns the
+ fraction of cells correctly predicted, and the timepoints of mistakes
+
+ Returns:
+ float: Fraction of cells correctly predicted
+ list of 2-sized tuples: list of tp id of errors and the mistaken cell
+
+ """
+ print("Processing trap {}, {}, {}".format(exp, pos, trap))
+ new_cids = self.predict_set(exp, pos, trap)
+
+ test_df = self.meta.loc(axis=0)[(exp, pos, trap)].copy()
+ test_df["pred_cellLabels"] = new_cids
+
+ orig = test_df["cellLabels"].values
+ new = test_df["pred_cellLabels"].values
+ local_indices = [[], []]
+
+ # Case just defines if it is the test or new set
+ # print("Making tp-wise comparison")
+ for i, case in enumerate(
+ (zip(orig[:-1], orig[1:]), zip(new[:-1], new[1:]))
+ ):
+ for prev_cells, pos_cells in case:
+ local_assignment = [
+ prev_cells.index(cell) if cell in prev_cells else -1
+ for cell in pos_cells
+ ]
+ local_indices[i] += local_assignment
+
+ # Flatten
+ if len(local_indices) > 2:
+ flt_test, flt_new = [
+ np.array([j for i in case for j in i])
+ for case in local_indices
+ ]
+ tp_list = np.array(
+ [i for i, vals in enumerate(local_indices[0]) for j in vals]
+ )
+ else:
+ flt_test, flt_new = [
+ np.array([i for i in case]) for case in local_indices
+ ]
+ # tp_list = np.array(
+ # [i for i, vals in enumerate(local_indices[0]) for j in vals])
+
+ correct = flt_test == flt_new
+ if len(local_indices) > 2:
+ error_list = tp_list[~correct]
+ error_cid = (
+ test_df.iloc[1:]["cellLabels"].explode().dropna()[~correct].values
+ )
+ frac_correct = np.mean(correct)
+
+ print("Fraction of correct predictions", frac_correct)
+ if len(local_indices) > 2:
+ return (frac_correct, list(zip(error_list, error_cid)))
+ else:
+ # print("Warning: Single set of tps for this position")
+ return (frac_correct, error_cid)
+
+ def predict_all(self):
+ """
+ Predict all datasets defined in self.traps_loc
+ """
+ stepsback = [2]
+ threshs = [0.9]
+ self.predictions = {}
+ for nstepsback in stepsback:
+ for thresh in threshs:
+ self.nstepsback = nstepsback
+ self.tracker.nstepsback = nstepsback
+ self.low_thresh = 1 - thresh
+ self.high_thresh = thresh
+ self.thresh = thresh * 5 / 8
+ for address in self.traps_loc:
+ self.predictions[
+ (nstepsback, thresh, address)
+ ] = self.predict_set(*address)
+
+ def calculate_errsum(self):
+ """
+ Calculate all errors, addresses of images with errors and error fractions.
+ """
+ frac_errs = {}
+ all_errs = {}
+ nerrs = {}
+ stepsback = list(range(1, 3))
+ threshs = [0.95]
+ for nstepsback in stepsback:
+ for thresh in threshs:
+ self.nstepsback = nstepsback
+ self.tracker.nstepsback = nstepsback
+ self.low_thresh = 1 - thresh
+ self.high_thresh = thresh
+ self.thresh = thresh * 5 / 8
+ all_errs[(thresh, nstepsback)] = {}
+ frac_errs[(thresh, nstepsback)] = []
+ nerrs[(thresh, nstepsback)] = []
+ for address in self.traps_loc:
+ fraction, errors = self.compare_traps(*address)
+ if len(errors):
+ all_errs[(thresh, nstepsback)][address] = errors
+ frac_errs[(thresh, nstepsback)].append(fraction)
+ nerrs[(thresh, nstepsback)].append(len(errors))
+ else:
+ nerrs[(thresh, nstepsback)].append(0)
+ frac_errs[(thresh, nstepsback)].append(1.0)
+
+ return (frac_errs, all_errs, nerrs)
+
+ def get_truth_matrix_from_pair(self, pair):
+ """
+ Requires self.meta
+
+ args:
+ :pair: tuple of size 4 (experimentID, position, trap (tp1, tp2))
+
+ returns
+
+ :truth_mat: boolean ndarray of shape (ncells(tp1) x ncells(tp2)
+ links cells in tp1 to cells in tp2
+ """
+
+ clabs1 = self.meta.loc[pair[:3] + (pair[3][0],), "cellLabels"]
+ clabs2 = self.meta.loc[pair[:3] + (pair[3][1],), "cellLabels"]
+
+ truth_mat = gen_boolmat_from_clabs(clabs1, clabs2)
+
+ return truth_mat
+
+ def get_mota_stats(self, pair):
+ true_mat = self.get_truth_matrix_from_pair(pair)
+ prob_mat = self.tracker.predict_proba_from_ndarray(
+ ndarray, *args, **kwargs
+ )
+ pred_mat = prob_mat > thresh
+
+ true_flat = true_mat.flatten()
+ pred_flat = pred_mat.flatten()
+
+ true_pos = np.sum(true_flat & pred_flat)
+ false_pos = np.sum(true_flat & ~pred_flat)
+
+ # TODO add identity switch
+
+ def gen_cm_stats(self, pair, thresh=0.7, *args, **kwargs):
+ """
+ Calculate confusion matrix for a pair of pos-timepoints
+ """
+
+ masks = [self.masks[i] for i in self.meta.loc[pair, "cont_list_index"]]
+ feats = [self.tracker.calc_feats_from_mask(mask) for mask in masks]
+ ndarray = self.tracker.calc_feat_ndarray(*feats)
+ self.tracker.low_thresh = 1 - thresh
+ self.tracker.high_thresh = thresh
+ prob_mat = self.tracker.predict_proba_from_ndarray(
+ ndarray, *args, **kwargs
+ )
+ pred_mat = prob_mat > thresh
+
+ true_mat = self.get_truth_matrix_from_pair(pair)
+ if not len(true_mat) and not len(pred_mat):
+ return (0, 0, 0, 0)
+
+ true_flat = true_mat.flatten()
+ pred_flat = pred_mat.flatten()
+
+ true_pos = np.sum(true_flat & pred_flat)
+ false_pos = np.sum(true_flat & ~pred_flat)
+ false_neg = np.sum(~true_flat & pred_flat)
+ true_neg = np.sum(~true_flat & ~pred_flat)
+
+ return (true_pos, false_pos, false_neg, true_neg)
+
+ def extract_pairs_from_trap(self, trap_loc):
+ subdf = self.meta[["list_index", "cellLabels"]].loc(axis=0)[trap_loc]
+ pairs = [
+ trap_loc + tuple((pair,))
+ for pair in zip(subdf.index[:-1], subdf.index[1:])
+ ]
+
+ return pairs
+
+ def gen_pairlist(self):
+ self.pairs = [
+ self.extract_pairs_from_trap(trap) for trap in self.traps_loc
+ ]
+
+ def gen_cm_from_pairs(self, thresh=0.5, *args, **kwargs):
+ con_mat = {}
+ con_mat["tp"] = 0
+ con_mat["fp"] = 0
+ con_mat["fn"] = 0
+ con_mat["tn"] = 0
+ for pairset in self.pairs:
+ for pair in pairset:
+ res = self.gen_cm_stats(pair, thresh=thresh, *args, **kwargs)
+ con_mat["tp"] += res[0]
+ con_mat["fp"] += res[1]
+ con_mat["fn"] += res[2]
+ con_mat["tn"] += res[3]
+ self._con_mat = con_mat
+ return self._con_mat
+
+ def get_frac_error_df(self):
+ """
+ Calculates the trap-wise error and averages across a position.
+ """
+ self.frac_errs, self.all_errs, self.nerrs = self.calculate_errsum()
+
+ # nerrs_df = pd.DataFrame(self.nerrs).melt()
+ frac_df = pd.DataFrame(self.frac_errs).melt()
+
+ return frac_df
+
+ def gen_errorplots(self):
+ """
+ Calculates the trap-wise error and averages across a position.
+ """
+ frac_df = self.get_frac_error_df()
+ import seaborn as sns
+ from matplotlib import pyplot as plt
+
+ # ax = sns.barplot(x='variable_0', y='value', data=frac_df)
+ ax = sns.barplot(
+ x="variable_1", y="value", hue="variable_0", data=frac_df
+ )
+ ax.set(
+ xlabel="Backtrace depth",
+ ylabel="Fraction of correct assignments",
+ ylim=(0.9, 1),
+ )
+ plt.legend(title="Threshold")
+ plt.savefig("tracker_benchmark_btdepth.png")
+ plt.show()
+
+ # def plot_pair(self, address)
+
+
+def gen_boolmat_from_clabs(clabs1, clabs2):
+ if not np.any(clabs1) and not np.any(clabs2):
+ return np.array([])
+
+ boolmat = np.zeros((len(clabs1), len(clabs2))).astype(bool)
+ for i, lab1 in enumerate(clabs1):
+ for j, lab2 in enumerate(clabs2):
+ if lab1 == lab2:
+ boolmat[i, j] = True
+
+ return boolmat
+
+
+def gen_stats_dict(results):
+ """
+ Generates a dictionary using results from different binary classification tasks,
+ for example, using different thresholds
+
+ output
+
+ dictionary containing the name of statistic as a key and a list
+ of that statistic for the data subsets.
+ """
+ funs = (get_precision, get_recall, get_tnr, get_balanced_acc)
+ names = ("precision", "recall", "TNR", "balanced_acc")
+ stats_dict = {
+ name: [fun(res) for res in results] for fun, name in zip(funs, names)
+ }
+
+ return stats_dict
+
+
+def get_precision(res_dict):
+ return (res_dict["tp"]) / (res_dict["tp"] + res_dict["fp"])
+
+
+def get_recall(res_dict):
+ return res_dict["tp"] / (res_dict["tp"] + res_dict["fn"])
+
+
+def get_tnr(res_dict):
+ return res_dict["tn"] / (res_dict["tn"] + res_dict["fp"])
+
+
+def get_balanced_acc(res_dict):
+ return (get_recall(res_dict) + get_tnr(res_dict)) / 2
diff --git a/src/aliby/track/cell_tracker.py b/src/aliby/track/cell_tracker.py
new file mode 100644
index 0000000000000000000000000000000000000000..b51674409e06df09a964a87c97e868735023c316
--- /dev/null
+++ b/src/aliby/track/cell_tracker.py
@@ -0,0 +1,672 @@
+#!/usr/bin/env python
+"""
+TrackerCoordinator class to coordinate cell tracking and bud assignment.
+"""
+import pickle
+import typing as t
+from collections import Counter
+from os.path import dirname, join
+from pathlib import Path, PosixPath
+
+import numpy as np
+from scipy.optimize import linear_sum_assignment
+from sklearn.ensemble import RandomForestClassifier
+from sklearn.svm import SVC
+
+from agora.track_abc import FeatureCalculator
+
+models_path = join(dirname(__file__), "../models")
+
+
+class CellTracker(FeatureCalculator):
+ """
+ Class used to manage cell tracking. You can call it using an existing model or
+ use the inherited CellTrainer to get a new one.
+
+ Initialization parameters:
+
+ :model: sklearn.ensemble.RandomForestClassifier object
+ :trapfeats: Features to manually calculate within a trap
+ :extrafeats: Additional features to calculate
+ :model: Model to use, if provided ignores all other args but threshs
+ :bak_model: Backup mode to use when prediction is unsure
+ :nstepsback: int Number of timepoints to go back
+ :thresh: float Cut-off value to assume a cell is not new
+ :low_thresh: Lower thresh for model switching
+ :high_thresh: Higher thresh for model switching.
+ Probabilities between these two thresholds summon the
+ backup model.
+
+ :aweights: area weight for barycentre calculations
+
+
+ # Feature order in array features
+ 1. basic features
+ 2. trap features (within a trap)
+ 3. extra features (process between imgs)
+
+ """
+
+ def __init__(
+ self,
+ feats2use=None,
+ trapfeats=None,
+ extrafeats=None,
+ model=None,
+ bak_model=None,
+ thresh=None,
+ low_thresh=None,
+ high_thresh=None,
+ nstepsback=None,
+ aweights=None,
+ red_fun=None,
+ max_distance=None,
+ **kwargs,
+ ):
+
+ if trapfeats is None:
+ trapfeats = ()
+
+ if extrafeats is None:
+ extrafeats = ()
+
+ if type(model) is str or type(model) is PosixPath:
+ with open(Path(model), "rb") as f:
+ model = pickle.load(f)
+
+ if type(bak_model) is str or type(bak_model) is PosixPath:
+ with open(Path(bak_model), "rb") as f:
+ bak_model = pickle.load(f)
+
+ if aweights is None:
+ self.aweights = None
+
+ if feats2use is None: # Ignore this block when training
+ if model is None:
+ model = self.load_model(models_path, "ct_XGBC_20220703_11.pkl")
+ if bak_model is None:
+ bak_model = self.load_model(
+ models_path, "ct_XGBC_20220703_10.pkl"
+ )
+ self.model = model
+ self.bak_model = bak_model
+
+ main_feats = model.all_ifeats
+ bak_feats = bak_model.all_ifeats
+ feats2use, trapfeats, extrafeats = [
+ tuple(sorted(set(main).union(bak)))
+ for main, bak in zip(main_feats, bak_feats)
+ ]
+
+ # Training AND non-training part
+ super().__init__(
+ feats2use, trapfeats=trapfeats, extrafeats=extrafeats, **kwargs
+ )
+
+ self.extrafeats = tuple(extrafeats)
+
+ self.all_ofeats = self.outfeats + trapfeats + extrafeats
+
+ self.noutfeats = len(self.all_ofeats)
+
+ if hasattr(self, "bak_model"): # Back to non-training only
+
+ self.mainof_ids = [
+ self.all_ofeats.index(f) for f in self.model.all_ofeats
+ ]
+
+ self.bakof_ids = [
+ self.all_ofeats.index(f) for f in self.bak_model.all_ofeats
+ ]
+
+ if nstepsback is None:
+ nstepsback = 3
+ self.nstepsback = nstepsback
+
+ if thresh is None:
+ thresh = 0.5
+ self.thresh = thresh
+ if low_thresh is None:
+ low_thresh = 0.4
+ if high_thresh is None:
+ high_thresh = 0.6
+ self.low_thresh, self.high_thresh = low_thresh, high_thresh
+
+ if red_fun is None:
+ red_fun = np.nanmax
+ self.red_fun = red_fun
+
+ if max_distance is None:
+ self.max_distance = max_distance
+
+ def get_feats2use(self):
+ """
+ Return feats to be used from a loaded random forest model model
+ """
+ nfeats = get_nfeats_from_model(self.model)
+ nfeats_bak = get_nfeats_from_model(self.bak_model)
+ # max_nfeats = max((nfeats, nfeats_bak))
+
+ return (switch_case_nfeats(nfeats), switch_case_nfeats(nfeats_bak))
+
+ def calc_feat_ndarray(self, prev_feats, new_feats):
+ """
+ Calculate feature ndarray using two ndarrays of features.
+ ---
+
+ input
+
+ :prev_feats: ndarray (ncells, nfeats) of timepoint 1
+ :new_feats: ndarray (ncells, nfeats) of timepoint 2
+
+ returns
+
+ :n3darray: ndarray (ncells_prev, ncells_new, nfeats) containing a
+ cell-wise substraction of the features in the input ndarrays.
+ """
+ if not (new_feats.any() and prev_feats.any()):
+ return np.array([])
+
+ n3darray = np.empty((len(prev_feats), len(new_feats), self.ntfeats))
+
+ # print('self: ', self, ' self.ntfeats: ', self.ntfeats, ' featsshape: ', new_feats.shape)
+ for i in range(self.ntfeats):
+ n3darray[..., i] = np.subtract.outer(
+ prev_feats[:, i], new_feats[:, i]
+ )
+
+ n3darray = self.calc_dtfeats(n3darray)
+
+ return n3darray
+
+ def calc_dtfeats(self, n3darray):
+ """
+ Calculates features obtained between timepoints, such as distance
+ for every pair of cells from t1 to t2.
+ ---
+
+ input
+
+ :n3darray: ndarray (ncells_prev, ncells_new, ntfeats) containing a
+ cell-wise substraction of the features in the input ndarrays.
+
+ returns
+
+ :newarray: 3d array taking the features specified in self.outfeats and self.trapfeats
+ and adding dtfeats
+ """
+ newarray = np.empty(n3darray.shape[:-1] + (self.noutfeats,))
+ newarray[..., : len(self.outfeats)] = n3darray[
+ ..., : len(self.outfeats)
+ ]
+ newarray[
+ ..., len(self.outfeats) : len(self.outfeats) + len(self.trapfeats)
+ ] = n3darray[..., len(self.out_merged) :]
+ for i, feat in enumerate(self.all_ofeats):
+ if feat == "distance":
+ newarray[..., i] = np.sqrt(
+ n3darray[..., self.xind] ** 2
+ + n3darray[..., self.yind] ** 2
+ )
+
+ return newarray
+
+ def assign_lbls(
+ self,
+ prob_backtrace: np.ndarray,
+ prev_lbls: t.List[t.List[int]],
+ red_fun=None,
+ ):
+ """Assign labels using a prediction matrix of nxmxl where n is the number
+ of cells in the previous image, m the number of steps back considered
+ and l in the new image. It assigns the
+ number zero if it doesn't find the cell.
+ ---
+ input
+
+ :prob_backtrace: Probability n x m x l array obtained as an output of rforest
+ :prev_labels: List of cell labels for previous timepoint to be compared.
+ :red_fun: Function used to collapse the previous timepoints into one.
+ If none provided it uses maximum and ignores np.nans.
+
+ returns
+
+ :new_lbls: ndarray of newly assigned labels obtained, new cells as
+ zero.
+ """
+ if red_fun is None:
+ red_fun = self.red_fun
+
+ new_lbls = np.zeros(prob_backtrace.shape[2], dtype=int)
+ pred_matrix = red_fun(prob_backtrace, axis=1)
+
+ if pred_matrix.any():
+ # assign available hits
+ row_ids, col_ids = linear_sum_assignment(-pred_matrix)
+ for i, j in zip(row_ids, col_ids):
+ if pred_matrix[i, j] > self.thresh:
+ new_lbls[j] = prev_lbls[i]
+
+ return new_lbls
+
+ def predict_proba_from_ndarray(
+ self,
+ array_3d: np.ndarray,
+ model: str = None,
+ boolean: bool = False,
+ max_distance: float = None,
+ ):
+ """
+
+ input
+
+ :array_3d: (ncells_tp1, ncells_tp2, out_feats) ndarray
+ :model: str, {'model', 'bak_model'} can force a unique model instead of an ensemble
+ :boolean: bool, if False returns probability, if True returns prediction
+ :max_distance: float Maximum distance (in um) to be considered. If None it uses the instance's value,
+ if zero it skips checking distances.
+
+ requires
+ :self.model:
+ :self.mainof_ids: list of indices corresponding to the main model's features
+ :self.bakof_ids: list of indices corresponding to the backup model's features
+
+ returns
+
+ (ncells_tp1, ncells_tp2) ndarray with probabilities or prediction
+ of cell identities depending on "boolean" arg.
+
+ """
+
+ if array_3d.size == 0:
+ return np.array([])
+
+ if model is None:
+ model2use = self.model
+ bak_model2use = self.bak_model
+ bakof_ids = self.bakof_ids
+ mainof_ids = self.mainof_ids
+ else:
+ model2use = getattr(self, "model")
+ bak_model2use = model2use
+ bakof_ids = [
+ self.all_ofeats.index(f) for f in model2use.all_ofeats
+ ]
+ mainof_ids = [
+ self.all_ofeats.index(f) for f in model2use.all_ofeats
+ ]
+
+ fun2use = "predict" if boolean else "predict_proba"
+ predict_fun = getattr(model2use, fun2use)
+ bak_pred_fun = getattr(bak_model2use, fun2use)
+
+ if max_distance is None:
+ max_distance = self.max_distance
+
+ orig_shape = array_3d.shape[:2]
+
+ # Ignore cells that are too far away to possibly be the same
+ cells_near = np.ones(orig_shape, dtype=bool)
+ if max_distance and set(self.all_ofeats).intersection(
+ ("distance", "centroid-0", "centroid-1")
+ ):
+ if "distance" in self.all_ofeats:
+ cells_near = (
+ array_3d[..., self.all_ofeats.index("distance")]
+ < max_distance
+ )
+ else: # Calculate euclidean distance
+ cells_near = (
+ np.sqrt(
+ array_3d[..., self.all_ofeats.index("centroid-0")]
+ + array_3d[..., self.all_ofeats.index("centroid-1")]
+ )
+ < max_distance
+ )
+
+ pred_matrix = np.zeros(orig_shape)
+ prob = np.zeros(orig_shape)
+ if cells_near.any():
+ prob = predict_fun(array_3d[cells_near][:, mainof_ids])[:, 1]
+ uncertain_dfeats = (self.low_thresh < prob) & (
+ prob < self.high_thresh
+ )
+ if uncertain_dfeats.any():
+ bak_prob = bak_pred_fun(
+ array_3d[cells_near][uncertain_dfeats][:, bakof_ids]
+ )[:, 1]
+ probs_compared = np.stack((prob[uncertain_dfeats], bak_prob))
+ most_confident_proba = abs((probs_compared - 0.5)).argmax(
+ axis=0
+ )
+ prob[uncertain_dfeats] = probs_compared[
+ most_confident_proba, range(most_confident_proba.shape[0])
+ ]
+
+ pred_matrix[cells_near] = prob
+
+ return pred_matrix
+
+ def get_new_lbls(
+ self,
+ new_img,
+ prev_lbls,
+ prev_feats,
+ max_lbl,
+ new_feats=None,
+ pixel_size=None,
+ **kwargs,
+ ):
+ """
+ Core function to calculate the new cell labels.
+
+ ----
+
+ input
+
+ :new_img: ndarray (len, width, ncells) containing the cell outlines
+ :max_lbl: int indicating the last assigned cell label
+ :prev_feats: list of ndarrays of size (ncells x nfeatures)
+ containing the features of previous timepoints
+ :prev_lbls: list of list of ints corresponding to the cell labels in
+ the previous timepoints
+ :new_feats: (optional) Directly give a feature ndarray. It ignores
+ new_img if given.
+ :kwargs: Additional keyword values passed to self.predict_proba_from_ndarray
+
+ returns
+
+ :new_lbls: list of labels assigned to new timepoint
+ :new_feats: list of ndarrays containing the updated features
+ :new_max: updated max cell label assigned
+
+ """
+
+ if new_feats is None:
+ new_feats = self.calc_feats_from_mask(new_img)
+
+ if new_feats.any():
+ if np.any([len(prev_feat) for prev_feat in prev_feats]):
+ counts = Counter(
+ [lbl for lbl_set in prev_lbls for lbl in lbl_set]
+ )
+ lbls_order = list(counts.keys())
+ probs = np.full(
+ (len(lbls_order), self.nstepsback, len(new_feats)), np.nan
+ )
+
+ for i, (lblset, prev_feat) in enumerate(
+ zip(prev_lbls, prev_feats)
+ ):
+ if len(prev_feat):
+ feats_3darray = self.calc_feat_ndarray(
+ prev_feat, new_feats
+ )
+
+ pred_matrix = self.predict_proba_from_ndarray(
+ feats_3darray,
+ **kwargs,
+ )
+
+ for j, lbl in enumerate(lblset):
+ probs[lbls_order.index(lbl), i, :] = pred_matrix[
+ j, :
+ ]
+
+ new_lbls = self.assign_lbls(probs, lbls_order)
+ new_cells_pos = new_lbls == 0
+ new_max = max_lbl + sum(new_cells_pos)
+ new_lbls[new_cells_pos] = [*range(max_lbl + 1, new_max + 1)]
+
+ # ensure that label output is consistently a list
+ new_lbls = new_lbls.tolist()
+
+ else:
+ new_lbls = [*range(max_lbl + 1, max_lbl + len(new_feats) + 1)]
+
+ new_max = max_lbl + len(new_feats)
+
+ else:
+ return ([], [], max_lbl)
+ return (new_lbls, new_feats, new_max)
+
+ def probabilities_from_impair(
+ self,
+ image_t1: np.ndarray,
+ image_t2: np.ndarray,
+ kwargs_feat_calc: dict = {},
+ **kwargs,
+ ):
+ """
+ Convenience function to test tracking between two time-points
+
+ :image_t1: np.ndarray containing mask of first time-point
+ :image_t2: np.ndarray containing mask of second time-point
+ :kwargs_feat_calc: are passed to self.calc_feats_from_mask calls
+ :kwargs: are passed to self.predict_proba_from_ndarray
+ """
+ feats_t1 = self.calc_feats_from_mask(image_t1, **kwargs_feat_calc)
+ feats_t2 = self.calc_feats_from_mask(image_t2, **kwargs_feat_calc)
+
+ probability_matrix = np.array([])
+ if feats_t1.any() and feats_t2.any():
+ feats_3darray = self.calc_feat_ndarray(feats_t1, feats_t2)
+
+ probability_matrix = self.predict_proba_from_ndarray(
+ feats_3darray, **kwargs
+ )
+
+ return probability_matrix
+
+ # Step CellTracker
+ def run_tp(
+ self,
+ masks: np.ndarray,
+ state: t.Dict[str, t.Union[int, list]] = None,
+ **kwargs,
+ ) -> t.Dict[str, t.Union[t.List[int], t.Dict[str, t.Union[int, list]]]]:
+ """Assign labels to new masks using a state dictionary.
+
+ Parameters
+ ----------
+ masks : np.ndarray
+ Cell masks to label.
+ state : t.Dict[str, t.Union[int, list]]
+ Dictionary containing maximum cell label, and previous cell labels
+ and features for those cells.
+ kwargs : keyword arguments
+ Keyword arguments passed to self.get_new_lbls
+
+ Returns
+ -------
+ t.Dict[str, t.Union[t.List[int], t.Dict[str, t.Union[int, list]]]]
+ New labels and new state dictionary.
+
+ Examples
+ --------
+ FIXME: Add example beyond the trivial one.
+
+ import numpy as np
+ from baby.tracker.core import celltracker
+ from tqdm import tqdm
+
+ # overlapping outlines are of shape (t,z,x,y)
+ masks = np.zeros((5, 3, 20, 20), dtype=bool)
+ masks[0, 0, 2:6, 2:6] = true
+ masks[1:, 0, 13:14, 13:14] = true
+ masks[:, 1, 8:12, 8:12] = true
+ masks[:, 2, 14:18, 14:18] = true
+
+ # 13um pixel size
+ ct = celltracker(pixel_size=0.185)
+ labels = []
+
+ state = none
+ for masks_tp in tqdm(masks):
+ new_labels, state = ct.run_tp(masks_tp, state=state)
+ labels.append(new_labels)
+
+ # should result in state['cell_lbls']
+ # [[1, 2, 3], [4, 2, 3], [4, 2, 3], [4, 2, 3], [4, 2, 3]]
+ """
+ if state is None:
+ state = {}
+
+ max_lbl = state.get("max_lbl", 0)
+ cell_lbls = state.get("cell_lbls", [])
+ prev_feats = state.get("prev_feats", [])
+
+ # Get features for cells at this time point
+ feats = self.calc_feats_from_mask(masks)
+
+ lastn_lbls = cell_lbls[-self.nstepsback :]
+ lastn_feats = prev_feats[-self.nstepsback :]
+
+ new_lbls, _, max_lbl = self.get_new_lbls(
+ masks, lastn_lbls, lastn_feats, max_lbl, feats, **kwargs
+ )
+
+ state = {
+ "max_lbl": max_lbl,
+ "cell_lbls": cell_lbls + [new_lbls],
+ "prev_feats": prev_feats + [feats],
+ }
+
+ return (new_lbls, state)
+
+
+# Helper functions
+
+
+def switch_case_nfeats(nfeats):
+ """
+ Convenience TEMPORAL function to determine whether to use distance/location
+ as a feature for tracking or not (nfeats=5 for no distance, 7 for distance)
+ input
+ number of feats
+
+ returns
+ list of main and extra feats based on the number of feats
+ """
+ main_feats = {
+ 4: [
+ ("area", "minor_axis_length", "major_axis_length", "bbox_area"),
+ (),
+ (),
+ ],
+ 5: [
+ (
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ "perimeter",
+ ),
+ (),
+ (),
+ ],
+ 6: [
+ (
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ "perimeter",
+ ),
+ (),
+ ("distance",),
+ ],
+ # Including centroid
+ # 7 : [('centroid', 'area', 'minor_axis_length', 'major_axis_length',
+ # 'bbox_area', 'perimeter'), () , ()],
+ 7: [
+ (
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ "perimeter",
+ ),
+ ("baryangle", "barydist"),
+ (),
+ ],
+ 8: [ # Minus centroid
+ (
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ # "eccentricity",
+ "equivalent_diameter",
+ # "solidity",
+ # "extent",
+ "orientation",
+ # "perimeter",
+ ),
+ ("baryangle", "barydist"),
+ (),
+ ],
+ 10: [ # Minus distance
+ (
+ "centroid",
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ # "eccentricity",
+ "equivalent_diameter",
+ # "solidity",
+ # "extent",
+ "orientation",
+ # "perimeter",
+ ),
+ ("baryangle", "barydist"),
+ (),
+ ],
+ 11: [ # Minus computationally-expensive features
+ (
+ "centroid",
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ # "eccentricity",
+ "equivalent_diameter",
+ # "solidity",
+ # "extent",
+ "orientation",
+ # "perimeter",
+ ),
+ ("baryangle", "barydist"),
+ ("distance",),
+ ],
+ 15: [ # All features
+ (
+ "centroid",
+ "area",
+ "minor_axis_length",
+ "major_axis_length",
+ "bbox_area",
+ "eccentricity",
+ "equivalent_diameter",
+ "solidity",
+ "extent",
+ "orientation",
+ "perimeter",
+ ),
+ ("baryangle", "barydist"),
+ ("distance",),
+ ],
+ }
+
+ assert nfeats in main_feats.keys(), "invalid nfeats"
+ return main_feats.get(nfeats, [])
+
+
+def get_nfeats_from_model(model) -> int:
+ if isinstance(model, SVC):
+ nfeats = model.support_vectors_.shape[-1]
+ elif isinstance(model, RandomForestClassifier):
+ nfeats = model.n_features_
+
+ return nfeats
diff --git a/src/aliby/track/utils.py b/src/aliby/track/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..36850f533a32d4621300b20cc04d6405ac8e986b
--- /dev/null
+++ b/src/aliby/track/utils.py
@@ -0,0 +1,136 @@
+# If you publish results that make use of this software or the Birth Annotator
+# for Budding Yeast algorithm, please cite:
+# Julian M J Pietsch, Alán Muñoz, Diane Adjavon, Ivan B N Clark, Peter S
+# Swain, 2021, Birth Annotator for Budding Yeast (in preparation).
+#
+#
+# The MIT License (MIT)
+#
+# Copyright (c) Julian Pietsch, Alán Muñoz and Diane Adjavon 2021
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to
+# deal in the Software without restriction, including without limitation the
+# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+# sell copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+# IN THE SOFTWARE.
+import typing as t
+
+import numpy as np
+
+
+# Calculate barycentre
+def calc_barycentre(centres, weights=None, **kwargs):
+ """
+ :centres: ndarray containing the (x,y) centres of each cell
+ :weights: (optional) list of weights to consider for each cell
+ """
+ if weights is None:
+ weights = np.ones_like(centres)
+
+ barycentre = np.average(centres, axis=0, weights=weights)
+ return barycentre
+
+
+# Calculate distance to center
+def calc_barydists(centres, bary, **kwargs):
+ """
+ Calculate distances to the barycentre
+ :centre: int (2,) tuple. Centre of cell
+ :bary: float (2,) tuple. Barycentre of image
+ """
+ vec2bary = centres - bary
+ dists = np.sqrt(np.sum(vec2bary**2, axis=1))
+
+ return dists
+
+
+# Calculate angle to center
+def calc_baryangles(centres, bary, areas=None, **kwargs):
+ """
+ Calculate angle using centre of cell and barycentre
+ :centre: int (2,) tuple. Centre of cell
+ :bary: float (2,) tuple. Barycentre of image
+ :anchor_cell: int Cell id to use as angle 0.
+ """
+
+ angles = []
+ vec2bary = centres - bary
+ angles = np.apply_along_axis(lambda x: np.arctan2(*x), 1, vec2bary)
+ if areas is not None:
+ anchor_cell = np.argmax(areas)
+ angles -= angles[anchor_cell]
+
+ return angles
+
+
+def pick_baryfun(key):
+ baryfuns = {"barydist": calc_barydists, "baryangle": calc_baryangles}
+ return baryfuns[key]
+
+
+## Tracking benchmark utils
+
+
+def lol_to_adj(cell_ids: t.List[t.List[int]]):
+ """
+ Convert a series list of lists with cell ids into a matrix
+ representing a graph.
+
+ Note that information is lost in the process, and a matrix can't be
+ turned back into a list of list by itself.
+
+ input
+
+ :lol: list of lists with cell ids
+
+ returns
+
+ :adj_matrix: (n, n) ndarray where n is the number of cells
+
+ """
+ n = len([y for x in cell_ids for y in x])
+ adj_mat = np.zeros((n, n))
+
+ prev = None
+ cur = 0
+ for c_ids_single_lst in cell_ids:
+ if not prev:
+ prev = c_ids_single_lst
+ else:
+ for i, el in enumerate(c_ids_single_lst):
+ prev_idx = prev.index(el) if el in prev else None
+ if prev_idx is not None:
+ adj_mat[cur + len(prev) + i, cur + prev_idx] = True
+ cur += len(c_ids_single_lst)
+
+ return adj_mat
+
+
+def compare_pred_truth_lols(prediction, truth):
+ """
+ input
+
+ :prediction: list of lists with predicted cell ids
+ :truth: list of lists with real cell ids
+
+ returns
+
+ number of diferences between equivalent truth matrices
+
+ """
+ adj_pred = lol_to_adj(prediction)
+ adj_truth = lol_to_adj(truth)
+
+ return int(((adj_pred - adj_truth) != 0).sum())
diff --git a/src/aliby/utils/argo.py b/src/aliby/utils/argo.py
index 24255065e089a4a6b091f7b821d2bf594d2f4b2e..b4d80c364ab7cf62396660e0f19bd92a04673332 100644
--- a/src/aliby/utils/argo.py
+++ b/src/aliby/utils/argo.py
@@ -282,10 +282,6 @@ class OmeroExplorer:
for k, v in self.cache.items()
}
- # @staticfunction
- # def number_of_X(logfile: str):
- # return re.findall("X", logfile)
-
def dset_count(
self,
dset: t.Union[int, _DatasetWrapper],
@@ -351,17 +347,6 @@ class Argo(OmeroExplorer):
super().__init__(*args, **kwargs)
-def get_creds():
- return (
- "upload",
- "***REMOVED***", # OMERO Password
- )
-
-
-def list_files(dset):
- return {x for x in dset.listAnnotations() if hasattr(x, "getFileName")}
-
-
def annot_from_dset(dset, kind):
v = [x for x in dset.listAnnotations() if hasattr(x, "getFileName")]
infname = kind if kind == "log" else kind.title()
@@ -374,7 +359,6 @@ def annot_from_dset(dset, kind):
except Exception as e:
print(f"Conversion from acquisition file failed: {e}")
return {}
-
return acq
@@ -393,6 +377,7 @@ def check_channels(acq, channels, _all=True):
def get_chs(exptype):
+ # TODO Documentation
exptypes = {
"dual_ph": ("GFP", "pHluorin405", "mCherry"),
"ph": ("GFP", "pHluorin405"),
@@ -403,6 +388,7 @@ def get_chs(exptype):
def load_annot_from_cache(exp_id, cache_dir="cache/"):
+ # TODO Documentation
if type(cache_dir) is not PosixPath:
cache_dir = Path(cache_dir)
@@ -428,16 +414,6 @@ def parse_annot(str_io, fmt):
return parser.parse(io.StringIO(str_io))
-def get_log_date(annot_sets):
- log = get_annot(annot_sets, "log")
- return log.get("date", None)
-
-
-def get_log_microscope(annot_sets):
- log = get_annot(annot_sets, "log")
- return log.get("microscope", None)
-
-
def get_annotsets(dset):
annot_files = [
annot.getFile()
@@ -457,12 +433,8 @@ def get_annotsets(dset):
return annot_sets
-# def has_tags(d, tags):
-# if set(tags).intersection(annot_from_dset(d, "log").get("omero_tags", [])):
-# return True
-
-
def load_acq(dset):
+ # TODO Documentation
try:
acq = annot_from_dset(dset, kind="acq")
return acq
@@ -472,6 +444,7 @@ def load_acq(dset):
def has_channels(dset, exptype):
+ # TODO Documentation
acq = load_acq(dset)
if acq:
return check_channels(acq, get_chs(exptype))
@@ -479,26 +452,8 @@ def has_channels(dset, exptype):
return
-# Custom functions
-def compare_dsets_voltages_exp(dsets):
- a = {}
- for d in dsets:
- try:
- acq = annot_from_dset(d, kind="acq")["channels"]
- a[d.getId()] = {
- k: (v, e)
- for k, v, e in zip(
- acq["channel"], acq["voltage"], acq["exposure"]
- )
- }
-
- except Exception as e:
- print(d, f"Data set voltage comparison did not work:{e}")
-
- return a
-
-
def get_logfile(dset):
+ # TODO Documentation
annot_file = [
annot.getFile()
for annot in dset.listAnnotations()
diff --git a/src/aliby/utils/cache.py b/src/aliby/utils/cache.py
deleted file mode 100644
index a256c6350932f495423fec4f67ade468b9d607a6..0000000000000000000000000000000000000000
--- a/src/aliby/utils/cache.py
+++ /dev/null
@@ -1,141 +0,0 @@
-"""
-Utility functions and classes
-"""
-import itertools
-import logging
-import operator
-from functools import partial, wraps
-from pathlib import Path
-from time import perf_counter
-from typing import Callable
-
-import cv2
-import h5py
-import imageio
-import numpy as np
-
-
-def repr_obj(obj, indent=0):
- """
- Helper function to display info about OMERO objects.
- Not all objects will have a "name" or owner field.
- """
- string = """%s%s:%s Name:"%s" (owner=%s)""" % (
- " " * indent,
- obj.OMERO_CLASS,
- obj.getId(),
- obj.getName(),
- obj.getAnnotation(),
- )
-
- return string
-
-
-def imread(path):
- return cv2.imread(str(path), -1)
-
-
-class ImageCache:
- """HDF5-based image cache for faster loading of the images once they've
- been read.
- """
-
- def __init__(self, file, name, shape, remote_fn):
- self.store = h5py.File(file, "a")
- # Create a dataset
- self.dataset = self.store.create_dataset(
- name, shape, dtype=np.float, fill_value=np.nan
- )
- self.remote_fn = remote_fn
-
- def __getitem__(self, item):
- cached = self.dataset[item]
- if np.any(np.isnan(cached)):
- full = self.remote_fn(item)
- self.dataset[item] = full
- return full
- else:
- return cached
-
-
-class Cache:
- """
- Fixed-length mapping to use as a cache.
- Deletes items in FIFO manner when maximum allowed length is reached.
- """
-
- def __init__(self, max_len=5000, load_fn: Callable = imread):
- """
- :param max_len: Maximum number of items in the cache.
- :param load_fn: The function used to load new items if they are not
- available in the Cache
- """
- self._dict = dict()
- self._queue = []
- self.load_fn = load_fn
- self.max_len = max_len
-
- def __getitem__(self, item):
- if item not in self._dict:
- self.load_item(item)
- return self._dict[item]
-
- def load_item(self, item):
- self._dict[item] = self.load_fn(item)
- # Clean up the queue
- self._queue.append(item)
- if len(self._queue) > self.max_len:
- del self._dict[self._queue.pop(0)]
-
- def clear(self):
- self._dict.clear()
- self._queue.clear()
-
-
-def accumulate(lst: list):
- lst = sorted(lst)
- it = itertools.groupby(lst, operator.itemgetter(0))
- for key, sub_iter in it:
- yield key, [x[1] for x in sub_iter]
-
-
-def get_store_path(save_dir, store, name):
- """Create a path to a position-specific store.
-
- This combines the name and the store's base name into a file path within save_dir.
- For example.
- >>> get_store_path('data', 'baby_seg.h5', 'pos001')
- Path(data/pos001baby_seg.h5')
-
- :param save_dir: The root directory in which to save the file, absolute
- path.
- :param store: The base name of the store
- :param name: The name of the position
- :return: Path(save_dir) / name+store
- """
- store = Path(save_dir) / store
- store = store.with_name(name + store.name)
- return store
-
-
-def parametrized(dec):
- def layer(*args, **kwargs):
- def repl(f):
- return dec(f, *args, **kwargs)
-
- return repl
-
- return layer
-
-
-@parametrized
-def timed(f, name=None):
- @wraps(f)
- def decorated(*args, **kwargs):
- t = perf_counter()
- res = f(*args, **kwargs)
- to_print = name or f.__name__
- logging.debug(f"Timing:{to_print}:{perf_counter() - t}s")
- return res
-
- return decorated
diff --git a/src/extraction/core/extractor.py b/src/extraction/core/extractor.py
index 25668e6a0bbd38aee9a6af8d01745352678aa3fe..eee3185a613eb19af8faa5b88566d76d1f197d53 100644
--- a/src/extraction/core/extractor.py
+++ b/src/extraction/core/extractor.py
@@ -1,4 +1,3 @@
-import logging
import typing as t
from time import perf_counter
from typing import List
@@ -6,7 +5,7 @@ from typing import List
import h5py
import numpy as np
import pandas as pd
-from agora.abc import ParametersABC, ProcessABC
+from agora.abc import ParametersABC, StepABC
from agora.io.cells import Cells
from agora.io.writer import Writer, load_attributes
@@ -89,7 +88,7 @@ class ExtractorParameters(ParametersABC):
return cls(**exparams_from_meta(meta))
-class Extractor(ProcessABC):
+class Extractor(StepABC):
"""
Apply a metric to cells identified in the tiles.
@@ -266,7 +265,7 @@ class Extractor(ProcessABC):
channel_ids = None
if z is None:
# gets the tiles data via tiler
- z: t.List[int] = list(range(self.tiler.shape[-1]))
+ z: t.List[int] = list(range(self.tiler.shape[-3]))
tiles = (
self.tiler.get_tiles_timepoint(
tp, channels=channel_ids, z=z, **kwargs
@@ -306,7 +305,7 @@ class Extractor(ProcessABC):
A two-tuple of a tuple of results and a tuple with the corresponding trap_id and cell labels
"""
if labels is None:
- print("Warning: No labels given. Sorting cells using index.")
+ self._log("No labels given. Sorting cells using index.")
cell_fun = True if metric in self._all_cell_funs else False
idx = []
@@ -455,6 +454,7 @@ class Extractor(ProcessABC):
An example is d["GFP"]["np_max"]["mean"][0], which gives a tuple of the calculated mean GFP fluorescence for all cells.
"""
+ # TODO Can we split the different extraction types into sub-methods to make this easier to read?
if tree is None:
# use default
tree: extraction_tree = self.params.tree
@@ -588,28 +588,24 @@ class Extractor(ProcessABC):
elif channel in self.img_bgsub:
return self.img_bgsub[channel]
- def run_tp(self, tp, **kwargs):
- """
- Wrapper to add compatibility with other steps of the pipeline.
- """
- return self.run(tps=[tp], **kwargs)
-
- def run(
+ def _run_tp(
self,
- tree=None,
tps: List[int] = None,
+ tree=None,
save=True,
**kwargs,
) -> dict:
"""
+ Wrapper to add compatibility with other steps of the pipeline.
+
Parameters
----------
- tree: dict
+ tps: list of int (optional)
+ Time points to include.
+ tree: dict (optional)
Nested dictionary indicating channels, reduction functions and
metrics to be used.
For example: {'general': {'None': ['area', 'volume', 'eccentricity']}}
- tps: list of int (optional)
- Time points to include.
save: boolean (optional)
If True, save results to h5 file.
kwargs: keyword arguments (optional)
@@ -624,6 +620,9 @@ class Extractor(ProcessABC):
tree = self.params.tree
if tps is None:
tps = list(range(self.meta["time_settings/ntimepoints"][0]))
+ elif isinstance(tps, int):
+ tps = [tps]
+
# store results in dict
d = {}
for tp in tps:
diff --git a/src/extraction/core/functions/io.py b/src/extraction/core/functions/io.py
deleted file mode 100644
index 3b377bdffe47c991264e262ca83a18c7dbd726b0..0000000000000000000000000000000000000000
--- a/src/extraction/core/functions/io.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from yaml import dump, load
-
-
-def dict_to_yaml(d, f):
- with open(f, "w") as f:
- dump(d, f)
-
-
-def add_attrs(hdfile, path, files):
- group = hdfile.create_group(path)
- for k, v in files:
- group.attrs[k] = v
diff --git a/src/extraction/core/functions/utils.py b/src/extraction/core/functions/utils.py
deleted file mode 100644
index dcae68c9bde33a86108cf5a6c5bda701c80420e5..0000000000000000000000000000000000000000
--- a/src/extraction/core/functions/utils.py
+++ /dev/null
@@ -1,19 +0,0 @@
-from collections import deque
-
-
-def depth(d):
- """
- Copied from https://stackoverflow.com/a/23499088
-
- Used to determine the depth of our config trees and fill them
- """
- queue = deque([(id(d), d, 1)])
- memo = set()
- while queue:
- id_, o, level = queue.popleft()
- if id_ in memo:
- continue
- memo.add(id_)
- if isinstance(o, dict):
- queue += ((id(v), v, level + 1) for v in o.values())
- return level
diff --git a/src/extraction/core/functions/versioning.py b/src/extraction/core/functions/versioning.py
index 40c1e8f9c96d0e17af4fb88c8452da9c70fbe065..4c77a9f8ebdc4ce65626c6ac433ddd85194f01eb 100644
--- a/src/extraction/core/functions/versioning.py
+++ b/src/extraction/core/functions/versioning.py
@@ -2,6 +2,7 @@ import git
def get_sha():
+ # FIXME Unused, but *should* be used...
repo = git.Repo(search_parent_directories=True)
sha = repo.head.object.hexsha
return sha
diff --git a/src/extraction/core/omero.py b/src/extraction/core/omero.py
deleted file mode 100644
index e69aa7473a493dff46bb720dd4e8387d0a46bca0..0000000000000000000000000000000000000000
--- a/src/extraction/core/omero.py
+++ /dev/null
@@ -1,35 +0,0 @@
-from omero.gateway import BlitzGateway
-from tqdm import tqdm
-
-
-# Helper funs
-def connect_omero():
- conn = BlitzGateway(*get_creds(), host="islay.bio.ed.ac.uk", port=4064)
- conn.connect()
- return conn
-
-
-def get_creds():
- return (
- "upload",
- "***REMOVED***", # OMERO Password
- )
-
-
-def download_file(f):
- """
- Download file in chunks using FileWrapper object
- """
- desc = (
- "Downloading "
- + f.getFileName()
- + " ("
- + str(round(f.getFileSize() / 1000**2, 2))
- + "Mb)"
- )
-
- down_file = bytearray()
- for c in tqdm(f.getFileInChunks(), desc=desc):
- down_file += c
-
- return down_file
diff --git a/src/logfile_parser/swainlab_parser.py b/src/logfile_parser/swainlab_parser.py
index a83cb87c43ee861aa3cb14df3790c08c1a7c0cf3..cb39750adcba447cc45b01cdf17b5d3a71bfd326 100644
--- a/src/logfile_parser/swainlab_parser.py
+++ b/src/logfile_parser/swainlab_parser.py
@@ -1,5 +1,5 @@
#!/usr/bin/env jupyter
-
+# TODO should this be merged to the regular logfile_parser structure?
"""
Description of new logfile:
@@ -32,6 +32,7 @@ New grammar
- Tables are assumed to end with an empty line.
"""
+import logging
import typing as t
from pathlib import PosixPath
@@ -229,7 +230,9 @@ def parse_from_grammar(filepath: str, grammar: t.Dict):
subkey = "_".join((key, subkey))
d[subkey] = parse_x(header, **subvalues)
except Exception as e:
- print(f"Parsing failed for key {key}")
+ logging.getLogger("aliby").critical(
+ f"Parsing failed for key {key}"
+ )
raise (e)
return d
@@ -298,18 +301,3 @@ def parse_x(string: str, type: str, **kwargs):
def parse_from_swainlab_grammar(filepath: t.Union[str, PosixPath]):
return parse_from_grammar(filepath, grammar)
-
-
-# test_file = "/home/alan/Documents/dev/skeletons/scripts/dev/C1_60x.log"
-# test_file = "/home/alan/Documents/dev/skeletons/scripts/dev/bak"
-# test_file = "/home/alan/Documents/dev/skeletons/scripts/dev/two_tables.log"
-# test_file = "/home/alan/Downloads/pH_med_to_low 1.log"
-# test_file = "/home/alan/Documents/dev/skeletons/scripts/data/577_2022_12_20_pHCalibrate6_7_00/pHCalibrate6_7.log"
-
-
-# d = parse_from_grammar(test_file, grammar)
-# print(d)
-
-# from logfile_parser.legacy import get_legacy_log_example_interface
-
-# original = get_legacy_log_example_interface()
diff --git a/src/postprocessor/benchmarks/post_processing.py b/src/postprocessor/benchmarks/post_processing.py
index d17651a8dc839d19abba7f38528437dbc22c3a03..d36e2afbae5c0e40e219dd322e8e4acf81a9836b 100644
--- a/src/postprocessor/benchmarks/post_processing.py
+++ b/src/postprocessor/benchmarks/post_processing.py
@@ -1,3 +1,4 @@
+# TODO remove/to snippets?
"""
Post-processing utilities
diff --git a/src/postprocessor/chainer.py b/src/postprocessor/chainer.py
index a5387ccb1138357e54c46fb3e2cf6abadc5b6424..4b101922ee882744ee7d1e311a3b35bbe7f9bdec 100644
--- a/src/postprocessor/chainer.py
+++ b/src/postprocessor/chainer.py
@@ -4,7 +4,6 @@ import re
import typing as t
from copy import copy
-import numpy as np
import pandas as pd
from agora.io.signal import Signal
@@ -17,61 +16,36 @@ from postprocessor.core.lineageprocess import LineageProcessParameters
class Chainer(Signal):
"""
Extend Signal by applying post-processes and allowing composite signals that combine basic signals.
+ It "chains" multiple processes upon fetching a dataset to produce the desired datasets.
Instead of reading processes previously applied, it executes
them when called.
"""
- # these no longer seem to be used
- #process_types = ("multisignal", "processes", "reshapers")
- #common_chains = {}
+ _synonyms = {
+ "m5m": ("extraction/GFP/max/max5px", "extraction/GFP/max/median")
+ }
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
- for channel in self.candidate_channels:
- # find first channel in h5 file that corresponds to a candidate_channel
- # but channel is redefined. why is there a loop over candidate channels?
- # what about capitals?
- try:
- channel = [
- ch for ch in self.channels if re.match("channel", ch)
- ][0]
- break
- except:
- # is this still a good idea?
- pass
- try:
- # what's this?
- # composite statistic comprising the quotient of two others
- equivalences = {
- "m5m": (
- f"extraction/{channel}/max/max5px",
- f"extraction/{channel}/max/median",
- ),
- }
+ def replace_path(path: str, bgsub: bool = ""):
# function to add bgsub to paths
- def replace_path(path: str, bgsub: str = ""):
- channel = path.split("/")[1]
- if "bgsub" in bgsub:
- # add bgsub to path
- path = re.sub(channel, f"{channel}_bgsub", path)
- return path
-
- # for composite statistics
- # add chain with and without bgsub
- self.common_chains = {
- alias
- + bgsub: lambda **kwargs: self.get(
- replace_url(denominator, alias + bgsub), **kwargs
- )
- / self.get(replace_path(numerator, alias + bgsub), **kwargs)
- for alias, (denominator, numerator) in equivalences.items()
- for bgsub in ("", "_bgsub")
- }
- except:
- # Is this still a good idea?
- pass
+ channel = path.split("/")[1]
+ suffix = "_bgsub" if bgsub else ""
+ path = re.sub(channel, f"{channel}{suffix}", path)
+ return path
+
+ # Add chain with and without bgsub for composite statistics
+ self.common_chains = {
+ alias
+ + bgsub: lambda **kwargs: self.get(
+ replace_path(denominator, alias + bgsub), **kwargs
+ )
+ / self.get(replace_path(numerator, alias + bgsub), **kwargs)
+ for alias, (denominator, numerator) in self.synonyms.items()
+ for bgsub in ("", "_bgsub")
+ }
def get(
self,
@@ -83,7 +57,6 @@ class Chainer(Signal):
**kwargs,
):
"""Load data from an h5 file."""
- 1/0
if dataset in self.common_chains:
# get dataset for composite chains
data = self.common_chains[dataset](**kwargs)
@@ -96,7 +69,7 @@ class Chainer(Signal):
# keep data only from early time points
data = self.get_retained(data, retain)
# data = data.loc[data.notna().sum(axis=1) > data.shape[1] * retain]
- if (stages and "stage" not in data.columns.names):
+ if stages and "stage" not in data.columns.names:
# return stages as additional column level
stages_index = [
x
@@ -149,7 +122,5 @@ class Chainer(Signal):
if process_type == "reshapers":
if process == "merger":
raise (NotImplementedError)
- merges = process.as_function(result, **params)
- result = self.apply_merges(result, merges)
self._intermediate_steps.append(result)
return result
diff --git a/src/postprocessor/core/abc.py b/src/postprocessor/core/abc.py
index a299e19f6f87a33e0d9a6c5c099b291f377811dc..ffb677b4d7e5ad7be5249e23df552b7355bdf1f6 100644
--- a/src/postprocessor/core/abc.py
+++ b/src/postprocessor/core/abc.py
@@ -16,6 +16,7 @@ class PostProcessABC(ProcessABC):
@classmethod
def as_function(cls, data, *extra_data, **kwargs):
+ # FIXME can this be a __call__ method instead?
# Find the parameter's default
parameters = cls.default_parameters(**kwargs)
return cls(parameters=parameters).run(data, *extra_data)
diff --git a/src/postprocessor/core/group.py b/src/postprocessor/core/group.py
deleted file mode 100644
index 5613e5de2a0fd4981f98e7771bab47eda9307802..0000000000000000000000000000000000000000
--- a/src/postprocessor/core/group.py
+++ /dev/null
@@ -1,137 +0,0 @@
-"""
-Class to group multiple positions into one using one different available criteria.
-"""
-
-import re
-from pathlib import Path
-
-import h5py
-import pandas as pd
-from agora.io.bridge import groupsort
-from agora.io.signal import Signal
-
-from postprocessor.core.abc import ParametersABC, ProcessABC
-
-
-class GroupParameters(ParametersABC):
- def __init__(self, by="name", processes=[], signals=[]):
- self.by = by
- self.signals = signals
- self.processes = processes
-
- @classmethod
- def default(cls):
- return cls.from_dict({"by": "name", "signals": [], "processes": []})
-
-
-class Group(ProcessABC):
- def __init__(self, parameters):
- super().__init__(parameters)
-
- def get_position_filenames(self, exp_root, poses):
- """
- Get filenames as a dictionary where the key is the position and value the filename.
- """
- central_store = Path(exp_root) / "store.h5"
- if central_store.exists():
- hdf = h5py.File(central_store, "r")
- self.filenames = [
- pos.attrs["filename"] for pos in hdf["/positions/"]
- ]
- hdf.close()
- else: # If no central store just list position files in expt root folder
- fullfiles = [x for x in Path(exp_root).glob("*store.h5")]
- files = [x.name for x in Path(exp_root).glob("*store.h5")]
- filenames = [False for _ in poses]
- for i, pos in enumerate(poses):
- matches = [
- True if re.match(pos + ".*.h5", fname) else False
- for fname in files
- ]
- if any(matches):
- assert sum(matches) == 1, "More than one match"
- filenames[i] = (pos, fullfiles[matches.index(True)])
-
- self.filenames = {
- fname[0]: fname[1] for fname in filenames if fname
- }
-
- self.positions = list(self.filenames.keys())
- return self.filenames
-
- def get_signals(self):
- # hdf = h5py.File(central_store, "r")
- # keys_d = groupsort(keys)
- self.signals = {pos: {} for pos in self.filenames.keys()}
- for pos, fname in self.filenames.items():
- for signal in self.parameters.signals:
- self.signals[pos][signal] = pd.read_hdf(fname, signal)
-
- return self.signals
-
- def gen_groups(self):
- if self.by == "group": # Use group names in metadata
- pass
- elif self.by == "name": # Infer groups from signal concatenation
- # Remove last four characters and find commonalities larger than 4
- # characters between posnames and group them that way.
- groupnames = list(set([x[:-3] for x in self.positions]))
- self.group_signal_tree = {group: [] for group in groupnames}
- self.poses_grouped = {group: [] for group in groupnames}
- for pos in self.positions:
- group = groupnames[groupnames.index(pos[:-3])]
- self.group_signal_tree[group].append(self.signals[pos])
- self.poses_grouped[group].append(pos)
-
- elif (
- type(self.by) == tuple
- ): # Manually give groups as tuple or list of positions
- pass
-
- def concat_signals(self):
- self.concated_signals = {group: {} for group in self.group_signal_tree}
- for k, group in self.group_signal_tree.items():
- for signal in self.parameters.signals:
- self.concated_signals[k][signal] = pd.concat(
- [g[signal] for g in group], keys=self.poses_grouped[k]
- )
-
- return self.concated_signals
-
- def process_signals(self, grouped_signals):
- pass
-
- def run(self, central_store, poses):
-
- self.get_position_filenames(central_store, poses)
- self.get_signals()
- self.gen_groups()
- self.concat_signals()
- # processed_signals = self.process_signals(grouped_signals)
-
- return self.concated_signals
- # return processed_signals
-
-
-poses = [
- x.name.split("store")[0]
- for x in Path(
- "/shared_libs/pipeline-core/scripts/data/ph_calibration_dual_phl_ura8_5_04_5_83_7_69_7_13_6_59__01"
- ).rglob("*")
- if x.name != "images.h5"
-]
-gr = Group(
- GroupParameters(
- signals=[
- "/extraction/general/None/area",
- "/extraction/mCherry/np_max/median",
- ]
- )
-)
-gr.run(
- central_store="/shared_libs/pipeline-core/scripts/data/ph_calibration_dual_phl_ura8_5_04_5_83_7_69_7_13_6_59__01",
- poses=poses,
-)
-signal = Signal(
- "/shared_libs/pipeline-core/scripts/data/ph_calibration_dual_phl_ura8_5_04_5_83_7_69_7_13_6_59__01/ph_5_04_001store.h5"
-)
diff --git a/src/postprocessor/core/lineageprocess.py b/src/postprocessor/core/lineageprocess.py
index 34ccf2349430c839ab1616c1d779ca0ff995dd27..f10d5b3e3fef47e23724115ab0f99669a5b6ad94 100644
--- a/src/postprocessor/core/lineageprocess.py
+++ b/src/postprocessor/core/lineageprocess.py
@@ -1,3 +1,4 @@
+# TODO Module docstring
import typing as t
from abc import abstractmethod
@@ -7,8 +8,6 @@ import pandas as pd
from agora.abc import ParametersABC
from postprocessor.core.abc import PostProcessABC
-# from agora.utils.lineage import group_matrix
-
class LineageProcessParameters(ParametersABC):
"""
@@ -47,18 +46,14 @@ class LineageProcess(PostProcessABC):
Overrides PostProcess.as_function classmethod.
Lineage functions require lineage information to be passed if run as function.
"""
- # if isinstance(lineage, np.ndarray):
- # lineage = group_matrix(lineage, n_keys=2)
-
parameters = cls.default_parameters(**kwargs)
return cls(parameters=parameters).run(
data, lineage=lineage, *extra_data
)
- # super().as_function(data, *extra_data, lineage=lineage, **kwargs)
def load_lineage(self, lineage):
"""
Reshape the lineage information if needed
"""
-
+ # TODO does this need to be a function?
self.lineage = lineage
diff --git a/src/postprocessor/core/processor.py b/src/postprocessor/core/processor.py
index 37efd0d4128d7ce6744fd0397cd0748ccfe30c52..4b9161c0fe7b8a1a70b60a162fb4d575dba3c494 100644
--- a/src/postprocessor/core/processor.py
+++ b/src/postprocessor/core/processor.py
@@ -56,26 +56,16 @@ class PostProcessorParameters(ParametersABC):
"/extraction/general/None/volume",
],
],
- # [
- # "savgol",
- # [
- # "/extraction/general/None/volume",
- # ],
- # ],
[
"dsignal",
[
"/extraction/general/None/volume",
- # "/postprocessing/savgol/extraction_general_None_volume",
],
],
[
"bud_metric",
[
"/extraction/general/None/volume",
- # "/postprocessing/dsignal/postprocessing_savgol_extraction_general_None_volume",
- # "/postprocessing/savgol/extraction_general_None_volume",
- # "/postprocessing/dsignal/extraction_general_None_volume",
],
],
[
@@ -84,15 +74,6 @@ class PostProcessorParameters(ParametersABC):
"/postprocessing/bud_metric/extraction_general_None_volume",
],
],
- # [
- # "aggregate",
- # [
- # [
- # "/extraction/general/None/volume",
- # "postprocessing/dsignal/extraction_general_None_volume",
- # ],
- # ],
- # ],
],
}
param_sets = {
@@ -105,22 +86,12 @@ class PostProcessorParameters(ParametersABC):
outpaths["aggregate"] = "/postprocessing/experiment_wide/aggregated/"
if "ph_batman" in kind:
- # targets["processes"]["bud_metric"].append(
- # [
- # [
- # "/extraction/em_ratio/np_max/mean",
- # "/extraction/em_ratio/np_max/median",
- # ],
- # ]
- # )
targets["processes"]["dsignal"].append(
[
"/extraction/em_ratio/np_max/mean",
"/extraction/em_ratio/np_max/median",
"/extraction/em_ratio_bgsub/np_max/mean",
"/extraction/em_ratio_bgsub/np_max/median",
- # "/postprocessing/bud_metric/extraction_em_ratio_np_max_mean",
- # "/postprocessing/bud_metric/extraction_em_ratio_np_max_median",
]
)
targets["processes"]["aggregate"].append(
@@ -132,10 +103,6 @@ class PostProcessorParameters(ParametersABC):
"/extraction/em_ratio_bgsub/np_max/median",
"/extraction/gsum/np_max/median",
"/extraction/gsum/np_max/mean",
- # "postprocessing/bud_metric/extraction_em_ratio_np_max_mean",
- # "postprocessing/bud_metric/extraction_em_ratio_np_max_median",
- # "postprocessing/dsignal/postprocessing_bud_metric_extraction_em_ratio_np_max_median",
- # "postprocessing/dsignal/postprocessing_bud_metric_extraction_em_ratio_np_max_mean",
]
],
)
@@ -178,6 +145,7 @@ class PostProcessor(ProcessABC):
self.targets = parameters["targets"]
def run_prepost(self):
+ # TODO Split function
"""Important processes run before normal post-processing ones"""
merge_events = self.merger.run(
@@ -301,12 +269,7 @@ class PostProcessor(ProcessABC):
return x
def run(self):
- # import cProfile
- # import pstats
-
- # profile = cProfile.Profile()
- # profile.enable()
-
+ # TODO Documentation :) + Split
self.run_prepost()
for process, datasets in tqdm(self.targets["processes"]):
diff --git a/src/postprocessor/core/reshapers/bud_metric.py b/src/postprocessor/core/reshapers/bud_metric.py
index 0893db667c5250fe20e2447d17712eb920d98e8d..c527134dc95ac69ef3e6b6d54e0b2e67b32fc41a 100644
--- a/src/postprocessor/core/reshapers/bud_metric.py
+++ b/src/postprocessor/core/reshapers/bud_metric.py
@@ -33,7 +33,6 @@ class bud_metric(LineageProcess):
mother_bud_ids: Dict[pd.Index, Tuple[pd.Index]] = None,
):
if mother_bud_ids is None:
- # filtered_lineage = self.filter_signal_cells(signal)
mother_bud_ids = mb_array_to_dict(self.lineage)
return self.get_bud_metric(signal, mother_bud_ids)
diff --git a/src/postprocessor/core/reshapers/buddings.py b/src/postprocessor/core/reshapers/buddings.py
index 0b01dad70c5eab067f26b430e5451fdf70bcb4e2..4b3fbba9db0a2a8b61c69903afe78664c1dcddca 100644
--- a/src/postprocessor/core/reshapers/buddings.py
+++ b/src/postprocessor/core/reshapers/buddings.py
@@ -25,16 +25,17 @@ class buddingsParameters(LineageProcessParameters):
FIXME: Add docs.
"""
-
_defaults = {"lineage_location": "postprocessing/lineage_merged"}
+# TODO Why not capitalized?
class buddings(LineageProcess):
"""
Calculate buddings in a trap assuming one mother per trap
returns a pandas series with the buddings
"""
+ # TODO might want to define "buddings" more scientifically
def __init__(self, parameters: buddingsParameters):
super().__init__(parameters)
diff --git a/src/postprocessor/core/reshapers/merger.py b/src/postprocessor/core/reshapers/merger.py
index 16d2b598ea720027a2d0ca566adfc1b73cc5c31a..1fbf4155fd04fe41b7a76e6aeb0253a6746bf2b2 100644
--- a/src/postprocessor/core/reshapers/merger.py
+++ b/src/postprocessor/core/reshapers/merger.py
@@ -3,7 +3,7 @@ from agora.abc import ParametersABC
from postprocessor.core.abc import PostProcessABC
from postprocessor.core.functions.tracks import get_joinable
-
+# TODO Why not capitalized?
class mergerParameters(ParametersABC):
"""
:param tol: float or int threshold of average (prediction error/std) necessary
@@ -22,6 +22,7 @@ class mergerParameters(ParametersABC):
}
+# TODO Why not capitalized?
class merger(PostProcessABC):
"""
Combines rows of tracklet that are likely to be the same.
@@ -34,8 +35,4 @@ class merger(PostProcessABC):
joinable = []
if signal.shape[1] > 4:
joinable = get_joinable(signal, tol=self.parameters.tolerance)
- # merged, _ = merge_tracks(signal) # , min_len=self.window + 1)
- # indices = (*zip(*merged.index.tolist()),)
- # names = merged.index.names
- # return {name: ids for name, ids in zip(names, indices)}
return joinable
diff --git a/src/postprocessor/core/reshapers/picker.py b/src/postprocessor/core/reshapers/picker.py
index bfabc76ad89152f7551a524d375506f1e83fd77d..b64bfcee9f4f795498b28e7920f35b041726e5e1 100644
--- a/src/postprocessor/core/reshapers/picker.py
+++ b/src/postprocessor/core/reshapers/picker.py
@@ -1,34 +1,20 @@
-# from abc import ABC, abstractmethod
-
-# from copy import copy
-# from itertools import groupby
-# from typing import List, Tuple, Union
import typing as t
-from typing import Union
-# import igraph as ig
import numpy as np
import pandas as pd
from agora.abc import ParametersABC
from agora.io.cells import Cells
-# from postprocessor.core.functions.tracks import max_nonstop_ntps, max_ntps
from agora.utils.association import validate_association
from postprocessor.core.lineageprocess import LineageProcess
-# from utils_find_1st import cmp_equal, find_1st
-
class pickerParameters(ParametersABC):
_defaults = {
"sequence": [
["lineage", "intersection", "families"],
- # ["condition", "intersection", "any_present", 0.7],
- # ["condition", "intersection", "growing", 80],
["condition", "intersection", "present", 7],
- # ["condition", "intersection", "mb_guess", 3, 0.7],
- # ("lineage", "intersection", "full_families"),
],
}
@@ -80,16 +66,8 @@ class picker(LineageProcess):
idx = idx[valid_indices]
mothers_daughters = mothers_daughters[valid_lineage]
-
- # return mothers_daughters, idx
return idx
- def loc_lineage(self, kymo: pd.DataFrame, how: str, lineage=None):
- _, valid_indices = self.pick_by_lineage(
- kymo, how, mothers_daughters=lineage
- )
- return kymo.loc[[tuple(x) for x in valid_indices]]
-
def pick_by_condition(self, signals, condition, thresh):
idx = self.switch_case(signals, condition, thresh)
return idx
@@ -122,7 +100,7 @@ class picker(LineageProcess):
indices = indices.intersection(new_indices)
else:
- print("WARNING:Picker: No lineage assignment")
+ self._log(f"No lineage assignment")
indices = np.array([])
return np.array(list(indices))
@@ -131,7 +109,7 @@ class picker(LineageProcess):
self,
signals: pd.DataFrame,
condition: str,
- threshold: Union[float, int, list],
+ threshold: t.Union[float, int, list],
):
if len(threshold) == 1:
threshold = [_as_int(*threshold, signals.shape[1])]
@@ -145,7 +123,7 @@ class picker(LineageProcess):
return set(signals.index[case_mgr[condition](signals, *threshold)])
-def _as_int(threshold: Union[float, int], ntps: int):
+def _as_int(threshold: t.Union[float, int], ntps: int):
if type(threshold) is float:
threshold = ntps * threshold
return threshold
diff --git a/src/postprocessor/grouper.py b/src/postprocessor/grouper.py
index 0bcfd8418ea92e1daa417a18093d754c32388f55..09bff52e5c22b54f5b6346a17a3d7e1a6cafc0af 100644
--- a/src/postprocessor/grouper.py
+++ b/src/postprocessor/grouper.py
@@ -15,11 +15,6 @@ import pandas as pd
import seaborn as sns
from pathos.multiprocessing import Pool
-from agora.utils.kymograph import (
- drop_level,
- get_mother_ilocs_from_daughters,
- intersection_matrix,
-)
from postprocessor.chainer import Chainer
diff --git a/tests/agora/example_test.py b/tests/agora/example_test.py
deleted file mode 100644
index 539b410aefd240a77593a2a0509c3e36117e82ba..0000000000000000000000000000000000000000
--- a/tests/agora/example_test.py
+++ /dev/null
@@ -1,25 +0,0 @@
-"""This is an example test file to show the structure."""
-import pytest
-
-from agora.utils.example import ExampleClass, example_function
-
-
-class TestExampleClass:
- x = ExampleClass(1)
-
- def test_add_one(self):
- assert self.x.add_one() == 2
-
- def test_add_n(self):
- assert self.x.add_n(10) == 11
-
-
-def test_example_function():
- x = example_function(1)
- assert isinstance(x, ExampleClass)
- assert x.parameter == 1
-
-
-def test_example_function_fail():
- with pytest.raises(ValueError):
- example_function("hello")
diff --git a/tests/aliby/network/test_baby_client.py b/tests/aliby/network/test_baby_client.py
deleted file mode 100644
index a8d131a04d3a06ce03f4634252c4e497111f986f..0000000000000000000000000000000000000000
--- a/tests/aliby/network/test_baby_client.py
+++ /dev/null
@@ -1,91 +0,0 @@
-import pytest
-
-pytest.mark.skip
-
-import json
-import time
-
-import numpy as np
-
-# from aliby.experiment import ExperimentLocal
-from aliby.baby_client import BabyClient
-from aliby.tile.tiler import Tiler
-
-
-@pytest.mark.skip(
- reason="No longer usable, requires local files. Kept until replaced."
-)
-def test_client():
- root_dir = (
- "/Users/s1893247/PhD/pipeline-core/data/glclvl_0"
- ".1_mig1_msn2_maf1_sfp1_dot6_03"
- )
-
- expt = ExperimentLocal(root_dir, finished=True)
- seg_expt = Tiler(expt, finished=True)
-
- print(seg_expt.positions)
- seg_expt.current_position = "pos007"
-
- config = {
- "camera": "evolve",
- "channel": "brightfield",
- "zoom": "60x",
- "n_stacks": "5z",
- }
-
- baby_client = BabyClient(expt, **config)
-
- print("The session is {}".format(baby_client.sessions["default"]))
-
- # Channel 0, 0, X,Y,Z all
- num_timepoints = 5
-
- traps_tps = [
- seg_expt.get_tiles_timepoint(
- tp, tile_size=81, channels=[0], z=[0, 1, 2, 3, 4]
- ).squeeze()
- for tp in range(num_timepoints)
- ]
-
- segmentations = []
- try:
- for i, timpoint in enumerate(traps_tps):
- print("Sending timepoint {};".format(i))
- status = baby_client.queue_image(
- timpoint,
- baby_client.sessions["default"],
- assign_mothers=True,
- return_baprobs=True,
- with_edgemasks=True,
- )
- while True:
- try:
- print("Loading.", end="")
- result = baby_client.get_segmentation(
- baby_client.sessions["default"]
- )
- except:
- print(".", end="")
- time.sleep(1)
- continue
- break
- print("Received timepoint {}".format(i))
- segmentations.append(result)
- except Exception as e:
- print(segmentations)
- raise e
-
- with open("segmentations.json", "w") as fd:
- json.dump(segmentations, fd)
-
- print("Done.")
- # print(len(segmentations[0]))
- # for i in range(5):
- # print("trap {}".format(i))
- # for k, v in segmentations[0][i].items():
- # print(k, v)
- #
- # import matplotlib.pyplot as plt
- # plt.imshow(np.squeeze(batches[0][0, ..., 0]))
- # plt.savefig('test_baby.pdf')
diff --git a/tests/aliby/network/test_post_processing.py b/tests/aliby/network/test_post_processing.py
index 3ac8ba403602778a92fa482f4c64395df59663b0..dbf9b8a9bb9e0bfc4703d79b93dd1f28d4f68016 100644
--- a/tests/aliby/network/test_post_processing.py
+++ b/tests/aliby/network/test_post_processing.py
@@ -7,14 +7,6 @@ import skimage.morphology as morph
from scipy import ndimage
from skimage import draw
-# from aliby.post_processing import (
-# circle_outline,
-# conical,
-# ellipse_perimeter,
-# union_of_spheres,
-# volume_of_sphere,
-# )
-
@pytest.mark.skip(
reason="No longer usable, post_processing unused inside aliby. Kept temporarily"
diff --git a/tests/aliby/network/test_tiler.py b/tests/aliby/network/test_tiler.py
index 52187b58ee4db91d11856bf7cdd6584adf4c1b6d..1ec1755c11ee519c9bdcd1556ed331dfb2d97e49 100644
--- a/tests/aliby/network/test_tiler.py
+++ b/tests/aliby/network/test_tiler.py
@@ -19,6 +19,12 @@ def define_parser():
return parser
+def initialise_dummy():
+ tiler_parameters = TilerParameters.default().to_dict()
+ dummy_tiler = Tiler.dummy(tiler_parameters)
+ return dummy_tiler
+
+
def initialise_objects(data_path, template=None):
image = ImageLocalOME(data_path)
tiler = Tiler.from_image(image, TilerParameters.default())
@@ -53,6 +59,8 @@ if __name__ == "__main__":
parser = define_parser()
args = parser.parse_args()
+ dummy_tiler = initialise_dummy()
+
tiler = initialise_objects(args.root_dir, template=args.template)
if args.position is not None:
diff --git a/tests/aliby/pipeline/test_image.py b/tests/aliby/pipeline/test_image.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a2114006c0bd2703bbfed159ee99201853bae05
--- /dev/null
+++ b/tests/aliby/pipeline/test_image.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python3
+import numpy as np
+import dask.array as da
+import pytest
+
+from aliby.io.image import ImageDummy
+
+tiler_parameters = {"tile_size": 117, "ref_channel": "Brightfield", "ref_z": 0}
+
+sample_da = da.from_array(np.array([[1, 2], [3, 4]]))
+# Make it 5-dimensional
+sample_da = da.reshape(
+ sample_da, (1, 1, 1, sample_da.shape[-2], sample_da.shape[-1])
+)
+
+
+@pytest.mark.parametrize("sample_da", [sample_da])
+@pytest.mark.parametrize("dim", [2])
+@pytest.mark.parametrize("n_empty_slices", [4])
+@pytest.mark.parametrize("image_position", [1])
+def test_pad_array(sample_da, dim, n_empty_slices, image_position):
+ """Test ImageDummy.pad_array() method"""
+ # create object
+ imgdmy = ImageDummy(tiler_parameters)
+ # pads array
+ padded_da = imgdmy.pad_array(
+ sample_da,
+ dim=dim,
+ n_empty_slices=n_empty_slices,
+ image_position=image_position,
+ )
+
+ # select which dimension to index the multidimensional array
+ indices = {dim: image_position}
+ ix = [
+ indices.get(dim, slice(None))
+ for dim in range(padded_da.compute().ndim)
+ ]
+
+ # Checks that original image array is there and is at the correct index
+ assert np.array_equal(padded_da.compute()[ix], sample_da.compute()[0])
+ # Checks that the additional axis is extended correctly
+ assert padded_da.compute().shape[dim] == n_empty_slices + 1
diff --git a/tests/postprocessor/test_interpolate.py b/tests/postprocessor/test_interpolate.py
new file mode 100644
index 0000000000000000000000000000000000000000..c9c993de2bf0833de6fa053b1993576588a3f7ee
--- /dev/null
+++ b/tests/postprocessor/test_interpolate.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+import numpy as np
+import pandas as pd
+from postprocessor.core.processes.interpolate import (
+ interpolate,
+ interpolateParameters,
+)
+
+
+def dummy_signal_array(n_cells, n_tps):
+ """Creates dummy signal array, i.e. increasing gradient"""
+ signal = np.array([np.linspace(1, 2, n_tps) for _ in range(n_cells)])
+ return signal
+
+
+def test_dummy_signal_array():
+ ds = dummy_signal_array(5, 10)
+ # Check dimensions
+ assert ds.shape[0] == 5
+ assert ds.shape[1] == 10
+
+
+def randomly_add_na(input_array, num_of_na):
+ """Randomly replaces a 2d numpy array with NaNs, number of NaNs specified"""
+ input_array.ravel()[
+ np.random.choice(input_array.size, num_of_na, replace=False)
+ ] = np.nan
+ return input_array
+
+
+def test_interpolate():
+ dummy_array = dummy_signal_array(5, 10)
+ # Poke holes so interpolate can fill
+ holey_array = randomly_add_na(dummy_array, 15)
+
+ dummy_signal = pd.DataFrame(dummy_array)
+ holey_signal = pd.DataFrame(holey_array)
+
+ interpolate_runner = interpolate(interpolateParameters.default())
+ interpolated_signal = interpolate_runner.run(holey_signal)
+
+ subtr = interpolated_signal - dummy_signal
+ # Check that interpolated values are the ones that exist in the dummy
+ assert np.nansum(subtr.to_numpy()) == 0
+ # TODO: Check that if there are NaNs remaining after interpolation, they
+ # are at the ends