rename the package as anscombe_numcodecs

.github/workflows/tests.yaml

@@ -16,17 +16,18 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: ["ubuntu-latest", "macos-latest"]
+        os: ["ubuntu-latest"]
+
     steps:
       - uses: actions/checkout@v2
       - uses: s-weigand/setup-conda@v1
-        with:
-          python-version: 3.9
+
       - name: Install dependencies
         run: |
           pip install -r requirements.txt
           pip install pytest
           pip install .
+
       - name: Test with pytest 
         run: |
           pytest -v

README.md

@@ -1,10 +1,10 @@
-[![PyPI version](https://badge.fury.io/py/poisson-numcodecs.svg)](https://badge.fury.io/py/poisson-numcodecs) ![tests](https://github.com/datajoint/poisson-numcodecs/actions/workflows/tests.yaml/badge.svg)
+[![PyPI version](https://badge.fury.io/py/anscombe-numcodecs.svg)](https://badge.fury.io/py/anscombe-numcodecs) ![tests](https://github.com/datajoint/anscombe-numcodecs/actions/workflows/tests.yaml/badge.svg)
 
-# Poisson - numcodecs implementation
+# Anscombe numcodecs
 
 This codec is designed for compressing movies with Poisson noise, which are produced by photon-limited modalities such multiphoton microscopy, radiography, and astronomy.
 
-The codec assumes that the video is linearly encoded with a potential offset (`zero_level`) and that the `photon_sensitivity` (the average increase in intensity per photon) is known or can be accurately estimated from the data.
+The codec assumes that the video is linearly encoded with a potential offset (`zero_level`) and that the `photon_sensitivity` (the average increase in intensity per photon) is either already known or can be accurately estimated from the data.
 
 The codec re-quantizes the grayscale efficiently with a square-root-like transformation to equalize the noise variance across the grayscale levels: the [Anscombe Transform](https://en.wikipedia.org/wiki/Anscombe_transform).
 This results in a smaller number of unique grayscale levels and significant improvements in the compressibility of the data without sacrificing signal accuracy.
@@ -20,16 +20,16 @@ The codec is used in Zarr as a filter prior to compression.
 Install via `pip`:
 
 ```
-pip install poisson-numcodecs
+pip install anscombe-numcodecs
 ```
 
 ### Developer installation
 
 ```
-conda create -n poisson_numcodecs python=3.xx
-conda activate poisson_numcodecs
-git clone https://github.com/AllenNeuralDynamics/poisson-numcodecs.git
-cd poisson-numcodecs
+conda create -n anscombe_numcodecs python=3.xx
+conda activate anscombe_numcodecs
+git clone https://github.com/datajoint/anscombe-numcodecs.git
+cd anscombe-numcodecs
 pip install -r requirements.txt
 pip install -e .
 ```

pyproject.toml

@@ -3,20 +3,37 @@ requires = ["setuptools>=61.0"]
 build-backend = "setuptools.build_meta"
 
 [project]
-name = "poisson_numcodecs"
-version = "0.1.0"
+name = "anscombe_numcodecs"
+version = "0.1.2"
 authors = [
-  { name="Jerome Lecoq", email="[email protected]" },
+  { name = "Jerome Lecoq", email = "[email protected]" },
+  { name = "Dimitri Yatsenko", email = "[email protected]" },
 ]
-description = "Numcodecs implementation of a poisson noise calibration."
+maintainers = [
+  { name = "Dimitri Yatsenko", email = "[email protected]" },
+]
+license = { file = "LICENSE" }
+description = "Numcodecs implementation for compressing photon-limited movies."
 readme = "README.md"
 requires-python = ">=3.8"
+keywords = ["video compression", "numcodecs", "zarr", "compression algorithms", "multimedia", "multiphoton microscopy"]
 classifiers = [
     "Programming Language :: Python :: 3",
     "License :: OSI Approved :: MIT License",
     "Operating System :: OS Independent",
+    "Topic :: Multimedia :: Video",
+    "Topic :: Software Development :: Libraries",
+]
+
+dependencies = [
+    "numpy>=2.0",
+    "numcodecs>=0.1",
+    "scikit-learn",
 ]
 
+# Dependencies will be managed via requirements.txt
+
 [project.urls]
-Homepage = "https://github.com/AllenNeuralDynamics/poisson-numcodecs"
-Issues = "https://github.com/AllenNeuralDynamics/poisson-numcodecs/issues"
+Homepage = "https://github.com/datajoint/anscombe-numcodecs"
+Source = "https://github.com/datajoint/anscombe-numcodecs"
+Issues = "https://github.com/datajoint/anscombe-numcodecs/issues"

requirements.txt

@@ -1,7 +1,9 @@
-numpy
-numcodecs
-zarr
-matplotlib
+zarr>=2.18
 scipy
-scikit-learn
 imageio
+matplotlib
+pytest>=7.0
+pytest-cov
+black
+flake8
+mypy

src/anscombe_numcodecs/__init__.py

@@ -0,0 +1,4 @@
+from .codec import AnscombeCodec
+from . import estimate 
+
+__all__ = ["AnscombeCodec", "estimate"]

src/poisson_numcodecs/codec.py → src/anscombe_numcodecs/codec.py

@@ -1,98 +1,109 @@
 """
 Numcodecs Codec implementation for Poisson noise calibration
 """
+
 import numpy as np
 import numcodecs
 from numcodecs.abc import Codec
 
+
 def make_anscombe_lookup(
-          sensitivity: float, 
-          input_max: int=0x7fff,
-          zero_level: int=0,
-          beta: float=0.5,
-          output_type='uint8'):
+    sensitivity: float,
+    input_max: int = 0x7FFF,
+    zero_level: int = 0,
+    beta: float = 0.5,
+    output_type="uint8",
+):
     """
     Compute the Anscombe lookup table.
-    The lookup converts a linear grayscale image into a uniform variance image. 
+    The lookup converts a linear grayscale image into a uniform variance image.
     :param sensitivity: the size of one photon in the linear input image.
     :param input_max: the maximum value in the input
     :param beta: the grayscale quantization step expressed in units of noise std dev
     """
-    xx = (np.r_[:input_max + 1] - zero_level) / sensitivity  # input expressed in photon rates
-    zero_slope = 1 / beta / np.sqrt(3/8)  # slope for negative values
+    xx = (
+        np.r_[: input_max + 1] - zero_level
+    ) / sensitivity  # input expressed in photon rates
+    zero_slope = 1 / beta / np.sqrt(3 / 8)  # slope for negative values
     offset = zero_level * zero_slope / sensitivity
-    lookup_table = np.round(offset +
-          (xx < 0) * (xx * zero_slope) +
-          (xx >= 0) * (2.0 / beta * (np.sqrt(np.maximum(0, xx) + 3/8) - np.sqrt(3/8))))
+    lookup_table = np.round(
+        offset
+        + (xx < 0) * (xx * zero_slope)
+        + (xx >= 0)
+        * (2.0 / beta * (np.sqrt(np.maximum(0, xx) + 3 / 8) - np.sqrt(3 / 8)))
+    )
     lookup = lookup_table.astype(output_type)
     assert np.diff(lookup_table).min() >= 0, "non-monotonic lookup generated"
     return lookup
 
 
-def make_inverse_lookup(lookup_table: np.ndarray, output_type='int16') -> np.ndarray:
+def make_inverse_lookup(lookup_table: np.ndarray, output_type="int16") -> np.ndarray:
     """Compute the inverse lookup table for a monotonic forward lookup table."""
     _, inv1 = np.unique(lookup_table, return_index=True)  # first entry
     _, inv2 = np.unique(lookup_table[::-1], return_index=True)  # last entry
-    inverse = (inv1 + lookup_table.size - 1 - inv2)/2
+    inverse = (inv1 + lookup_table.size - 1 - inv2) / 2
     return inverse.astype(output_type)
 
 
 def lookup(movie: np.ndarray, lookup_table: np.ndarray) -> np.ndarray:
     """Apply lookup table to movie"""
-    return lookup_table[np.maximum(0, np.minimum(movie, lookup_table.size-1))]
+    return lookup_table[np.maximum(0, np.minimum(movie, lookup_table.size - 1))]
 
 
-class PoissonCodec(Codec):
+class AnscombeCodec(Codec):
     """Codec for 3-dimensional Filter. The codec assumes that input data are of shape:
     (time, x, y).
 
     Parameters
     ----------
     zero_level : float
-        Signal level when no photons are recorded. 
+        Signal level when no photons are recorded.
         This should pre-computed or measured directly on the instrument.
     photon_sensitivity : float
         Conversion scalor to convert the measure signal into absolute photon numbers.
         This should pre-computed or measured directly on the instrument.
     """
+
     codec_id = "poisson"
 
-    def __init__(self, 
-                 zero_level, 
-                 photon_sensitivity, 
-                 encoded_dtype='int8', 
-                 decoded_dtype='int16',
-                 ): 
+    def __init__(
+        self,
+        zero_level,
+        photon_sensitivity,
+        encoded_dtype="int8",
+        decoded_dtype="int16",
+    ):
         self.zero_level = zero_level
         self.photon_sensitivity = photon_sensitivity
         self.encoded_dtype = encoded_dtype
         self.decoded_dtype = decoded_dtype
 
     def encode(self, buf: np.ndarray) -> np.ndarray:
         lookup_table = make_anscombe_lookup(
-            self.photon_sensitivity, 
+            self.photon_sensitivity,
             output_type=self.encoded_dtype,
             zero_level=self.zero_level,
         )
         encoded = lookup(buf, lookup_table)
         shape = [encoded.ndim] + list(encoded.shape)
-        shape = np.array(shape, dtype='uint8')
+        shape = np.array(shape, dtype="uint32")
         return shape.tobytes() + encoded.astype(self.encoded_dtype).tobytes()
 
     def decode(self, buf: bytes, out=None) -> np.ndarray:
         lookup_table = make_anscombe_lookup(
-            self.photon_sensitivity, 
+            self.photon_sensitivity,
             output_type=self.encoded_dtype,
             zero_level=self.zero_level,
         )
         inverse_table = make_inverse_lookup(
-            lookup_table,
-            output_type=self.decoded_dtype
+            lookup_table, output_type=self.decoded_dtype
         )
-        ndims = int(buf[0])
-        shape = [int(_) for _ in buf[1:ndims+1]]
-        decoded = np.frombuffer(buf[ndims+1:], dtype=self.encoded_dtype).reshape(shape)
+        ndims = np.frombuffer(buf[:4], "uint32")[0]
+        shape = np.frombuffer(buf[4 : 4 * (ndims + 1)], "uint32")
+        decoded = np.frombuffer(
+            buf[(ndims + 1) * 4 :], dtype=self.encoded_dtype
+        ).reshape(shape)
         return lookup(decoded, inverse_table).astype(self.decoded_dtype)
 
 
-numcodecs.register_codec(PoissonCodec)
+numcodecs.register_codec(AnscombeCodec)

src/poisson_numcodecs/estimate.py → src/anscombe_numcodecs/estimate.py

@@ -17,12 +17,12 @@ def _longest_run(bool_array: np.ndarray) -> slice:
     return slice(on[i], off[i])
 
 
-def compute_sensitivity(movie: np.array, count_weight_gamma: float=0.2) -> dict:
+def compute_sensitivity(movie: np.array, count_weight_gamma: float = 0.2) -> dict:
     """Calculate photon sensitivity
 
     Args:
         movie (np.array):  A movie in the format (time, height, width).
-        count_weight_gamma: 0.00001=weigh each intensity level equally, 
+        count_weight_gamma: 0.00001=weigh each intensity level equally,
             1.0=weigh each intensity in proportion to pixel counts.
 
     Returns:
@@ -42,12 +42,14 @@ def compute_sensitivity(movie: np.array, count_weight_gamma: float=0.2) -> dict:
     intensity = (movie[:-1, :, :] + movie[1:, :, :] + 1) // 2
     difference = movie[:-1, :, :].astype(np.float32) - movie[1:, :, :]
 
-    select = intensity > 0   # discard non-positive values
+    select = intensity > 0  # discard non-positive values
     intensity = intensity[select]
     difference = difference[select]
 
     counts = np.bincount(intensity.flatten())
-    bins = _longest_run(counts > 0.01 * counts.mean())  # consider only bins with at least 1% of mean counts 
+    bins = _longest_run(
+        counts > 0.01 * counts.mean()
+    )  # consider only bins with at least 1% of mean counts
     bins = slice(max(bins.stop * 3 // 100, bins.start), bins.stop)
     assert (
         bins.stop - bins.start > 100
@@ -63,9 +65,9 @@ def compute_sensitivity(movie: np.array, count_weight_gamma: float=0.2) -> dict:
         / counts
     )
     model = Regressor()
-    model.fit(np.c_[bins], variance, counts ** count_weight_gamma)
+    model.fit(np.c_[bins], variance, counts**count_weight_gamma)
     sensitivity = model.coef_[0]
-    zero_level = - model.intercept_ / model.coef_[0]
+    zero_level = -model.intercept_ / model.coef_[0]
 
     return dict(
         model=model,

tests/test_poisson_codec.py → tests/test_codec.py

@@ -1,4 +1,4 @@
-from poisson_numcodecs import PoissonCodec
+from anscombe_numcodecs import AnscombeCodec
 import numpy as np
 import pytest
 
@@ -24,7 +24,7 @@ def test_data(dtype="int16"):
     return [test2d, test2d_long]
 
 def test_poisson_encode_decode(test_data):
-    poisson_codec = PoissonCodec(
+    poisson_codec = AnscombeCodec(
         zero_level=0,
         photon_sensitivity=sensitivity,
         encoded_dtype='uint8',