Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • swain-lab/aliby/aliby-mirror
  • swain-lab/aliby/alibylite
2 results
Show changes
Commits on Source (1025)
925 additional commits have been omitted to prevent performance issues.
Hide whitespace changes
Inline Side-by-side
Showing
with 390 additions and 1846 deletions
.flake8 0 → 100644
View file @ 95958192
[flake8]
ignore = E203, E266, E501, W503, F403, F401
max-line-length = 79
select = B,C,E,F,W,T4,B9
exclude =
# No need to traverse our git directory
.git,
# There's no value in checking cache directories
__pycache__,
# Ignore virtual environment contents
.venv
# The conf file is mostly autogenerated, ignore it
docs/source/conf.py,
# The old directory contains Flake8 2.0
old,
# This contains our built documentation
build,
# This contains builds of flake8 that we don't want to check
dist,
# Any data produced inside the folder during development
data/,
# Temporarily ignore tests
tests/
max-complexity = 18
.gitignore
View file @ 95958192
......@@ -110,7 +110,6 @@ venv.bak/
omero_py/omeroweb/
omero_py/pipeline/
**.ipynb
data/
notebooks/
*.pdf
*.h5
......
.gitlab-ci.yml
View file @ 95958192
image: python:3.7
image: python:3.8
cache:
key: "project-${CI_JOB_NAME}"
......@@ -12,8 +12,8 @@ variables:
TRIGGER_PYPI_NAME: ""
stages:
- test
- check
- tests
- checks
# - release
before_script:
......@@ -25,28 +25,49 @@ before_script:
# - git remote rm origin && git remote add origin https://${ACCESS_TOKEN_NAME}:${ACCESS_TOKEN}@${CI_SERVER_HOST}/${CI_PROJECT_PATH}.git
# - git config pull.rebase false
# - git pull origin HEAD:master
- rm -rf ~/.cache/pypoetry
- if [ ${var+TRIGGER_PYPI_NAME} ]; then echo "Pipeline triggered by ${TRIGGER_PYPI_NAME}"; poetry add ${TRIGGER_PYPI_NAME}@latest; fi
- poetry install -vv
# - rm -rf ~/.cache/pypoetry
# - if [ ${var+TRIGGER_PYPI_NAME} ]; then echo "Pipeline triggered by ${TRIGGER_PYPI_NAME}"; poetry add ${TRIGGER_PYPI_NAME}@latest; fi
# - export WITHOUT="docs,network";
- export ARGS="--with test,dev";
- if [[ "$CI_STAGE_NAME" == "tests" ]]; then echo "Installing system dependencies for ${CI_STAGE_NAME}"; apt update && apt install -y ffmpeg libsm6 libxext6; fi
- if [[ "$CI_JOB_NAME" == "Static Type" ]]; then echo "Activating development group"; export ARGS="${ARGS},dev"; fi
- if [[ "$CI_JOB_NAME" == "Network Tools Tests" ]]; then echo "Setting flag to compile zeroc-ice"; export ARGS="${ARGS} --all-extras"; fi
- poetry install -vv $ARGS
Unit test:
stage: test
Local Tests:
stage: tests
script:
- apt update && apt install ffmpeg libsm6 libxext6 -y
- poetry run pytest ./tests/
# - poetry install -vv
- 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
Python Code Lint:
stage: check
Network Tools Tests:
stage: tests
script:
- poetry run black .
- poetry run pytest ./tests/aliby/network
- DIRNAME="test_datasets"
- curl https://zenodo.org/record/7513194/files/test_datasets.tar.gz\?download\=1 -o "test_datasets.tar.gz"
- mkdir -p $DIRNAME
- tar xvf test_datasets.tar.gz -C $DIRNAME
- poetry run pytest -s tests/aliby/pipeline --file $DIRNAME/560_2022_11_30_pypipeline_unit_test_reconstituted_00
Static Type:
stage: check
stage: checks
allow_failure: true
script:
- poetry run black .
- poetry run isort .
- poetry run mypy . --exclude 'setup\.py$'
# We can remove the flag once this is resolved https://github.com/pypa/setuptools/issues/2345
# TODO add more tests before activating auto-release
# Bump_release:
# stage: release
# script:
......
.gitlab/issue_templates/bug.md
View file @ 95958192
## Summary
(Summarize the bug encountered concisely)
{Summarize the bug encountered concisely}
I confirm that I have (if relevant):
- [ ] Read the troubleshooting guide: https://gitlab.com/aliby/aliby/-/wikis/Troubleshooting-(basic)
- [ ] Updated aliby and aliby-baby.
- [ ] Tried the unit test.
- [ ] Tried a scaled-down version of my experiment (distributed=0, filter=0, tps=10)
- [ ] Tried re-postprocessing.
## Steps to reproduce
(How one can reproduce the issue - this is very important)
{How one can reproduce the issue - this is very important}
- aliby version: 0.1.{...}, or if development/unreleased version, commit SHA: {...}
- platform(s):
- [ ] Jura
- [ ] Other Linux, please specify distribution and version: {...}
- [ ] MacOS, please specify version: {...}
- [ ] Windows, please specify version: {...}
- experiment ID: {...}
- Any special things you need to know about this experiment: {...}
## What is the current bug behavior?
......@@ -19,7 +35,12 @@
(Paste any relevant logs - please use code blocks (```) to format console output, logs, and code, as
it's very hard to read otherwise.)
```
{PASTE YOUR ERROR MESSAGE HERE!!}
```
## Possible fixes
(If you can, link to the line of code that might be responsible for the problem)
.pre-commit-config.yaml 0 → 100644
View file @ 95958192
exclude: '^examples/'
repos:
- repo: https://github.com/pre-commit/pre-commit-hooks
rev: v4.3.0
hooks:
# - id: check-yaml
- id: end-of-file-fixer
- id: trailing-whitespace
args: [--markdown-linebreak-ext=md]
- id: check-added-large-files
- id: check-docstring-first
- id: name-tests-test
args: [--pytest-test-first]
- repo: https://github.com/psf/black
rev: 22.6.0
hooks:
- id: black
- repo: https://github.com/PyCQA/isort
rev: 5.10.1
hooks:
- id: isort
- repo: https://github.com/PyCQA/flake8
rev: 4.0.1
hooks:
- id: flake8
exclude: ^tests/ # Temporarily exclude tests
.readthedocs.yaml
View file @ 95958192
......@@ -9,7 +9,7 @@ version: 2
build:
os: ubuntu-20.04
tools:
python: "3.7"
python: "3.8"
# Build documentation in the docs/ directory with Sphinx
sphinx:
......
CONTRIBUTING.md
View file @ 95958192
# 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
......
README.md
View file @ 95958192
# ALIBY (Analyser of Live-cell Imaging for Budding Yeast)
[![docs](https://readthedocs.org/projects/aliby/badge/?version=master)](https://aliby.readthedocs.io/en/latest)
[![PyPI version](https://badge.fury.io/py/aliby.svg)](https://badge.fury.io/py/aliby)
[![readthedocs](https://readthedocs.org/projects/aliby/badge/?version=latest)](https://aliby.readthedocs.io/en/latest)
[![pipeline status](https://git.ecdf.ed.ac.uk/swain-lab/aliby/aliby/badges/master/pipeline.svg)](https://git.ecdf.ed.ac.uk/swain-lab/aliby/aliby/-/pipelines)
[![pipeline](https://gitlab.com/aliby/aliby/badges/master/pipeline.svg?key_text=master)](https://gitlab.com/aliby/aliby/-/pipelines)
[![dev pipeline](https://gitlab.com/aliby/aliby/badges/dev/pipeline.svg?key_text=dev)](https://gitlab.com/aliby/aliby/-/commits/dev)
[![coverage](https://gitlab.com/aliby/aliby/badges/dev/coverage.svg)](https://gitlab.com/aliby/aliby/-/commits/dev)
The core classes and methods for the python microfluidics, microscopy, data analysis and reporting.
### Installation
See [INSTALL.md](./INSTALL.md) for installation instructions.
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
### 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.
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.
To stop them, in the same directory, run:
```shell script
docker-compose stop
```
To analyse local data
```bash
pip install aliby
```
Add any of the optional flags `omero` and `utils` (e.g., `pip install aliby[omero, utils]`). `omero` provides tools to connect with an OMERO server and `utils` provides visualisation, user interface and additional deep learning tools.
See our [installation instructions]( https://aliby.readthedocs.io/en/latest/INSTALL.html ) for more details.
### CLI
If installed via poetry, you have access to a Command Line Interface (CLI)
```bash
aliby-run --expt_id EXPT_PATH --distributed 4 --tps None
```
And to run Omero servers, the basic arguments are shown:
```bash
aliby-run --expt_id XXX --host SERVER.ADDRESS --user USER --password PASSWORD
```
The output is a folder with the original logfiles and a set of hdf5 files, one with the results of each multidimensional inside.
For more information, including available options, see the page on [running the analysis pipeline](https://aliby.readthedocs.io/en/latest/PIPELINE.html)
### Raw data access
## Using specific components
### Access raw data
ALIBY's tooling can also be used as an interface to OMERO servers, for example, to fetch a brightfield channel.
```python
from aliby.io.dataset import Dataset
from aliby.io.image import Image
from aliby.io.omero import Dataset, Image
server_info= {
"host": "host_address",
......@@ -59,7 +60,7 @@ with Image(list(image_ids.values())[0], **server_info) as image:
imgs = dimg[tps, image.metadata["channels"].index("Brightfield"), 2, ...].compute()
# tps timepoints, Brightfield channel, z=2, all x,y
```
### Tiling the raw data
A `Tiler` object performs trap registration. It may be built in different ways but the simplest one is using an image and a the default parameters set.
......@@ -72,43 +73,40 @@ with Image(list(image_ids.values())[0], **server_info) as image:
```
The initialisation should take a few seconds, as it needs to align the images
in time.
in time.
It fetches the metadata from the Image object, and uses the TilerParameters values (all Processes in aliby depend on an associated Parameters class, which is in essence a dictionary turned into a class.)
#### Get a timelapse for a given trap
#### Get a timelapse for a given tile (remote connection)
```python
fpath = "h5/location"
trap_id = 9
trange = list(range(0, 30))
tile_id = 9
trange = range(0, 10)
ncols = 8
riv = remoteImageViewer(fpath)
trap_tps = riv.get_trap_timepoints(trap_id, trange, ncols)
```
This can take several seconds at the moment.
For a speed-up: take fewer z-positions if you can.
trap_tps = [riv.tiler.get_tiles_timepoint(tile_id, t) for t in trange]
If you're not sure what indices to use:
```python
seg_expt.channels # Get a list of channels
channel = 'Brightfield'
ch_id = seg_expt.get_channel_index(channel)
# You can also access labelled traps
m_ts = riv.get_labelled_trap(tile_id=0, tps=[0])
n_traps = seg_expt.n_traps # Get the number of traps
# And plot them directly
riv.plot_labelled_trap(trap_id=0, channels=[0, 1, 2, 3], trange=range(10))
```
#### Get the traps for a given time point
Depending on the network speed can take several seconds at the moment.
For a speed-up: take fewer z-positions if you can.
#### Get the tiles for a given time point
Alternatively, if you want to get all the traps at a given timepoint:
```python
timepoint = 0
seg_expt.get_traps_timepoints(timepoint, tile_size=96, channels=None,
timepoint = (4,6)
tiler.get_tiles_timepoint(timepoint, channels=None,
z=[0,1,2,3,4])
```
### Contributing
See [CONTRIBUTING.md](./CONTRIBUTING.md) for installation instructions.
See [CONTRIBUTING](https://aliby.readthedocs.io/en/latest/INSTALL.html) on how to help out or get involved.
aliby/io/dataset.py deleted 100644 → 0
View file @ e8e2a3bc
#!/usr/bin/env python3
import shutil
from typing import Union
from pathlib import Path, PosixPath
import omero
from aliby.io.omero import Argo
from aliby.io.image import ImageLocal
class DatasetLocal:
"""Load a dataset from a folder
We use a given image of a dataset to obtain the metadata, for we cannot expect folders to contain it straight away.
"""
def __init__(self, dpath: Union[str, PosixPath], *args, **kwargs):
self.fpath = Path(dpath)
assert len(self.get_images()), "No tif files found"
def __enter__(self):
return self
def __exit__(self, *exc):
return False
@property
def dataset(self):
return self.fpath
@property
def name(self):
return self.fpath.name
@property
def unique_name(self):
return self.fpath.name
@property
def date(self):
return ImageLocal(list(self.get_images().values())[0]).date
def get_images(self):
return {f.name: str(f) for f in self.fpath.glob("*.tif")}
@property
def files(self):
if not hasattr(self, "_files"):
self._files = {
f: f for f in self.fpath.rglob("*") if str(f).endswith(".txt")
}
return self._files
def cache_logs(self, root_dir):
for name, annotation in self.files.items():
shutil.copy(annotation, root_dir / name.name)
return True
class Dataset(Argo):
def __init__(self, expt_id, **server_info):
super().__init__(**server_info)
self.expt_id = expt_id
self._files = None
@property
def dataset(self):
return self.conn.getObject("Dataset", self.expt_id)
@property
def name(self):
return self.dataset.getName()
@property
def date(self):
return self.dataset.getDate()
@property
def unique_name(self):
return "_".join(
(
str(self.expt_id),
self.date.strftime("%Y_%m_%d").replace("/", "_"),
self.name,
)
)
def get_images(self):
return {im.getName(): im.getId() for im in self.dataset.listChildren()}
@property
def files(self):
if self._files is None:
self._files = {
x.getFileName(): x
for x in self.dataset.listAnnotations()
if isinstance(x, omero.gateway.FileAnnotationWrapper)
}
if not len(self._files):
raise Exception("exception:metadata: experiment has no annotation files.")
return self._files
@property
def tags(self):
if self._tags is None:
self._tags = {
x.getname(): x
for x in self.dataset.listAnnotations()
if isinstance(x, omero.gateway.TagAnnotationWrapper)
}
return self._tags
def cache_logs(self, root_dir):
for name, annotation in self.files.items():
filepath = root_dir / annotation.getFileName().replace("/", "_")
if str(filepath).endswith("txt") 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
aliby/io/image.py deleted 100644 → 0
View file @ e8e2a3bc
#!/usr/bin/env python3
import typing as t
from datetime import datetime
from pathlib import Path, PosixPath
import dask.array as da
import xmltodict
from dask.array.image import imread
from tifffile import TiffFile
from aliby.io.omero import Argo, get_data_lazy
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):
instatiator = Image
elif isinstance(source, dict) or (
isinstance(source, (str, PosixPath)) and Path(source).is_dir()
):
instatiator = ImageDirectory
elif isinstance(source, str) and Path(source).is_file():
instatiator = ImageLocal
else:
raise ("Invalid data source at {}".format(source))
return instatiator
class ImageLocal:
def __init__(self, path: str, dimorder=None, *args, **kwargs):
self.path = path
self.image_id = str(path)
meta = dict()
try:
with TiffFile(path) as f:
self.meta = xmltodict.parse(f.ome_metadata)["OME"]
for dim in self.dimorder:
meta["size_" + dim.lower()] = int(
self.meta["Image"]["Pixels"]["@Size" + dim]
)
meta["channels"] = [
x["@Name"] for x in self.meta["Image"]["Pixels"]["Channel"]
]
meta["name"] = self.meta["Image"]["@Name"]
meta["type"] = self.meta["Image"]["Pixels"]["@Type"]
except Exception as e:
print("Metadata not found: {}".format(e))
assert "dims" != None, "No dimensional info provided."
# Mark non-existent dimensions for padding
base = "TCZXY"
self.base = base
self.ids = [base.index(i) for i in dimorder]
self._dimorder = dimorder
self._meta = meta
def __enter__(self):
return self
def __exit__(self, *exc):
return False
@property
def name(self):
return self._meta["name"]
@property
def data(self):
return self.get_data_lazy_local()
@property
def date(self):
date_str = [
x
for x in self.meta["StructuredAnnotations"]["TagAnnotation"]
if x["Description"] == "Date"
][0]["Value"]
return datetime.strptime(date_str, "%d-%b-%Y")
@property
def dimorder(self):
"""Order of dimensions in image"""
if not hasattr(self, "_dimorder"):
self._dimorder = self.meta["Image"]["Pixels"]["@DimensionOrder"]
return self._dimorder
@dimorder.setter
def dimorder(self, order: str):
self._dimorder = order
return self._dimorder
@property
def metadata(self):
return self._meta
def get_data_lazy_local(self) -> da.Array:
"""Return 5D dask array. For lazy-loading local multidimensional tiff files"""
if not hasattr(self, "formatted_img"):
if not hasattr(self, "ids"): # Standard dimension order
img = (imread(str(self.path))[0],)
else: # Custom dimension order, we rearrange the axes for compatibility
img = imread(str(self.path))[0]
for i, d in enumerate(self._dimorder):
self._meta["size_" + d.lower()] = img.shape[i]
target_order = (
*self.ids,
*[
i
for i, d in enumerate(self.base)
if d not in self.dimorder
],
)
reshaped = da.reshape(
img,
shape=(
*img.shape,
*[1 for _ in range(5 - len(self.dimorder))],
),
)
img = da.moveaxis(
reshaped, range(len(reshaped.shape)), target_order
)
self._formatted_img = da.rechunk(
img,
chunks=(1, 1, 1, self._meta["size_y"], self._meta["size_x"]),
)
return self._formatted_img
class ImageDirectory(ImageLocal):
"""
Image class for case where 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.
Assumptions:
- Assumes individual folders for individual channels. If only one path provided it assumes it to be brightfield.
- Assumes that images are flat.
- Provides Dimorder as TCZYX
"""
def __init__(self, path: t.Union[str, t.Dict[str, str]], *args, **kwargs):
if isinstance(path, str):
path = {"Brightfield": path}
self.path = path
self.image_id = str(path)
self._meta = dict(channels=path.keys(), name=list(path.values())[0])
# Parse name if necessary
# Build lazy-loading array using dask?
def get_data_lazy_local(self) -> da.Array:
"""Return 5D dask array. For lazy-loading local multidimensional tiff files"""
img = da.stack([imread(v) for v in self.path.values()])
if (
img.ndim < 5
): # Files do not include z-stack: Add and swap with time dimension.
img = da.stack((img,)).swapaxes(0, 2)
# TODO check whether x and y swap is necessary
# Use images to redefine axes
for i, dim in enumerate(("t", "c", "z", "y", "x")):
self._meta["size_" + dim] = img.shape[i]
self._formatted_img = da.rechunk(
img,
chunks=(1, 1, 1, self._meta["size_y"], self._meta["size_x"]),
)
return self._formatted_img
class Image(Argo):
"""
Loads images from OMERO and gives access to the data and metadata.
"""
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__(**server_info)
self.image_id = image_id
# images from OMERO
self._image_wrap = None
@property
def image_wrap(self):
"""
Get images from OMERO
"""
if self._image_wrap is None:
# get images using OMERO
self._image_wrap = self.conn.getObject("Image", self.image_id)
return self._image_wrap
@property
def name(self):
return self.image_wrap.getName()
@property
def data(self):
return get_data_lazy(self.image_wrap)
@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.image_wrap.getSizeX()
meta["size_y"] = self.image_wrap.getSizeY()
meta["size_z"] = self.image_wrap.getSizeZ()
meta["size_c"] = self.image_wrap.getSizeC()
meta["size_t"] = self.image_wrap.getSizeT()
meta["channels"] = self.image_wrap.getChannelLabels()
meta["name"] = self.image_wrap.getName()
return meta
aliby/io/matlab.py deleted 100644 → 0
View file @ e8e2a3bc
"""Read and convert MATLAB files from Swain Lab platform.
TODO: Information that I need from lab members esp J and A
* Lots of examples to try
* Any ideas on what these Map objects are?
TODO: Update Swain Lab wiki
All credit to Matt Bauman for
the reverse engineering at https://nbviewer.jupyter.org/gist/mbauman/9121961
"""
import re
import struct
import sys
from collections.abc import Iterable
from io import BytesIO
import h5py
import numpy as np
import pandas as pd
import scipy
from numpy.compat import asstr
# TODO only use this if scipy>=1.6 or so
from scipy.io import matlab
from scipy.io.matlab.mio5 import MatFile5Reader
from scipy.io.matlab.mio5_params import mat_struct
from aliby.io.utils import read_int, read_string, read_delim
def read_minimat_vars(rdr):
rdr.initialize_read()
mdict = {"__globals__": []}
i = 0
while not rdr.end_of_stream():
hdr, next_position = rdr.read_var_header()
name = asstr(hdr.name)
if name == "":
name = "var_%d" % i
i += 1
res = rdr.read_var_array(hdr, process=False)
rdr.mat_stream.seek(next_position)
mdict[name] = res
if hdr.is_global:
mdict["__globals__"].append(name)
return mdict
def read_workspace_vars(fname):
fp = open(fname, "rb")
rdr = MatFile5Reader(fp, struct_as_record=True, squeeze_me=True)
vars = rdr.get_variables()
fws = vars["__function_workspace__"]
ws_bs = BytesIO(fws.tostring())
ws_bs.seek(2)
rdr.mat_stream = ws_bs
# Guess byte order.
mi = rdr.mat_stream.read(2)
rdr.byte_order = mi == b"IM" and "<" or ">"
rdr.mat_stream.read(4) # presumably byte padding
mdict = read_minimat_vars(rdr)
fp.close()
return mdict
class matObject:
"""A python read-out of MATLAB objects
The objects pulled out of the
"""
def __init__(self, filepath):
self.filepath = filepath # For record
self.classname = None
self.object_name = None
self.buffer = None
self.version = None
self.names = None
self.segments = None
self.heap = None
self.attrs = dict()
self._init_buffer()
self._init_heap()
self._read_header()
self.parse_file()
def __getitem__(self, item):
return self.attrs[item]
def keys(self):
"""Returns the names of the available properties"""
return self.attrs.keys()
def get(self, item, default=None):
return self.attrs.get(item, default)
def _init_buffer(self):
fp = open(self.filepath, "rb")
rdr = MatFile5Reader(fp, struct_as_record=True, squeeze_me=True)
vars = rdr.get_variables()
self.classname = vars["None"]["s2"][0].decode("utf-8")
self.object_name = vars["None"]["s0"][0].decode("utf-8")
fws = vars["__function_workspace__"]
self.buffer = BytesIO(fws.tostring())
fp.close()
def _init_heap(self):
super_data = read_workspace_vars(self.filepath)
elem = super_data["var_0"][0, 0]
if isinstance(elem, mat_struct):
self.heap = elem.MCOS[0]["arr"]
else:
self.heap = elem["MCOS"][0]["arr"]
def _read_header(self):
self.buffer.seek(248) # the start of the header
version = read_int(self.buffer)
n_str = read_int(self.buffer)
offsets = read_int(self.buffer, n=6)
# check that the next two are zeros
reserved = read_int(self.buffer, n=2)
assert all(
[x == 0 for x in reserved]
), "Non-zero reserved header fields: {}".format(reserved)
# check that we are at the right place
assert self.buffer.tell() == 288, "String elemnts begin at 288"
hdrs = []
for i in range(n_str):
hdrs.append(read_string(self.buffer))
self.names = hdrs
self.version = version
# The offsets are actually STARTING FROM 248 as well
self.segments = [x + 248 for x in offsets] # list(offsets)
return
def parse_file(self):
# Get class attributes from segment 1
self.buffer.seek(self.segments[0])
classes = self._parse_class_attributes(self.segments[1])
# Get first set of properties from segment 2
self.buffer.seek(self.segments[1])
props1 = self._parse_properties(self.segments[2])
# Get the property description from segment 3
self.buffer.seek(self.segments[2])
object_info = self._parse_prop_description(classes, self.segments[3])
# Get more properties from segment 4
self.buffer.seek(self.segments[3])
props2 = self._parse_properties(self.segments[4])
# Check that the last segment is empty
self.buffer.seek(self.segments[4])
seg5_length = (self.segments[5] - self.segments[4]) // 8
read_delim(self.buffer, seg5_length)
props = (props1, props2)
self._to_attrs(object_info, props)
def _to_attrs(self, object_info, props):
"""Re-organise the various classes and subclasses into a nested
dictionary.
:return:
"""
for pkg_clss, indices, idx in object_info:
pkg, clss = pkg_clss
idx = max(indices)
which = indices.index(idx)
obj = flatten_obj(props[which][idx])
subdict = self.attrs
if pkg != "":
subdict = self.attrs.setdefault(pkg, {})
if clss in subdict:
if isinstance(subdict[clss], list):
subdict[clss].append(obj)
else:
subdict[clss] = [subdict[clss]]
subdict[clss].append(obj)
else:
subdict[clss] = obj
def describe(self):
describe(self.attrs)
def _parse_class_attributes(self, section_end):
"""Read the Class attributes = the first segment"""
read_delim(self.buffer, 4)
classes = []
while self.buffer.tell() < section_end:
package_index = read_int(self.buffer) - 1
package = self.names[package_index] if package_index > 0 else ""
name_idx = read_int(self.buffer) - 1
name = self.names[name_idx] if name_idx > 0 else ""
classes.append((package, name))
read_delim(self.buffer, 2)
return classes
def _parse_prop_description(self, classes, section_end):
"""Parse the description of each property = the third segment"""
read_delim(self.buffer, 6)
object_info = []
while self.buffer.tell() < section_end:
class_idx = read_int(self.buffer) - 1
class_type = classes[class_idx]
read_delim(self.buffer, 2)
indices = [x - 1 for x in read_int(self.buffer, 2)]
obj_id = read_int(self.buffer)
object_info.append((class_type, indices, obj_id))
return object_info
def _parse_properties(self, section_end):
"""
Parse the actual values of the attributes == segments 2 and 4
"""
read_delim(self.buffer, 2)
props = []
while self.buffer.tell() < section_end:
n_props = read_int(self.buffer)
d = parse_prop(n_props, self.buffer, self.names, self.heap)
if not d: # Empty dictionary
break
props.append(d)
# Move to next 8-byte aligned offset
self.buffer.seek(self.buffer.tell() + self.buffer.tell() % 8)
return props
def to_hdf(self, filename):
f = h5py.File(filename, mode="w")
save_to_hdf(f, "/", self.attrs)
def describe(d, indent=0, width=4, out=None):
for key, value in d.items():
print(f'{"": <{width * indent}}' + str(key), file=out)
if isinstance(value, dict):
describe(value, indent + 1, out=out)
elif isinstance(value, np.ndarray):
print(
f'{"": <{width * (indent + 1)}} {value.shape} array '
f"of type {value.dtype}",
file=out,
)
elif isinstance(value, scipy.sparse.csc.csc_matrix):
print(
f'{"": <{width * (indent + 1)}} {value.shape} '
f"sparse matrix of type {value.dtype}",
file=out,
)
elif isinstance(value, Iterable) and not isinstance(value, str):
print(
f'{"": <{width * (indent + 1)}} {type(value)} of len ' f"{len(value)}",
file=out,
)
else:
print(f'{"": <{width * (indent + 1)}} {value}', file=out)
def parse_prop(n_props, buff, names, heap):
d = dict()
for i in range(n_props):
name_idx, flag, heap_idx = read_int(buff, 3)
if flag not in [0, 1, 2] and name_idx == 0:
n_props = flag
buff.seek(buff.tell() - 1) # go back on one byte
d = parse_prop(n_props, buff, names, heap)
else:
item_name = names[name_idx - 1]
if flag == 0:
d[item_name] = names[heap_idx]
elif flag == 1:
d[item_name] = heap[heap_idx + 2] # Todo: what is the heap?
elif flag == 2:
assert 0 <= heap_idx <= 1, (
"Boolean flag has a value other " "than 0 or 1 "
)
d[item_name] = bool(heap_idx)
else:
raise ValueError(
"unknown flag {} for property {} with heap "
"index {}".format(flag, item_name, heap_idx)
)
return d
def is_object(x):
"""Checking object dtype for structured numpy arrays"""
if x.dtype.names is not None and len(x.dtype.names) > 1: # Complex obj
return all(x.dtype[ix] == np.object for ix in range(len(x.dtype)))
else: # simple object
return x.dtype == np.object
def flatten_obj(arr):
# TODO turn structured arrays into nested dicts of lists rather that
# lists of dicts
if isinstance(arr, np.ndarray):
if arr.dtype.names:
arrdict = dict()
for fieldname in arr.dtype.names:
arrdict[fieldname] = flatten_obj(arr[fieldname])
arr = arrdict
elif arr.dtype == np.object and arr.ndim == 0:
arr = flatten_obj(arr[()])
elif arr.dtype == np.object and arr.ndim > 0:
try:
arr = np.stack(arr)
if arr.dtype.names:
d = {k: flatten_obj(arr[k]) for k in arr.dtype.names}
arr = d
except:
arr = [flatten_obj(x) for x in arr.tolist()]
elif isinstance(arr, dict):
arr = {k: flatten_obj(v) for k, v in arr.items()}
elif isinstance(arr, list):
try:
arr = flatten_obj(np.stack(arr))
except:
arr = [flatten_obj(x) for x in arr]
return arr
def save_to_hdf(h5file, path, dic):
"""
Saving a MATLAB object to HDF5
"""
if isinstance(dic, list):
dic = {str(i): v for i, v in enumerate(dic)}
for key, item in dic.items():
if isinstance(item, (int, float, str)):
h5file[path].attrs.create(key, item)
elif isinstance(item, list):
if len(item) == 0 and path + key not in h5file: # empty list empty group
h5file.create_group(path + key)
if all(isinstance(x, (int, float, str)) for x in item):
if path not in h5file:
h5file.create_group(path)
h5file[path].attrs.create(key, item)
else:
if path + key not in h5file:
h5file.create_group(path + key)
save_to_hdf(
h5file, path + key + "/", {str(i): x for i, x in enumerate(item)}
)
elif isinstance(item, scipy.sparse.csc.csc_matrix):
try:
h5file.create_dataset(
path + key, data=item.todense(), compression="gzip"
)
except Exception as e:
print(path + key)
raise e
elif isinstance(item, (np.ndarray, np.int64, np.float64)):
if item.dtype == np.dtype("<U1"): # Strings to 'S' type for HDF5
item = item.astype("S")
try:
h5file.create_dataset(path + key, data=item, compression="gzip")
except Exception as e:
print(path + key)
raise e
elif isinstance(item, dict):
if path + key not in h5file:
h5file.create_group(path + key)
save_to_hdf(h5file, path + key + "/", item)
elif item is None:
continue
else:
raise ValueError(f"Cannot save {type(item)} type at key {path + key}")
## NOT YET FULLY IMPLEMENTED!
class _Info:
def __init__(self, info):
self.info = info
self._identity = None
def __getitem__(self, item):
val = self.info[item]
if val.shape[0] == 1:
val = val[0]
if 0 in val[1].shape:
val = val[0]
if isinstance(val, scipy.sparse.csc.csc_matrix):
return np.asarray(val.todense())
if val.dtype == np.dtype("O"):
# 3d "sparse matrix"
if all(isinstance(x, scipy.sparse.csc.csc_matrix) for x in val):
val = np.array([x.todense() for x in val])
# TODO: The actual object data
equality = val[0] == val[1]
if isinstance(equality, scipy.sparse.csc.csc_matrix):
equality = equality.todense()
if equality.all():
val = val[0]
return np.squeeze(val)
@property
def categories(self):
return self.info.dtype.names
class TrapInfo(_Info):
def __init__(self, info):
"""
The information on all of the traps in a given position.
:param info: The TrapInfo structure, can be found in the heap of
the CTimelapse at index 7
"""
super().__init__(info)
class CellInfo(_Info):
def __init__(self, info):
"""
The extracted information of all cells in a given position.
:param info: The CellInfo structure, can be found in the heap
of the CTimelapse at index 15.
"""
super().__init__(info)
@property
def identity(self):
if self._identity is None:
self._identity = pd.DataFrame(
zip(self["trapNum"], self["cellNum"]), columns=["trapNum", "cellNum"]
)
return self._identity
def index(self, trapNum, cellNum):
query = "trapNum=={} and cellNum=={}".format(trapNum, cellNum)
try:
result = self.identity.query(query).index[0]
except Exception as e:
print(query)
raise e
return result
@property
def nucEstConv1(self):
return np.asarray(self.info["nuc_est_conv"][0][0].todense())
@property
def nucEstConv2(self):
return np.asarray(self.info["nuc_est_conv"][0][1].todense())
@property
def mothers(self):
return np.where((self["births"] != 0).any(axis=1))[0]
def daughters(self, mother_index):
"""
Get daughters of cell with index `mother_index`.
:param mother_index: the index of the mother within the data. This is
different from the mother's cell/trap identity.
"""
daughter_ids = np.unique(self["daughterLabel"][mother_index]).tolist()
daughter_ids.remove(0)
mother_trap = self.identity["trapNum"].loc[mother_index]
daughters = [self.index(mother_trap, cellNum) for cellNum in daughter_ids]
return daughters
def _todict(matobj):
"""
A recursive function which constructs from matobjects nested dictionaries
"""
if not hasattr(matobj, "_fieldnames"):
return matobj
d = {}
for strg in matobj._fieldnames:
elem = matobj.__dict__[strg]
if isinstance(elem, matlab.mio5_params.mat_struct):
d[strg] = _todict(elem)
elif isinstance(elem, np.ndarray):
d[strg] = _toarray(elem)
else:
d[strg] = elem
return d
def _toarray(ndarray):
"""
A recursive function which constructs ndarray from cellarrays
(which are loaded as numpy ndarrays), recursing into the elements
if they contain matobjects.
"""
if ndarray.dtype != "float64":
elem_list = []
for sub_elem in ndarray:
if isinstance(sub_elem, matlab.mio5_params.mat_struct):
elem_list.append(_todict(sub_elem))
elif isinstance(sub_elem, np.ndarray):
elem_list.append(_toarray(sub_elem))
else:
elem_list.append(sub_elem)
return np.array(elem_list)
else:
return ndarray
from pathlib import Path
class Strain:
"""The cell info for all the positions of a strain."""
def __init__(self, origin, strain):
self.origin = Path(origin)
self.files = [x for x in origin.iterdir() if strain in str(x)]
self.cts = [matObject(x) for x in self.files]
self.cinfos = [CellInfo(x.heap[15]) for x in self.cts]
self._identity = None
def __getitem__(self, item):
try:
return np.concatenate([c[item] for c in self.cinfos])
except ValueError: # If first axis is the channel
return np.concatenate([c[item] for c in self.cinfos], axis=1)
@property
def categories(self):
return set.union(*[set(c.categories) for c in self.cinfos])
@property
def identity(self):
if self._identity is None:
identities = []
for pos_id, cinfo in enumerate(self.cinfos):
identity = cinfo.identity
identity["position"] = pos_id
identities.append(identity)
self._identity = pd.concat(identities, ignore_index=True)
return self._identity
def index(self, posNum, trapNum, cellNum):
query = "position=={} and trapNum=={} and cellNum=={}".format(
posNum, trapNum, cellNum
)
try:
result = self.identity.query(query).index[0]
except Exception as e:
raise e
return result
@property
def mothers(self):
# At least two births are needed to be considered a mother cell
return np.where(np.count_nonzero(self["births"], axis=1) > 3)[0]
def daughters(self, mother_index):
"""
Get daughters of cell with index `mother_index`.
:param mother_index: the index of the mother within the data. This is
different from the mother's pos/trap/cell identity.
"""
daughter_ids = np.unique(self["daughterLabel"][mother_index]).tolist()
if 0 in daughter_ids:
daughter_ids.remove(0)
mother_pos_trap = self.identity[["position", "trapNum"]].loc[mother_index]
daughters = []
for cellNum in daughter_ids:
try:
daughters.append(self.index(*mother_pos_trap, cellNum))
except IndexError:
continue
return daughters
aliby/io/omero.py deleted 100644 → 0
View file @ e8e2a3bc
import numpy as np
import dask.array as da
from dask import delayed
from omero.gateway import BlitzGateway
from omero.model import enums as omero_enums
# 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 Argo:
"""
Base class to interact with OMERO.
See
https://docs.openmicroscopy.org/omero/5.6.0/developers/Python.html
"""
def __init__(
self,
host="islay.bio.ed.ac.uk",
username="upload",
password="***REMOVED***",
):
"""
Parameters
----------
host : string
web address of OMERO host
username: string
password : string
"""
self.conn = None
self.host = host
self.username = username
self.password = password
# standard method required for Python's with statement
def __enter__(self):
self.conn = BlitzGateway(
host=self.host, username=self.username, passwd=self.password
)
self.conn.connect()
self.conn.c.enableKeepAlive(60)
return self
# standard method required for Python's with statement
def __exit__(self, *exc):
self.conn.close()
return 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)
aliby/io/utils.py deleted 100644 → 0
View file @ e8e2a3bc
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"
aliby/multiexperiment.py deleted 100644 → 0
View file @ e8e2a3bc
"""Experimental methods to deal with multiple experiments and data aggregation."""
from pathos.multiprocessing import Pool
from aliby.pipeline import PipelineParameters, Pipeline
class MultiExp:
"""
Manages cases when you need to segment several different experiments with a single
position (e.g. pH calibration).
"""
def __init__(self, expt_ids, npools=8, *args, **kwargs):
self.expt_ids = expt_ids
def run(self):
run_expt = lambda expt: Pipeline(
PipelineParameters.default(general={"expt_id": expt, "distributed": 0})
).run()
with Pool(npools) as p:
results = p.map(lambda x: self.create_pipeline(x), self.exp_ids)
@classmethod
def default(self):
return cls(expt_ids=list(range(20448, 20467 + 1)))
aliby/post_processing.py deleted 100644 → 0
View file @ e8e2a3bc
"""
Post-processing utilities
Notes: I don't have statistics on ranges of radii for each of the knots in
the radial spline representation, but we regularly extract the average of
these radii for each cell. So, depending on camera/lens, we get:
* 60x evolve: mean radii of 2-14 pixels (and measured areas of 30-750
pixels^2)
* 60x prime95b: mean radii of 3-24 pixels (and measured areas of 60-2000
pixels^2)
And I presume that for a 100x lens we would get an ~5/3 increase over those
values.
In terms of the current volume estimation method, it's currently only
implemented in the AnalysisToolbox repository, but it's super simple:
mVol = 4/3*pi*sqrt(mArea/pi).^3
where mArea is simply the sum of pixels for that cell.
"""
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy import ndimage
from skimage.morphology import erosion, ball
from skimage import measure, draw
def my_ball(radius):
"""Generates a ball-shaped structuring element.
This is the 3D equivalent of a disk.
A pixel is within the neighborhood if the Euclidean distance between
it and the origin is no greater than radius.
Parameters
----------
radius : int
The radius of the ball-shaped structuring element.
Other Parameters
----------------
dtype : data-type
The data type of the structuring element.
Returns
-------
selem : ndarray
The structuring element where elements of the neighborhood
are 1 and 0 otherwise.
"""
n = 2 * radius + 1
Z, Y, X = np.mgrid[-radius:radius:n * 1j,
-radius:radius:n * 1j,
-radius:radius:n * 1j]
X **= 2
Y **= 2
Z **= 2
X += Y
X += Z
# s = X ** 2 + Y ** 2 + Z ** 2
return X <= radius * radius
def circle_outline(r):
return ellipse_perimeter(r, r)
def ellipse_perimeter(x, y):
im_shape = int(2*max(x, y) + 1)
img = np.zeros((im_shape, im_shape), dtype=np.uint8)
rr, cc = draw.ellipse_perimeter(int(im_shape//2), int(im_shape//2),
int(x), int(y))
img[rr, cc] = 1
return np.pad(img, 1)
def capped_cylinder(x, y):
max_size = (y + 2*x + 2)
pixels = np.zeros((max_size, max_size))
rect_start = ((max_size-x)//2, x + 1)
rr, cc = draw.rectangle_perimeter(rect_start, extent=(x, y),
shape=(max_size, max_size))
pixels[rr, cc] = 1
circle_centres = [(max_size//2 - 1, x),
(max_size//2 - 1, max_size - x - 1 )]
for r, c in circle_centres:
rr, cc = draw.circle_perimeter(r, c, (x + 1)//2,
shape=(max_size, max_size))
pixels[rr, cc] = 1
pixels = ndimage.morphology.binary_fill_holes(pixels)
pixels ^= erosion(pixels)
return pixels
def volume_of_sphere(radius):
return 4 / 3 * np.pi * radius**3
def plot_voxels(voxels):
verts, faces, normals, values = measure.marching_cubes_lewiner(
voxels, 0)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts[faces])
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)
ax.set_xlim(0, voxels.shape[0])
ax.set_ylim(0, voxels.shape[1])
ax.set_zlim(0, voxels.shape[2])
plt.tight_layout()
plt.show()
# Volume estimation
def union_of_spheres(outline, shape='my_ball', debug=False):
filled = ndimage.binary_fill_holes(outline)
nearest_neighbor = ndimage.morphology.distance_transform_edt(
outline == 0) * filled
voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
c_z = voxels.shape[2] // 2
for x,y in zip(*np.where(filled)):
radius = nearest_neighbor[(x,y)]
if radius > 0:
if shape == 'ball':
b = ball(radius)
elif shape == 'my_ball':
b = my_ball(radius)
else:
raise ValueError(f"{shape} is not an accepted value for "
f"shape.")
centre_b = ndimage.measurements.center_of_mass(b)
I,J,K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
K + int(c_z - centre_b[2])] += b
if debug:
plot_voxels(voxels)
return voxels.astype(bool).sum()
def improved_uos(outline, shape='my_ball', debug=False):
filled = ndimage.binary_fill_holes(outline)
nearest_neighbor = ndimage.morphology.distance_transform_edt(
outline == 0) * filled
voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
c_z = voxels.shape[2] // 2
while np.any(nearest_neighbor != 0):
radius = np.max(nearest_neighbor)
x, y = np.argwhere(nearest_neighbor == radius)[0]
if shape == 'ball':
b = ball(np.ceil(radius))
elif shape == 'my_ball':
b = my_ball(np.ceil(radius))
else:
raise ValueError(f"{shape} is not an accepted value for shape")
centre_b = ndimage.measurements.center_of_mass(b)
I, J, K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
K + int(c_z - centre_b[2])] += b
# Use the central disk of the ball from voxels to get the circle
# = 0 if nn[x,y] < r else nn[x,y]
rr, cc = draw.circle(x, y, np.ceil(radius), nearest_neighbor.shape)
nearest_neighbor[rr, cc] = 0
if debug:
plot_voxels(voxels)
return voxels.astype(bool).sum()
def conical(outline, debug=False):
nearest_neighbor = ndimage.morphology.distance_transform_edt(
outline == 0) * ndimage.binary_fill_holes(outline)
if debug:
hf = plt.figure()
ha = hf.add_subplot(111, projection='3d')
X, Y = np.meshgrid(np.arange(nearest_neighbor.shape[0]),
np.arange(nearest_neighbor.shape[1]))
ha.plot_surface(X, Y, nearest_neighbor)
plt.show()
return 4 * nearest_neighbor.sum()
def volume(outline, method='spheres'):
if method=='conical':
return conical(outline)
elif method=='spheres':
return union_of_spheres(outline)
else:
raise ValueError(f"Method {method} not implemented.")
def circularity(outline):
pass
\ No newline at end of file
aliby/results.py deleted 100644 → 0
View file @ e8e2a3bc
"""Pipeline results classes and utilities"""
class SegmentationResults:
"""
Object storing the data from the Segmentation pipeline.
Everything is stored as an `AttributeDict`, which is a `defaultdict` where
you can get elements as attributes.
In addition, it implements:
- IO functionality (read from file, write to file)
"""
def __init__(self, raw_expt):
pass
class CellResults:
"""
Results on a set of cells TODO: what set of cells, how?
Contains:
* cellInf describing which cells are taken into account
* annotations on the cell
* segmentation maps of the cell TODO: how to define and save this?
* trapLocations TODO: why is this not part of cellInf?
"""
def __init__(self, cellInf=None, annotations=None, segmentation=None,
trapLocations=None):
self._cellInf = cellInf
self._annotations = annotations
self._segmentation = segmentation
self._trapLocations = trapLocations
aliby/timelapse.py deleted 100644 → 0
View file @ e8e2a3bc
"""
Classes to handle multidimensional images, both remotely and local.
"""
import itertools
import logging
import h5py
import numpy as np
from pathlib import Path
from tqdm import tqdm
import cv2
from aliby.io.matlab import matObject
from agora.io.utils import Cache, imread, get_store_path
logger = logging.getLogger(__name__)
def parse_local_fs(pos_dir, tp=None):
"""
Local file structure:
- pos_dir
-- exptID_{timepointID}_{ChannelID}_{z_position_id}.png
:param pos_dirs:
:return: Image_mapper
"""
pos_dir = Path(pos_dir)
img_mapper = dict()
def channel_idx(img_name):
return img_name.stem.split("_")[-2]
def tp_idx(img_name):
return int(img_name.stem.split("_")[-3]) - 1
def z_idx(img_name):
return img_name.stem.split("_")[-1]
if tp is not None:
img_list = [img for img in pos_dir.iterdir() if tp_idx(img) in tp]
else:
img_list = [img for img in pos_dir.iterdir()]
for tp, group in itertools.groupby(sorted(img_list, key=tp_idx), key=tp_idx):
img_mapper[int(tp)] = {
channel: {i: item for i, item in enumerate(sorted(grp, key=z_idx))}
for channel, grp in itertools.groupby(
sorted(group, key=channel_idx), key=channel_idx
)
}
return img_mapper
class Timelapse:
"""
Timelapse class contains the specifics of one position.
"""
def __init__(self):
self._id = None
self._name = None
self._channels = []
self._size_c = 0
self._size_t = 0
self._size_x = 0
self._size_y = 0
self._size_z = 0
self.image_cache = None
self.annotation = None
def __repr__(self):
return self.name
def full_mask(self):
return np.full(self.shape, False)
def __getitem__(self, item):
cached = self.image_cache[item]
# Check if there are missing values, if so reload
# TODO only reload missing
mask = np.isnan(cached)
if np.any(mask):
full = self.load_fn(item)
shape = self.image_cache[
item
].shape # TODO speed this up by recognising the shape from the item
self.image_cache[item] = np.reshape(full, shape)
return full
return cached
def get_hypercube(self):
pass
def load_fn(self, item):
"""
The hypercube is ordered as: C, T, X, Y, Z
:param item:
:return:
"""
def parse_slice(s):
step = s.step if s.step is not None else 1
if s.start is None and s.stop is None:
return None
elif s.start is None and s.stop is not None:
return range(0, s.stop, step)
elif s.start is not None and s.stop is None:
return [s.start]
else: # both s.start and s.stop are not None
return range(s.start, s.stop, step)
def parse_subitem(subitem, kw):
if isinstance(subitem, (int, float)):
res = [int(subitem)]
elif isinstance(subitem, list) or isinstance(subitem, tuple):
res = list(subitem)
elif isinstance(subitem, slice):
res = parse_slice(subitem)
else:
res = subitem
# raise ValueError(f"Cannot parse slice {kw}: {subitem}")
if kw in ["x", "y"]:
# Need exactly two values
if res is not None:
if len(res) < 2:
# An int was passed, assume it was
res = [res[0], self.size_x]
elif len(res) > 2:
res = [res[0], res[-1] + 1]
return res
if isinstance(item, int):
return self.get_hypercube(
x=None, y=None, z_positions=None, channels=[item], timepoints=None
)
elif isinstance(item, slice):
return self.get_hypercube(channels=parse_slice(item))
keywords = ["channels", "timepoints", "x", "y", "z_positions"]
kwargs = dict()
for kw, subitem in zip(keywords, item):
kwargs[kw] = parse_subitem(subitem, kw)
return self.get_hypercube(**kwargs)
@property
def shape(self):
return (self.size_c, self.size_t, self.size_x, self.size_y, self.size_z)
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def size_z(self):
return self._size_z
@property
def size_c(self):
return self._size_c
@property
def size_t(self):
return self._size_t
@property
def size_x(self):
return self._size_x
@property
def size_y(self):
return self._size_y
@property
def channels(self):
return self._channels
def get_channel_index(self, channel):
return self.channels.index(channel)
def load_annotation(filepath: Path):
try:
return matObject(filepath)
except Exception as e:
raise (
"Could not load annotation file. \n"
"Non MATLAB files currently unsupported"
) from e
class TimelapseOMERO(Timelapse):
"""
Connected to an Image object which handles database I/O.
"""
def __init__(self, image, annotation, cache, **kwargs):
super(TimelapseOMERO, self).__init__()
self.image = image
# Pre-load pixels
self.pixels = self.image.getPrimaryPixels()
self._id = self.image.getId()
self._name = self.image.getName()
self._size_x = self.image.getSizeX()
self._size_y = self.image.getSizeY()
self._size_z = self.image.getSizeZ()
self._size_c = self.image.getSizeC()
self._size_t = self.image.getSizeT()
self._channels = self.image.getChannelLabels()
# Check whether there are file annotations for this position
if annotation is not None:
self.annotation = load_annotation(annotation)
# Get an HDF5 dataset to use as a cache.
compression = kwargs.get("compression", None)
self.image_cache = cache.require_dataset(
self.name,
self.shape,
dtype=np.float16,
fillvalue=np.nan,
compression=compression,
)
def get_hypercube(
self, x=None, y=None, z_positions=None, channels=None, timepoints=None
):
if x is None and y is None:
tile = None # Get full plane
elif x is None:
ymin, ymax = y
tile = (None, ymin, None, ymax - ymin)
elif y is None:
xmin, xmax = x
tile = (xmin, None, xmax - xmin, None)
else:
xmin, xmax = x
ymin, ymax = y
tile = (xmin, ymin, xmax - xmin, ymax - ymin)
if z_positions is None:
z_positions = range(self.size_z)
if channels is None:
channels = range(self.size_c)
if timepoints is None:
timepoints = range(self.size_t)
z_positions = z_positions or [0]
channels = channels or [0]
timepoints = timepoints or [0]
zcttile_list = [
(z, c, t, tile)
for z, c, t in itertools.product(z_positions, channels, timepoints)
]
planes = list(self.pixels.getTiles(zcttile_list))
order = (
len(z_positions),
len(channels),
len(timepoints),
planes[0].shape[-2],
planes[0].shape[-1],
)
result = np.stack([x for x in planes]).reshape(order)
# Set to C, T, X, Y, Z order
result = np.moveaxis(result, -1, -2)
return np.moveaxis(result, 0, -1)
def cache_set(self, save_dir, timepoints, expt_name, quiet=True):
# TODO deprecate when this is default
pos_dir = save_dir / self.name
if not pos_dir.exists():
pos_dir.mkdir()
for tp in tqdm(timepoints, desc=self.name):
for channel in tqdm(self.channels, disable=quiet):
for z_pos in tqdm(range(self.size_z), disable=quiet):
ch_id = self.get_channel_index(channel)
image = self.get_hypercube(
x=None,
y=None,
channels=[ch_id],
z_positions=[z_pos],
timepoints=[tp],
)
im_name = "{}_{:06d}_{}_{:03d}.png".format(
expt_name, tp + 1, channel, z_pos + 1
)
cv2.imwrite(str(pos_dir / im_name), np.squeeze(image))
# TODO update positions table to get the number of timepoints?
return list(itertools.product([self.name], timepoints))
def run(self, keys, store, save_dir="./", **kwargs):
"""
Parse file structure and get images for the timepoints in keys.
"""
save_dir = Path(save_dir)
if keys is None:
# TODO save final metadata
return None
store = save_dir / store
# A position specific store
store = store.with_name(self.name + store.name)
# Create store if it does not exist
if not store.exists():
# The first run, add metadata to the store
with h5py.File(store, "w") as pos_store:
# TODO Add metadata to the store.
pass
# TODO check how sensible the keys are with what is available
# if some of the keys don't make sense, log a warning and remove
# them so that the next steps of the pipeline make sense
return keys
def clear_cache(self):
self.image_cache.clear()
class TimelapseLocal(Timelapse):
def __init__(
self, position, root_dir, finished=True, annotation=None, cache=None, **kwargs
):
"""
Linked to a local directory containing the images for one position
in an experiment.
Can be a still running experiment or a finished one.
:param position: Name of the position
:param root_dir: Root directory
:param finished: Whether the experiment has finished running or the
class will be used as part of a pipeline, mostly with calls to `run`
"""
super(TimelapseLocal, self).__init__()
self.pos_dir = Path(root_dir) / position
assert self.pos_dir.exists()
self._id = position
self._name = position
if finished:
self.image_mapper = parse_local_fs(self.pos_dir)
self._update_metadata()
else:
self.image_mapper = dict()
self.annotation = None
# Check whether there are file annotations for this position
if annotation is not None:
self.annotation = load_annotation(annotation)
compression = kwargs.get("compression", None)
self.image_cache = cache.require_dataset(
self.name,
self.shape,
dtype=np.float16,
fillvalue=np.nan,
compression=compression,
)
def _update_metadata(self):
self._size_t = len(self.image_mapper)
# Todo: if cy5 is the first one it causes issues with getting x, y
# hence the sorted but it's not very robust
self._channels = sorted(
list(set.union(*[set(tp.keys()) for tp in self.image_mapper.values()]))
)
self._size_c = len(self._channels)
# Todo: refactor so we don't rely on there being any images at all
self._size_z = max([len(self.image_mapper[0][ch]) for ch in self._channels])
single_img = self.get_hypercube(
x=None, y=None, z_positions=None, channels=[0], timepoints=[0]
)
self._size_x = single_img.shape[2]
self._size_y = single_img.shape[3]
def get_hypercube(
self, x=None, y=None, z_positions=None, channels=None, timepoints=None
):
xmin, xmax = x if x is not None else (None, None)
ymin, ymax = y if y is not None else (None, None)
if z_positions is None:
z_positions = range(self.size_z)
if channels is None:
channels = range(self.size_c)
if timepoints is None:
timepoints = range(self.size_t)
def z_pos_getter(z_positions, ch_id, t):
default = np.zeros((self.size_x, self.size_y))
names = [
self.image_mapper[t][self.channels[ch_id]].get(i, None)
for i in z_positions
]
res = [imread(name) if name is not None else default for name in names]
return res
# nested list of images in C, T, X, Y, Z order
ctxyz = []
for ch_id in channels:
txyz = []
for t in timepoints:
xyz = z_pos_getter(z_positions, ch_id, t)
txyz.append(np.dstack(list(xyz))[xmin:xmax, ymin:ymax])
ctxyz.append(np.stack(txyz))
return np.stack(ctxyz)
def clear_cache(self):
self.image_cache.clear()
def run(self, keys, store, save_dir="./", **kwargs):
"""
Parse file structure and get images for the time points in keys.
"""
if keys is None:
return None
elif isinstance(keys, int):
keys = [keys]
self.image_mapper.update(parse_local_fs(self.pos_dir, tp=keys))
self._update_metadata()
# Create store if it does not exist
store = get_store_path(save_dir, store, self.name)
if not store.exists():
# The first run, add metadata to the store
with h5py.File(store, "w") as pos_store:
# TODO Add metadata to the store.
pass
# TODO check how sensible the keys are with what is available
# if some of the keys don't make sense, log a warning and remove
# them so that the next steps of the pipeline make sense
return keys
docs/requirements.txt
View file @ 95958192
numpydoc>=1.3.1
aliby>=0.1.25
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
docs/source/INSTALL.md
View file @ 95958192
# Installation
Tested on: Mac OSX Mojave and Ubuntu 20.04
# Installation
## Requirements
We strongly recommend installing within a python environment as there are many dependencies that you may not want polluting your regular python environment.
We strongly recommend installing within a python environment as there are many dependencies that you may not want polluting your regular python environment.
Make sure you are using python 3.
An environment can be created with using the conda package manager:
An environment can be created using [Anaconda](https://www.anaconda.com/):
$ conda create --name <env>
$ conda activate <env>
Which you can deactivate with:
Which you can deactivate with:
$ conda deactivate
Or using virtualenv:
Or using virtualenv:
$ python -m virtualenv /path/to/venv/
$ source /path/to/venv/bin/activate
This will download all of your packages under `/path/to/venv` and then activate it.
Deactivate using
This will download all of your packages under `/path/to/venv` and then activate it.
Deactivate using
$ deactivate
You will also need to make sure you have a recent version of pip.
In your local environment, run:
In your local environment, run:
$ pip install --upgrade pip
Or using [pyenv](https://github.com/pyenv/pyenv) with pyenv-virtualenv:
$ pyenv install 3.7.9
$ pyenv virtualenv 3.7.9 aliby
$ pyenv install 3.8.14
$ pyenv virtualenv 3.8.14 aliby
$ pyenv local aliby
## Pipeline installation
### Pip version
Once you have created your local environment, run:
Once you have created and activated your virtual environment, run:
If you are not using an OMERO server setup:
$ pip install aliby
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).
$ cd aliby
$ pip install -e ./
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
For local access and processing, follow the same instructions as Linux. Remote access to OMERO servers depends on some issues in one of our depedencies being solved (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@gitlab.com/aliby/aliby.git
$ cd aliby
and then either
$$ poetry install --all-extras
for everything, including tools to access OMERO servers, or
$$ poetry install
for a version with only local access, or
$$ poetry install --with dev
to install with compatible versions of the development tools we use, such as black.
These commands 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
Segmentation has been tested on: Mac OSX Mojave, Ubuntu 20.04 and Arch Linux.
Data processing has been tested on all the above and Windows 11.
### Detailed Windows installation
#### Create environment
Open anaconda powershell as administrator
```shell script
conda create -n devaliby2 -c conda-forge python=3.8 omero-py
conda activate devaliby2
```
#### Install poetry
You may have to specify the python executable to get this to work :
```shell script
(Invoke-WebRequest -Uri https://install.python-poetry.org -UseBasicParsing).Content | C:\Users\USERNAME\Anaconda3\envs\devaliby2\python.exe -
``` Also specify full path when running poetry (there must be a way to sort this)
- Clone the repository (Assuming you have ssh properly set up)
```shell script
git clone git@gitlab.com:aliby/aliby.git
cd aliby
poetry install --all-extras
```
You may need to run the full poetry path twice - first time gave an error message, worked second time
```shell script
C:\Users\v1iclar2\AppData\Roaming\Python\Scripts\poetry install --all-extras
```
In case you want to have local versions (usually for development) the main three aliby dependencies you must install them in a specific order:
confirm installation of aliby - python...import aliby - get no error message
$ git clone git@git.ecdf.ed.ac.uk:swain-lab/aliby/aliby.git
$ git clone git@git.ecdf.ed.ac.uk:swain-lab/aliby/postprocessor.git
$ git clone git@git.ecdf.ed.ac.uk:swain-lab/aliby/agora.git
#### Access the virtual environment from the IDE (e.g., PyCharm)
New project
In location - navigate to the aliby folder (eg c::/Users/Public/Repos/aliby
$ cd aliby && poetry install
$ cd ../postprocessor && poetry install
$ cd ../agora && poetry install
- Select the correct python interpreter
click the interpreter name at the bottom right
click add local interpreter
on the left click conda environment
click the 3 dots to the right of the interpreter path and navigate to the python executable from the environment created above (eg C:\Users\v1iclar2\Anaconda3\envs\devaliby2\python.exe)
And that should install all three main dependencies in an editable mode. The same process can be used for [BABY](https://git.ecdf.ed.ac.uk/swain-lab/aliby/baby)
#### Potential Windows issues
- Sometimes the library pywin32 gives trouble, just install it using pip or conda
docs/source/PIPELINE.md 0 → 100644
View file @ 95958192
# Running the analysis pipeline
You can run the analysis pipeline either via the command line interface (CLI) or using a script that incorporates the `aliby.pipeline.Pipeline` object.
## CLI
On a CLI, you can use the `aliby-run` command. This command takes options as follows:
- `--host`: Address of image-hosting server.
- `--username`: Username to access image-hosting server.
- `--password`: Password to access image-hosting server.
- `--expt_id`: Number ID of experiment stored on host server.
- `--distributed`: Number of distributed cores to use for segmentation and signal processing. If 0, there is no parallelisation.
- `--tps`: Optional. Number of time points from the beginning of the experiment to use. If not specified, the pipeline processes all time points.
- `--directory`: Optional. Parent directory to save the data files (HDF5) generated, `./data` by default; the files will be stored in a child directory whose name is the name of the experiment.
- `--filter`: Optional. List of positions to use for analysis. Alternatively, a regex (regular expression) or list of regexes to search for positions. **Note: for the CLI, currently it is not able to take a list of strings as input.**
- `--overwrite`: Optional. Whether to overwrite an existing data directory. True by default.
- `--override_meta`: Optional. Whether to overwrite an existing data directory. True by default.
Example usage:
```bash
aliby-run --expt_id EXPT_PATH --distributed 4 --tps None
```
And to run Omero servers, the basic arguments are shown:
```bash
aliby-run --expt_id XXX --host SERVER.ADDRESS --user USER --password PASSWORD
```
## Script
Use the `aliby.pipeline.Pipeline` object and supply a dictionary, following the example below. The meaning of the parameters are the same as described in the CLI section above.
```python
#!/usr/bin/env python3
from aliby.pipeline import Pipeline, PipelineParameters
# Specify experiment IDs
ids = [101, 102]
for i in ids:
print(i)
try:
params = PipelineParameters.default(
# Create dictionary to define pipeline parameters.
general={
"expt_id": i,
"distributed": 6,
"host": "INSERT ADDRESS HERE",
"username": "INSERT USERNAME HERE",
"password": "INSERT PASSWORD HERE",
# Ensure data will be overwriten
"override_meta": True,
"overwrite": True,
}
)
# Fine-grained control beyond general parameters:
# change specific leaf in the extraction tree.
# This example tells the pipeline to additionally compute the
# nuc_est_conv quantity, which is a measure of the degree of
# localisation of a signal in a cell.
params = params.to_dict()
leaf_to_change = params["extraction"]["tree"]["GFP"]["np_max"]
leaf_to_change.add("nuc_est_conv")
# Regenerate PipelineParameters
p = Pipeline(PipelineParameters.from_dict(params))
# Run pipeline
p.run()
# Error handling
except Exception as e:
print(e)
```
This example code can be the contents of a `run.py` file, and you can run it via
```bash
python run.py
```
in the appropriate virtual environment.
Alternatively, the example code can be the contents of a cell in a jupyter notebook.