Updated structure #13

jeromelecoq · 2024-03-02T04:08:16Z

This PR follow up on your great Hackathonwork and introduces a few useful things:

2 calibrating classes that share the same interface to facilitate comparison
A standardization of variable names throughout
2 notebooks to evaluate both compression approaches
Tests that compare the calibration directly to make sure we are starting from the same place
A modification of the encoding format to better match Zarr philosophy
Simplification of the code throughout

Fix buffer decoding to convert buffer to nd.array

Fix decoding formula

THis function was in the core and not necessary

Move to beta convention

Try to use the full unbiased backward form

typo

dimitri-yatsenko

I am marking these as a comment for now.

dimitri-yatsenko · 2024-03-02T07:52:42Z

src/poisson_numcodecs/poisson.py

-        return decoded.astype(self.decoded_dtype)
+        if self.use_lookup:
+            # produce anscombe lookup_tables
+            input_max = np.iinfo(self.decoded_dtype).max


oof, what if it's int64?

dimitri-yatsenko · 2024-03-02T07:56:41Z

src/poisson_numcodecs/poisson.py


+            # https://en.wikipedia.org/wiki/Anscombe_transform for the forward


We subtract sqrt(⅜) to map zero in the input to zero in the output. Otherwise. 0 will become 1 when rounding.

dimitri-yatsenko · 2024-03-02T08:01:49Z

README.md

@@ -4,12 +4,12 @@

 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 (`dark_signal`) and that the `photon_sensitivity` (the average increase in intensity per photon) is known or can be accurately estimated from the data.


Why dark_signal? "Signal" implies a time series. I would use a term that indicate a scalar value, e.g. "dark_level" or "dark_value".

dimitri-yatsenko · 2024-03-02T08:03:35Z

src/poisson_numcodecs/calibrate.py

            slope = found_fits[i, 0] 
            offset = found_fits[i, 1]


Suggested change

slope = found_fits[i, 0]

offset = found_fits[i, 1]

slope, offset = found_fits[i, :]

dimitri-yatsenko · 2024-03-02T08:04:49Z

src/poisson_numcodecs/poisson.py

+        if self.use_lookup:
+            decoded = self._lookup(dec,  self.inverse_table)
+        else:
+            # We first unapply beta


Suggested change

# We first unapply beta

dimitri-yatsenko · 2024-03-02T08:05:43Z

src/poisson_numcodecs/poisson.py

+            # We convert back to arbitrary pixels
+            dec = dec * self.photon_sensitivity + self.dark_signal
+
+            # We have to go back to integers


The code explains itself.

Suggested change

# We have to go back to integers

# We have to go back to integers

dimitri-yatsenko · 2024-03-02T08:10:59Z

src/poisson_numcodecs/poisson.py

+            # JEROME: I do not understand why we need to subtract np.sqrt(3/8) here
+            # Also, shouldn't te +3/8 be inside the maximum function as xx is a float?
+            forward = 2.0 / self.beta * (np.sqrt(np.maximum(0, xx) + 3/8) - np.sqrt(3/8))
+
+            # JEROME : I think this might be better syntax?
+            forward = np.round(forward).astype(self.encoded_dtype)
+            self.forward_table = forward


We subtracted sqrt(⅜) to map zero in the input to zero in the output.

Suggested change

# JEROME: I do not understand why we need to subtract np.sqrt(3/8) here

# Also, shouldn't te +3/8 be inside the maximum function as xx is a float?

forward = 2.0 / self.beta * (np.sqrt(np.maximum(0, xx) + 3/8) - np.sqrt(3/8))

# JEROME : I think this might be better syntax?

forward = np.round(forward).astype(self.encoded_dtype)

self.forward_table = forward

forward = 2.0 / self.beta * (np.sqrt(np.maximum(0, xx) + 3/8) - np.sqrt(3/8))

self.forward_table = forward.astype(self.encoded_dtype)

As long as the reverse is computed correctly from this forward function, then it does not make a big difference.
We just want to use the range of values starting from zero.

dimitri-yatsenko · 2024-03-02T08:16:29Z

src/poisson_numcodecs/poisson.py

+    def _lookup(self, movie, LUT):
+        """
+        Apply lookup table LUT to input movie
+        """
+        return LUT[np.maximum(0, np.minimum(movie, LUT.size-1))]


Suggested change

def _lookup(self, movie, LUT):

"""

Apply lookup table LUT to input movie

"""

return LUT[np.maximum(0, np.minimum(movie, LUT.size-1))]

def _lookup(self, movie, lookup_table):

"""

Apply lookup table to input movie

"""

return lookup_table[np.maximum(0, np.minimum(movie, lookup_table.size-1))]

dimitri-yatsenko · 2024-03-02T11:29:30Z

src/poisson_numcodecs/calibrate.py

+
+        return [self.photon_sensitivity, self.dark_signal]
+
+    def plot_poisson_curve(self):


In the literature, this is known as the Photon Transfer Curve.

dimitri-yatsenko · 2024-03-02T11:32:48Z

src/poisson_numcodecs/calibrate.py


-class RasterCalibratePhotons():
+class CalibratePhotons():


Do we really need these classes? All the functionality here can be implemented as pure functions with greater readability and maintainability. Class hierarchies require a lot of engineering and should be avoided when a function will do.

"calibrate" is not quite the right word. We calibrate a device. What we are doing here is more like estimate_photon_sensitivity.

dimitri-yatsenko · 2024-03-02T11:44:28Z

src/poisson_numcodecs/calibrate.py

+        self.fitted_pixels_mean = None
+        self.fitted_model = None
+
+    def _longest_run(self, bool_array: np.ndarray) -> slice:


This does not use self. This is a pure function. Does not need to be in a class. You could use @staticmethod but really we don't need these classes here.

dimitri-yatsenko · 2024-03-02T11:51:22Z

src/poisson_numcodecs/calibrate.py

+    def subsample_and_crop_video(self, crop, start_frame=0, end_frame=-1):
+        """Subsample and crop a video, cache results. Also functions as a data_pointer load.
+
+        Args:
+            crop:  A tuple (px_y, px_x) specifying the number of pixels to remove
+            start_frame:  The index of the first desired frame
+            end_frame:  The index of the last desired frame
+
+        Returns:
+            The resultant array.
+        """
+
+        # We first reset the saved data
+        self.mean_image = None
+        self.std_image = None
+        self.photon_sensitivity = None
+        self.dark_signal = None
+
+        _shape = self.data_array_movie.shape
+        px_y_start, px_x_start = crop
+        px_y_end = _shape[1] - px_y_start
+        px_x_end = _shape[2] - px_x_start
+
+        if start_frame == _shape[0] - 1 and (end_frame == -1 or end_frame == _shape[0]):
+            cropped_video = self.data_array_movie[
+                start_frame:_shape[0], px_y_start:px_y_end, px_x_start:px_x_end
+            ]
+        else:
+            cropped_video = self.data_array_movie[
+                start_frame:end_frame, px_y_start:px_y_end, px_x_start:px_x_end
+            ]
+        self.data_array_movie = cropped_video


Why does this method belong in this library?

Suggested change

def subsample_and_crop_video(self, crop, start_frame=0, end_frame=-1):

"""Subsample and crop a video, cache results. Also functions as a data_pointer load.

Args:

crop: A tuple (px_y, px_x) specifying the number of pixels to remove

start_frame: The index of the first desired frame

end_frame: The index of the last desired frame

Returns:

The resultant array.

"""

# We first reset the saved data

self.mean_image = None

self.std_image = None

self.photon_sensitivity = None

self.dark_signal = None

_shape = self.data_array_movie.shape

px_y_start, px_x_start = crop

px_y_end = _shape[1] - px_y_start

px_x_end = _shape[2] - px_x_start

if start_frame == _shape[0] - 1 and (end_frame == -1 or end_frame == _shape[0]):

cropped_video = self.data_array_movie[

start_frame:_shape[0], px_y_start:px_y_end, px_x_start:px_x_end

]

else:

cropped_video = self.data_array_movie[

start_frame:end_frame, px_y_start:px_y_end, px_x_start:px_x_end

]

self.data_array_movie = cropped_video

dimitri-yatsenko · 2024-03-02T12:01:43Z

src/poisson_numcodecs/poisson.py

+                        inverse.size * (inverse[-1] - inverse[-2])/2).astype(self.decoded_dtype)
+            self.inverse_table = inverse
+
+    def _lookup(self, movie, LUT):


This function does not need self.

dimitri-yatsenko · 2024-03-02T12:02:12Z

src/poisson_numcodecs/poisson.py

+        if self.use_lookup:
+            encoded = self._lookup(buf,  self.forward_table)
+        else:
+            # We convert to photons


Suggested change

# We convert to photons

# convert to photons

dimitri-yatsenko · 2024-03-02T12:03:04Z

src/poisson_numcodecs/poisson.py

+
+            # We have to go to integers in a clean way


Here the code explains itself better than the comment.

Suggested change

# We have to go to integers in a clean way

# We have to go to integers in a clean way

dimitri-yatsenko · 2024-03-02T12:04:06Z

src/poisson_numcodecs/poisson.py

+            # https://en.wikipedia.org/wiki/Anscombe_transform for the inverse without bias
+            dec = dec**2 / 4.0 - 1/8
+
+            # We convert back to arbitrary pixels


Suggested change

# We convert back to arbitrary pixels

# restore original grayscale

dimitri-yatsenko · 2024-03-02T12:04:30Z

tests/test_poisson_calibrate.py

-            [photon_gain, photon_offset]=calibrator.get_photon_gain_parameters(perc_min=0, perc_max=100)
-            print(photon_gain, photon_offset)
+            [photon_sensitivity, dark_signal]=calibrator.get_photon_sensitivity_parameters(perc_min=0, perc_max=100)
+            print(photon_sensitivity, dark_signal)

            # We check that the gain and offset are within X% of the true value


Suggested change

# We check that the gain and offset are within X% of the true value

# check that the gain and offset are within X% of the true value

dimitri-yatsenko · 2024-03-02T12:41:53Z

src/poisson_numcodecs/calibrate.py

-        return photon_flux
-
-    def get_photon_gain_parameters(self, max_pixel_range=2**15, n_groups=1, perc_min=3, perc_max=90):
+    def get_photon_sensitivity_parameters(self, max_pixel_range=2**15, n_groups=1, perc_min=3, perc_max=90):


max_pixel_range you probably mean "value", not "range". It should be 2**15-1 == 32767 == 0x7FFF. 2**15 overflows in int16.

dimitri-yatsenko · 2024-03-02T12:44:37Z

src/poisson_numcodecs/poisson.py

+
+    def _lookup(self, movie, LUT):
+        """
+        Apply lookup table LUT to input movie


Suggested change

Apply lookup table LUT to input movie

Apply lookup table to movie

jeromelecoq added 27 commits February 27, 2024 13:07

Update Poisson.py

ba7df37

Fix buffer decoding to convert buffer to nd.array

Update Poisson.py

0507adb

Fix decoding formula

Update Poisson.py

5946cb5

Update Poisson.py

50f18fc

THis function was in the core and not necessary

Update Poisson.py

f1056ee

Move to beta convention

Update README.md

ca2e394

Update test_poisson_codec.py

fe51a82

Update Poisson.py

8ca9a01

Update Poisson.py

08be71f

Update Poisson.py

9afee3f

Stick to Anscombe

dccfdca

Try to use the full unbiased backward form

Update Poisson.py

016a289

typo

Update Poisson.py

6ff5c51

Update Poisson.py

822e384

Fix to use without bytes

c483810

Edit poisson for simplicity

59e6c56

update naming

11bf76e

Fix typo

5129109

Add non lookup test

14ca150

Add framework with both calibration

6192ca4

Update notebooks

0c26c26

Add 2 notebooks

5c9ec50

Merge branch 'main' into datajoint/main

a47cc05

Update raster notebook

00dc316

rename file

cff4010

Rename

6ef1920

fix path

3b27416

dimitri-yatsenko reviewed Mar 2, 2024

View reviewed changes

Update notebooks with comparison with errors

745b3b5

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Updated structure #13

Updated structure #13

jeromelecoq commented Mar 2, 2024

dimitri-yatsenko left a comment

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading


		# https://en.wikipedia.org/wiki/Anscombe_transform for the forward

	slope = found_fits[i, 0]
	offset = found_fits[i, 1]
	slope, offset = found_fits[i, :]

	# We have to go back to integers
	# We have to go back to integers


		return [self.photon_sensitivity, self.dark_signal]

		def plot_poisson_curve(self):


	# We have to go to integers in a clean way

	# We have to go to integers in a clean way

	# We convert back to arbitrary pixels
	# restore original grayscale

	# We check that the gain and offset are within X% of the true value
	# check that the gain and offset are within X% of the true value

	Apply lookup table LUT to input movie
	Apply lookup table to movie


		class RasterCalibratePhotons():
		class CalibratePhotons():

Updated structure #13

Are you sure you want to change the base?

Updated structure #13

Conversation

jeromelecoq commented Mar 2, 2024

dimitri-yatsenko left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 • edited Loading

Choose a reason for hiding this comment

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading

dimitri-yatsenko Mar 2, 2024 •

edited

Loading