Move from 2d to 3d array operations

CHANGELOG

@@ -1,3 +1,5 @@
+0.4.0
+ - feat: move from 2d to 3d array operations (#12)
 0.3.2
  - maintenance release
 0.3.1

docs/README.md

@@ -6,6 +6,7 @@ To install the requirements for building the documentation, run
 
 To compile the documentation, run
 
+    cd docs
     sphinx-build . _build
 
 

docs/requirements.txt

@@ -1,4 +1,3 @@
-sphinx==4.3.0
-sphinxcontrib.bibtex>=2.0
-sphinx_rtd_theme==1.0
-
+sphinx
+sphinxcontrib.bibtex
+sphinx_rtd_theme

docs/sec_array_layout.rst

@@ -0,0 +1,62 @@
+Data Array Layouts
+==================
+
+.. _sec_doc_array_layout:
+
+Since version 0.4.0, `qpretrieve` accepts 3D (z,y,x) arrays as input.
+Additionally, it **always** returns data as the 3D array layout.
+
+We use the term "*array layout*" to define the different ways to represent data.
+The currently accepted array layouts are: 2D, RGB, RGBA, 3D.
+
+Summary of allowed Array Layouts::
+
+    Input         ->  Output
+    2D   (y,x)    ->  3D (1,y,x)
+    RGB  (y,x,3)  ->  3D (1,y,x)
+    RGBA (y,x,4)  ->  3D (1,y,x)
+    3D   (z,y,x)  ->  3D (z,y,x)
+
+
+Notes on RGB/RGBA Array Layouts
+-------------------------------
+
+**Inputting RGB(A)**: See the Notes section of
+:func:`.convert_data_to_3d_array_layout` for extra information on RGB(A)
+array layouts.
+
+**Outputting RGB(A)**: See the Notes section of
+:meth:`.OffAxisHologram.get_data_with_input_layout` or
+:meth:`.QLSInterferogram.get_data_with_input_layout` for information on
+outputting of RGB(A) array layouts.
+
+
+Converting to and from Array Layouts
+------------------------------------
+
+`qpretrieve` will automatically handle the above allowed array layouts.
+In other words, if you provide any 2D, RGB, RGBA, or 3D data as input to
+:class:`.OffAxisHologram` or :class:`.QLSInterferogram`
+the class will handle everything for you.
+
+However, if you want to have your processed data in the same array layout as when
+you inputted it, then you can use the convenience method
+:meth:`get_data_with_input_layout` to do just that. For example, if
+your input data was a 2D array, you can get the processed field, phase,
+amplitude etc like so:
+
+.. code-block:: python
+
+	# 2D data inputted
+	oah = qpretrieve.OffAxisHologram(hologram_2d)
+	# do some processing
+	...
+	# get your data as a 2D array layout
+	field_2d = oah.get_data_with_input_layout("field")
+	phase_2d = oah.get_data_with_input_layout("phase")
+	amplitude_2d = oah.get_data_with_input_layout("amplitude")
+
+	# you can also use the class attributes
+	field_2d = oah.get_data_with_input_layout(oah.field)
+	phase_2d = oah.get_data_with_input_layout(oah.phase)
+	amplitude_2d = oah.get_data_with_input_layout(oah.amplitude)

docs/sec_basic_use.rst

@@ -0,0 +1,78 @@
+Basic Use of qpretrieve
+=======================
+
+To run this code locally, you can download the example dataset from the
+`qpretrieve repository <https://github.com/RI-imaging/qpretrieve/blob/main/examples/data/hologram_cell.npz>`_
+
+
+.. code-block:: python
+
+	import qpretrieve
+
+	# load your experimental data (data from the qpretrieve repository)
+	edata = np.load("qpretrieve/examples/data/hologram_cell.npz")
+	data = edata["data"]
+	data_bg = edata["bg_data"]
+
+	# create an off-axis hologram object, to process the holography data
+	oah = qpretrieve.OffAxisHologram(data)
+	oah.run_pipeline(filter_name=filter_name, filter_size=filter_size)
+	# process background hologram
+	oah_bg = qpretrieve.OffAxisHologram(data=data_2d_bg)
+	oah_bg.process_like(oah)
+
+	print(f"Original Hologram shape: {data.shape}")
+	print(f"Processed Field shape: {oah.field.shape}")
+
+	# Now you can look at the phase data
+	phase_corrected = oah.phase - oah_bg.phase
+
+
+
+How qpretrieve now handles data
+-------------------------------
+
+In older versions of qpretrieve, only a single image could be processed at a time.
+Since version (0.4.0), it accepts 3D array layouts, and always returns
+the data as 3D, regardless of the input array layout.
+
+Click here for more :ref:`information on array layouts <sec_doc_array_layout>`.
+There are speed benchmarks given for different inputs in the
+:ref:`examples <sec_examples>` section.
+
+
+New version example code and output:
+....................................
+
+.. code-block:: python
+
+	import numpy as np
+	import qpretrieve
+
+	hologram = np.ones(shape=(256, 256))
+	oah = qpretrieve.OffAxisHologram(hologram)
+	oah.run_pipeline()
+	assert oah.field.shape == (1, 256, 256)  # <- now a 3D array is returned
+	# if you want the original array layout (2d)
+	field_2d = oah.get_data_with_input_layout("field")
+	assert field_2d.shape == (256, 256)
+
+	# this means you can input 3D arrays
+	hologram_3d = np.ones(shape=(50, 256, 256))
+	oah = qpretrieve.OffAxisHologram(hologram_3d)
+	oah.run_pipeline()
+	assert oah.field.shape == (50, 256, 256)  # <- always a 3D array
+
+
+Old version example code and output:
+....................................
+
+.. code-block:: python
+
+	import numpy as np
+	import qpretrieve  # versions older than 0.4.0
+
+	hologram = np.ones(shape=(256, 256))
+	oah = qpretrieve.OffAxisHologram(hologram)
+	oah.run_pipeline()
+	assert oah.field.shape == hologram.shape  # <- old version returned a 2D array

docs/sec_code_reference.rst

@@ -60,3 +60,12 @@ Quadriwave lateral shearing interferometry (QLSI)
 .. automodule:: qpretrieve.interfere.if_qlsi
    :members:
    :inherited-members:
+
+Data Array Layout
+=================
+
+.. _sec_code_array_layout:
+
+.. automodule:: qpretrieve.data_array_layout
+   :members:
+   :inherited-members:

docs/sec_examples.rst

@@ -12,3 +12,5 @@ Examples
 .. fancy_include:: filter_visualization.py
 
 .. fancy_include:: fourier_scale.py
+
+.. fancy_include:: fft_batch_speeds.py

docs/sec_getting_started.rst

@@ -6,4 +6,6 @@ Getting started
   :maxdepth: 2
 
   sec_installation
+  sec_basic_use
+  sec_array_layout
   sec_userapi

docs/sec_userapi.rst

@@ -21,4 +21,4 @@ Then your analysis could be as simple as
 With ``dhm``, an instance of :class:`.OffAxisHologram`, you now have full access
 to all intermediate computation results. You can pass additional keyword
 arguments during instantiation or pass them to
-:func:`.OffAxisHologram.run_pipeline`.
+:meth:`.OffAxisHologram.run_pipeline`.

examples/fft_batch_speeds.py

@@ -0,0 +1,91 @@
+"""Fourier Transform speed benchmarks for OAH
+
+This example visualizes the speed for different batch sizes for
+the available FFT Filters. The y-axis shows the average speed of a pipeline
+run for the Off-Axis Hologram class :class:`.OffAxisHologram`, including
+background data processing. Therefore, four FFTs are run per pipeline.
+
+- Optimum batch size is between 64 and 256 for 256x256pix imgs (incl padding),
+  but will be limited by your computer's RAM.
+- Here, batch size is the size of the 3D stack in z.
+- Note that each pipeline runs 4 FFTs. For example, batch 8 runs 8*4=32 FFTs.
+
+"""
+import time
+import matplotlib.pylab as plt
+import numpy as np
+import qpretrieve
+from qpretrieve.data_array_layout import convert_data_to_3d_array_layout
+from qpretrieve.fourier import FFTFilterNumpy, FFTFilterPyFFTW
+
+# load the experimental data
+edata = np.load("./data/hologram_cell.npz")
+
+n_transforms_list = [8, 16, 32, 64, 128, 256]
+subtract_mean = True
+padding = True
+# we take the PyFFTW speeds from the second run
+fft_interfaces = [FFTFilterNumpy, FFTFilterPyFFTW, FFTFilterPyFFTW]
+filter_name = "disk"
+filter_size = 1 / 2
+speed_norms = {}
+
+# load and prep the data
+data_2d = edata["data"].copy()
+data_2d_bg = edata["bg_data"].copy()
+data_3d_prep, _ = convert_data_to_3d_array_layout(data_2d)
+data_3d_bg_prep, _ = convert_data_to_3d_array_layout(data_2d_bg)
+
+for fft_interface in fft_interfaces:
+    results = {}
+    for n_transforms in n_transforms_list:
+        print(f"Running {n_transforms} transforms with "
+              f"{fft_interface.__name__}")
+
+        # create batches
+        data_3d = np.repeat(data_3d_prep, repeats=n_transforms, axis=0)
+        data_3d_bg = np.repeat(data_3d_bg_prep, repeats=n_transforms, axis=0)
+
+        assert data_3d.shape == data_3d_bg.shape == (n_transforms,
+                                                     edata["data"].shape[0],
+                                                     edata["data"].shape[1])
+
+        t0 = time.time()
+        holo = qpretrieve.OffAxisHologram(data=data_3d,
+                                          fft_interface=fft_interface,
+                                          subtract_mean=subtract_mean,
+                                          padding=padding)
+        holo.run_pipeline(filter_name=filter_name, filter_size=filter_size)
+        bg = qpretrieve.OffAxisHologram(data=data_3d_bg)
+        bg.process_like(holo)
+        t1 = time.time()
+        results[n_transforms] = t1 - t0
+
+    speed_norm = [timing / b_size for b_size, timing in results.items()]
+    # the initial PyFFTW run (incl wisdom calc is overwritten here)
+    speed_norms[fft_interface.__name__] = speed_norm
+
+# setup plot
+fig, axes = plt.subplots(1, 1, figsize=(8, 5))
+ax1 = axes
+width = 0.25  # the width of the bars
+multiplier = 0
+x_pos = np.arange(len(n_transforms_list))
+colors = ["darkmagenta", "lightseagreen"]
+
+for (name, speed), color in zip(speed_norms.items(), colors):
+    offset = width * multiplier
+    ax1.bar(x_pos + offset, speed, width, label=name,
+            color=color, edgecolor='k')
+    multiplier += 1
+
+ax1.set_xticks(x_pos + (width / 2), labels=n_transforms_list)
+ax1.set_xlabel("Input hologram batch size")
+ax1.set_ylabel("OAH processing time [Time / batch size] (s)")
+ax1.legend(loc='upper right', fontsize="large")
+
+plt.suptitle("Batch processing time for Off-Axis Hologram\n"
+             "(data+bg_data)")
+plt.tight_layout()
+# plt.show()
+plt.savefig("fft_batch_speeds.png", dpi=150)

examples/requirements.txt

@@ -0,0 +1 @@
+matplotlib