Move from 2d to 3d array operations

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,55 @@
+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 in 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
+-------------------------------
+
+If you give either a RGB or RGBA array layout as input, then the first
+channel is taken as the image to process. In other words, it is assumed that
+all channels contain the same information, so the first channel is used.
+
+If you use the `oah.get_array_with_input_layout("phase")` method for
+the RGBA array layout, then the alpha (A) channel will be an array of ones.
+
+3D RGB/RGBA array layouts, such as (50, 256, 256, 3), are not allowed (yet).
+
+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_array_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_array_with_input_layout("field")
+	phase_2d = oah.get_array_with_input_layout("phase")
+	amplitude_2d = oah.get_array_with_input_layout("amplitude")

docs/sec_basic_use.rst

@@ -0,0 +1,75 @@
+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 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_array_with_input_layout("field")
+
+	# 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 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_examples.rst

@@ -12,3 +12,7 @@ Examples
 .. fancy_include:: filter_visualization.py
 
 .. fancy_include:: fourier_scale.py
+
+.. fancy_include:: fft_options.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,86 @@
+"""Fourier Transform speed benchmarks
+
+This example visualizes the speed for different batch sizes for
+the available FFT Filters. The y-axis shows the average speed of a single
+FFT for the corresponding batch size.
+
+- 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.
+
+"""
+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
+
+# 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
+fft_interfaces = qpretrieve.fourier.get_available_interfaces()
+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()]
+    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 = ["lightseagreen", "darkmagenta"]
+
+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("Fourier transform batch size")
+ax1.set_ylabel("Time for single FFT [Time / batch size] (s)")
+ax1.legend(loc='upper right')
+
+plt.suptitle("Batch processing time for FFT Filters")
+plt.tight_layout()
+# plt.show()
+plt.savefig("fft_batch_speeds.png", dpi=150)

examples/fft_options.py

@@ -0,0 +1,122 @@
+"""Available Fourier Transform backends (`FFTFilter`)
+
+This example visualizes the different backends and packages available to the
+user for performing Fourier transforms.
+
+- PyFFTW is initially slow, but over many FFTs is very quick.
+
+"""
+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
+
+# load the experimental data
+edata = np.load("./data/hologram_cell.npz")
+
+# get the available fft interfaces
+interfaces_available = qpretrieve.fourier.get_available_interfaces()
+num_interfaces = len(interfaces_available)
+fft_interfaces = [interf.__name__ for interf in interfaces_available]
+
+n_transforms = 200
+subtract_mean = True
+padding = True
+filter_name = "disk"
+filter_size = 1 / 2
+
+# 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)
+
+print("Running single transform...")
+# one transform
+results_1 = {}
+for fft_interface in interfaces_available:
+    t0 = time.time()
+    holo = qpretrieve.OffAxisHologram(data=data_2d,
+                                      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_2d_bg)
+    bg.process_like(holo)
+    t1 = time.time()
+    results_1[fft_interface.__name__] = t1 - t0
+num_interfaces = len(results_1)
+
+# multiple transforms repeated in 2D
+print(f"Running {n_transforms} transforms in a loop...")
+results_2d = {}
+for fft_interface in interfaces_available:
+
+    t0 = time.time()
+    for _ in range(n_transforms):
+        assert data_3d_prep.shape == data_3d_bg_prep.shape == (
+            1, edata["data"].shape[0], edata["data"].shape[1])
+
+        holo = qpretrieve.OffAxisHologram(data=data_3d_prep,
+                                          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_prep)
+        bg.process_like(holo)
+    t1 = time.time()
+    results_2d[fft_interface.__name__] = t1 - t0
+
+# multiple transforms in 3D
+print(f"Running {n_transforms} transforms at once...")
+results_3d = {}
+for fft_interface in interfaces_available:
+    # create batch
+    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_3d[fft_interface.__name__] = t1 - t0
+
+speed_1 = list(results_1.values())
+speed_2d = list(results_2d.values())
+speed_3d = list(results_3d.values())
+
+fig, axes = plt.subplots(1, 3, figsize=(12, 5))
+ax1, ax2, ax3 = axes
+labels = [fftstr[9:] for fftstr in fft_interfaces]
+colors = ["lightseagreen", "darkmagenta"]
+
+ax1.bar(range(num_interfaces), height=speed_1, color=colors, edgecolor='k')
+ax1.set_xticks(range(num_interfaces), labels=labels)
+ax1.set_ylabel("Time (s)")
+ax1.set_title("1 Transform")
+
+ax2.bar(range(num_interfaces), height=speed_2d, color=colors, edgecolor='k')
+ax2.set_xticks(range(num_interfaces), labels=labels)
+ax2.set_ylabel("Time (s)")
+ax2.set_title(f"{n_transforms} Transforms (one-by-one)")
+ax2.set_ylim(0, 4.5)
+
+ax3.bar(range(num_interfaces), height=speed_3d, color=colors, edgecolor='k')
+ax3.set_xticks(range(num_interfaces), labels=labels)
+ax3.set_ylabel("Time (s)")
+ax3.set_title(f"{n_transforms} Transforms (single 3D batch)")
+ax3.set_ylim(0, 4.5)
+
+plt.suptitle("Processing time for FFT Filters")
+plt.tight_layout()
+# plt.show()
+plt.savefig("fft_options.png", dpi=150)

examples/requirements.txt

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