docs: update docstrings and docs example output image

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_code_reference.rst

@@ -23,21 +23,34 @@ Fourier transform methods
 =========================
 
 .. _sec_code_fourier_numpy:
-
 Numpy
 -----
 .. automodule:: qpretrieve.fourier.ff_numpy
    :members:
    :inherited-members:
 
 .. _sec_code_fourier_pyfftw:
-
 PyFFTW
 ------
 .. automodule:: qpretrieve.fourier.ff_pyfftw
    :members:
    :inherited-members:
 
+.. _sec_code_fourier_scipy:
+Scipy
+------
+.. automodule:: qpretrieve.fourier.ff_scipy
+   :members:
+   :inherited-members:
+
+.. _sec_code_fourier_cupy:
+Cupy
+----
+.. automodule:: qpretrieve.fourier.ff_cupy
+   :members:
+   :inherited-members:
+
+
 .. _sec_code_ifer:
 
 Interference image analysis

docs/sec_examples.rst

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

examples/fft_options.py

@@ -1,51 +1,73 @@
-"""Fourier Transform options available
+"""Fourier Transform interfaces available
 
 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.
+- CuPy using CUDA can be very fast, but is currently limited because we are
+  transferring one image at a time to the GPU.
+
 """
+import time
 import matplotlib.pylab as plt
 import numpy as np
 import qpretrieve
-from skimage.restoration import unwrap_phase
 
 # load the experimental data
 edata = np.load("./data/hologram_cell.npz")
 
 # get the available fft interfaces
 interfaces_available = qpretrieve.fourier.get_available_interfaces()
 
-prange = (-1, 5)
-frange = (0, 12)
-
-results = {}
+n_transforms = 100
 
+# one transform
+results_1 = {}
 for fft_interface in interfaces_available:
+    t0 = time.time()
     holo = qpretrieve.OffAxisHologram(data=edata["data"],
                                       fft_interface=fft_interface)
-    holo.run_pipeline(filter_name="disk", filter_size=1/2)
+    holo.run_pipeline(filter_name="disk", filter_size=1 / 2)
     bg = qpretrieve.OffAxisHologram(data=edata["bg_data"])
     bg.process_like(holo)
-    phase = unwrap_phase(holo.phase - bg.phase)
-    mask = np.log(1 + np.abs(holo.fft_filtered))
-    results[fft_interface.__name__] = mask, phase
+    t1 = time.time()
+    results_1[fft_interface.__name__] = t1 - t0
+num_interfaces = len(results_1)
 
-num_filters = len(results)
+# multiple transforms (should see speed increase for PyFFTW)
+results = {}
+for fft_interface in interfaces_available:
+    t0 = time.time()
+    for _ in range(n_transforms):
+        holo = qpretrieve.OffAxisHologram(data=edata["data"],
+                                          fft_interface=fft_interface)
+        holo.run_pipeline(filter_name="disk", filter_size=1 / 2)
+        bg = qpretrieve.OffAxisHologram(data=edata["bg_data"])
+        bg.process_like(holo)
+    t1 = time.time()
+    results[fft_interface.__name__] = t1 - t0
+num_interfaces = len(results)
 
-# plot the properties of `qpi`
-fig = plt.figure(figsize=(8, 22))
+fft_interfaces = list(results.keys())
+speed_1 = list(results_1.values())
+speed = list(results.values())
 
-for row, name in enumerate(results):
-    ax1 = plt.subplot(num_filters, 2, 2*row+1)
-    ax1.set_title(name, loc="left")
-    ax1.imshow(results[name][0], vmin=frange[0], vmax=frange[1])
+fig, axes = plt.subplots(1, 2, figsize=(8, 5))
+ax1, ax2 = axes
+labels = [fftstr[9:] for fftstr in fft_interfaces]
 
-    ax2 = plt.subplot(num_filters, 2, 2*row+2)
-    map2 = ax2.imshow(results[name][1], cmap="coolwarm",
-                      vmin=prange[0], vmax=prange[1])
-    plt.colorbar(map2, ax=ax2, fraction=.046, pad=0.02, label="phase [rad]")
+ax1.bar(range(num_interfaces), height=speed_1, color='lightseagreen')
+ax1.set_xticks(range(num_interfaces), labels=labels,
+               rotation=45)
+ax1.set_ylabel("Speed (s)")
+ax1.set_title("1 Transform")
 
-    ax1.axis("off")
-    ax2.axis("off")
+ax2.bar(range(num_interfaces), height=speed, color='lightseagreen')
+ax2.set_xticks(range(num_interfaces), labels=labels,
+               rotation=45)
+ax2.set_ylabel("Speed (s)")
+ax2.set_title(f"{n_transforms} Transforms")
 
+plt.suptitle("Speed of FFT Interfaces")
 plt.tight_layout()
 plt.show()

qpretrieve/filter.py

@@ -38,7 +38,7 @@ def get_filter_array(filter_name, filter_size, freq_pos, fft_shape):
         and must be between 0 and `max(fft_shape)/2`
     freq_pos: tuple of floats
         The position of the filter in frequency coordinates as
-        returned by :func:`nunpy.fft.fftfreq`.
+        returned by :func:`numpy.fft.fftfreq`.
     fft_shape: tuple of int
         The shape of the Fourier transformed image for which the
         filter will be applied. The shape must be squared (two

qpretrieve/fourier/__init__.py

@@ -9,6 +9,11 @@
 except ImportError:
     FFTFilterPyFFTW = None
 
+try:
+    from .ff_cupy import FFTFilterCupy
+except ImportError:
+    FFTFilterCupy = None
+
 PREFERRED_INTERFACE = None
 
 
@@ -18,6 +23,7 @@ def get_available_interfaces():
         FFTFilterPyFFTW,
         FFTFilterNumpy,
         FFTFilterScipy,
+        FFTFilterCupy,
     ]
     interfaces_available = []
     for interface in interfaces:

qpretrieve/fourier/base.py

@@ -178,7 +178,7 @@ def filter(self, filter_name: str, filter_size: float,
             and must be between 0 and `max(fft_shape)/2`
         freq_pos: tuple of floats
             The position of the filter in frequency coordinates as
-            returned by :func:`nunpy.fft.fftfreq`.
+            returned by :func:`numpy.fft.fftfreq`.
         scale_to_filter: bool or float
             Crop the image in Fourier space after applying the filter,
             effectively removing surplus (zero-padding) data and

qpretrieve/fourier/ff_cupy.py

@@ -0,0 +1,40 @@
+import scipy as sp
+import cupy as cp
+import cupyx.scipy.fft as cufft
+
+from .base import FFTFilter
+
+
+class FFTFilterCupy(FFTFilter):
+    """Wraps the cupy Fourier transform and uses it via the scipy backend
+    """
+    is_available = True
+    # sp.fft.set_backend(cufft)
+
+    def _init_fft(self, data):
+        """Perform initial Fourier transform of the input data
+
+        Parameters
+        ----------
+        data: 2d real-valued np.ndarray
+            Input field to be refocused
+
+        Returns
+        -------
+        fft_fdata: 2d complex-valued ndarray
+            Fourier transform `data`
+        """
+        data_gpu = cp.asarray(data)
+        # likely an inefficiency here, could use `set_global_backend`
+        with sp.fft.set_backend(cufft):
+            fft_gpu = sp.fft.fft2(data_gpu)
+        fft_cpu = fft_gpu.get()
+        return fft_cpu
+
+    def _ifft(self, data):
+        """Perform inverse Fourier transform"""
+        data_gpu = cp.asarray(data)
+        with sp.fft.set_backend(cufft):
+            ifft_gpu = sp.fft.ifft2(data_gpu)
+        ifft_cpu = ifft_gpu.get()
+        return ifft_cpu

qpretrieve/interfere/base.py

@@ -22,6 +22,13 @@ def __init__(self, data, fft_interface: FFTFilter = None,
         """
         Parameters
         ----------
+        fft_interface: FFTFilter
+            A Fourier transform interface.
+            See :func:`qpretrieve.fourier.get_available_interfaces`
+            to get a list of implemented interfaces.
+            Default is None, which will use
+            :func:`qpretrieve.fourier.get_best_interface`. This is in line
+            with old behaviour.
         subtract_mean: bool
             If True, remove the mean of the hologram before performing
             the Fourier transform. This setting is recommended as it

setup.py

@@ -2,7 +2,6 @@
 from setuptools import setup, find_packages
 import sys
 
-
 author = u"Paul Müller"
 authors = [author]
 description = 'library for phase retrieval from holograms'
@@ -27,8 +26,13 @@
         "numpy>=1.9.0",
         "scikit-image>=0.11.0",
         "scipy>=0.18.0",
-        ],
-    extras_require={"FFTW": "pyfftw>=0.12.0"},
+    ],
+    extras_require={
+        "FFTW": "pyfftw>=0.12.0",
+        # manually install 'cupy-cuda11x' if you have older CUDA.
+        # See https://cupy.dev/
+        "CUPY": "cupy-cuda12x",
+    },
     python_requires='>=3.10, <4',
     keywords=["digital holographic microscopy",
               "optics",
@@ -41,6 +45,6 @@
         'Operating System :: OS Independent',
         'Programming Language :: Python :: 3',
         'Intended Audience :: Science/Research'
-                 ],
+    ],
     platforms=['ALL'],
-    )
+)