Feature/monogenic (#38)

* Enhance documentation for version management in conf.py and .readthedocs.yaml

- Updated conf.py to clarify version retrieval from source and its implications for stable documentation builds.
- Added detailed comments in .readthedocs.yaml regarding dynamic version generation and stable documentation configurations.

* Update pyproject.toml to support Python 3.13 and 3.14, and enhance pytest configuration

- Added support for Python versions 3.13 and 3.14 in the classifiers section.
- Updated filterwarnings and markers sections for improved pytest configuration consistency.

* Remove deprecated Python 3.14t classifier from pyproject.toml to streamline version support.

* Refactor UDCT and tests to remove Meyer wavelet support

- Removed Meyer wavelet mode from UDCT class and associated tests, simplifying the implementation.
- Updated documentation to reflect the removal of Meyer wavelet references.
- Adjusted tests to ensure compatibility with the new structure, focusing on curvelet and wavelet modes only.
- Enhanced error handling and validation for remaining modes to improve robustness.

* Refactor curvelet vs wavelet example to remove Meyer wavelet references

- Eliminated all code and print statements related to Meyer wavelet, simplifying the example.
- Updated visualizations to focus solely on curvelet and wavelet transforms, enhancing clarity.
- Adjusted subplot configurations to ensure proper display of remaining windows and coefficients.
- Improved documentation to reflect the changes and clarify the example's focus on curvelet and wavelet comparisons.

* Enhance monogenic curvelet transform implementation

- Introduced new functions for processing monogenic transforms, including `_process_wedge_monogenic` and `_apply_forward_transform_monogenic`, to handle scalar and Riesz components.
- Updated the UDCT class to include `forward_monogenic` and `backward_monogenic` methods for applying and reconstructing monogenic coefficients.
- Enhanced type annotations and documentation for clarity, ensuring compatibility with real-valued inputs and providing detailed parameter descriptions.
- Refactored existing typing definitions to include new type aliases for improved type safety and readability.

* Refactor UDCT and MeyerWavelet classes for improved functionality and error handling

- Enhanced the UDCT class with a new _CoefficientsList class to manage highpass bands, ensuring thread safety and preventing data corruption during transformations.
- Updated forward and backward methods to validate highpass band presence and improve error handling.
- Implemented validation for odd dimensions in the MeyerWavelet class, raising appropriate errors for invalid inputs.
- Enhanced dtype handling in the MeyerWavelet class for upsampled arrays, ensuring correct preservation of data types.
- Added comprehensive tests to cover new features and validate error paths, improving overall test coverage and robustness.

* Enhance UDCT and MeyerWavelet classes with improved error handling and validation

- Introduced a new _CoefficientsList class in the UDCT to manage highpass bands, ensuring thread safety and preventing data corruption during transformations.
- Updated forward and backward methods to validate highpass band presence and improve error handling for incomplete coefficients.
- Implemented validation for odd dimensions in the MeyerWavelet class, raising appropriate errors for invalid inputs.
- Enhanced dtype handling in the MeyerWavelet class for upsampled arrays, ensuring correct preservation of data types.
- Added comprehensive tests to cover new features and validate error paths, improving overall test coverage and robustness.

* Enhance monogenic curvelet visualization and consistency checks

- Updated the monogenic curvelet example to improve clarity in mathematical definitions and visualizations.
- Introduced a smaller window overlap parameter in the UDCT transform for better reconstruction accuracy.
- Added comprehensive round-trip consistency checks to validate the monogenic transform against direct computations.
- Enhanced visual comparisons of scalar and Riesz components, including detailed statistics on differences and ratios.
- Improved frequency domain analysis to facilitate understanding of discrepancies between direct and round-trip transformations.

* Update monogenic curvelet example for improved reconstruction accuracy and consistency checks

- Increased window overlap parameter in the UDCT transform to 0.1 for better reconstruction accuracy.
- Simplified the monogenic transformation process by removing redundant small overlap methods.
- Enhanced visualizations and statistics for scalar and Riesz component comparisons.
- Improved print statements for clarity in output comparisons and error analysis.

* Refine tolerances in monogenic tests for improved accuracy

- Adjusted tolerance levels for scalar and Riesz component comparisons in the monogenic round-trip tests to account for numerical precision.
- Removed the skip marker from the test for monogenic matches with backward_monogenic, indicating readiness for execution after addressing known issues.

* Remove extensive round-trip consistency checks and visualizations from monogenic curvelet example

- Eliminated detailed component comparisons and frequency domain analyses to streamline the example.
- Removed print statements and visualizations related to scalar and Riesz component comparisons.
- Simplified the overall structure of the example for better clarity and focus on core functionalities.

* Enhance monogenic curvelet transform documentation and implementation

- Updated documentation to clarify the 2D limitation of the monogenic curvelet transform, emphasizing that it is mathematically defined only for 2D signals and that only the first two Riesz components (R_1 and R_2) are computed for higher dimensions.
- Introduced a new method in the UDCT class to compute the monogenic signal directly, providing a quaternion-like representation and simplifying the transformation process.
- Added detailed parameter descriptions and examples to improve usability and understanding of the monogenic functionality.

* Update monogenic curvelet example for improved clarity and accuracy

- Enhanced documentation with clearer mathematical definitions and references.
- Adjusted window overlap parameter for the UDCT transform to optimize reconstruction accuracy.
- Improved component comparison visualizations and statistics for scalar and Riesz components.
- Streamlined the layout of plots and print statements for better readability and output clarity.
- Removed redundant comments and sections to focus on core functionalities and results.

* Add optimal window overlap computation to UDCT class

- Introduced a new static method to compute the optimal window overlap based on Nguyen & Chauris (2010) constraint, enhancing reconstruction quality.
- Updated the UDCT methods to utilize the computed window overlap when not provided, replacing hardcoded defaults for improved flexibility.
- Enhanced documentation with detailed explanations, parameters, and examples for the new method, ensuring better usability and understanding.

* Refactor monogenic curvelet transform to support variable-length Riesz components

- Updated the implementation to handle an arbitrary number of Riesz components, enhancing flexibility for ND signals.
- Modified type annotations and return structures to reflect the new list format for coefficients, replacing fixed tuples.
- Improved documentation to clarify the extension of the monogenic transform to arbitrary dimensions and the handling of Riesz components.
- Adjusted tests to verify the correct number of components returned and ensure consistency with the new implementation.

* Refactor UDCT transform initialization and add tests for optimal window overlap

- Removed hardcoded window_overlap parameter in UDCT initialization for improved flexibility.
- Added a comprehensive test suite for the UDCT._compute_optimal_window_overlap method, validating its correctness and constraints.
- Ensured that auto-computed window_overlap is utilized when not explicitly provided, enhancing usability.

* Refactor UDCT and related components for improved clarity and functionality

- Renamed parameters and adjusted type annotations for consistency, particularly in the UDCT class, enhancing code readability.
- Updated documentation to clarify the extension of the monogenic curvelet transform to N-D signals and improved descriptions of parameters.
- Streamlined the implementation of forward and backward monogenic transforms, ensuring better handling of Riesz components.
- Enhanced test cases for monogenic transformations, ensuring correct component counts and improved error messages for assertions.
- Removed redundant comments and improved code structure for better maintainability.

* Refactor backward transform functions to improve code clarity and maintainability

- Added pylint disable comments for nested blocks and duplicate code to enhance readability.
- Streamlined the handling of wavelet and curvelet modes in the backward transform implementations.
- Removed unnecessary imports and adjusted type annotations for consistency across the UDCT module.

* Update type annotations in _typing.py for compatibility with Python 3.10

- Adjusted type annotations for UDCTCoefficients, UDCTWindows, and MUDCTCoefficients to utilize modern list and tuple syntax.
- Updated import statements to reflect changes in typing for improved clarity and maintainability.
- Ensured compatibility with Python versions 3.10 and above by modifying conditional imports.

* Add MUDCTCoefficients to documentation and update imports

- Included MUDCTCoefficients in the Sphinx nitpick_ignore list for improved documentation clarity.
- Updated the __all__ list in the numpy module to include MUDCTCoefficients, ensuring it is publicly accessible.
- Adjusted import statements in the numpy module to reflect the addition of MUDCTCoefficients, enhancing type consistency across the codebase.

* Update type annotations in _typing.py for Python 3.10 compatibility

- Adjusted the import statement for Union to enhance clarity and maintain consistency with modern typing practices.
- Updated the type annotations for UDCTCoefficients and UDCTWindows to utilize the latest list syntax, ensuring compatibility with Python 3.10 and above.

* Add comprehensive tests for MeyerWavelet class

- Introduced a new test suite for the MeyerWavelet class, covering 1D and multi-dimensional forward and inverse transforms.
- Implemented tests for various input types, including real and complex signals, ensuring shape and dtype preservation.
- Added round-trip accuracy tests for different data types and dimensions.
- Included tests for internal methods and edge cases in the transform process, enhancing overall test coverage.

* Refactor test cases in test_udct_from_sparse.py for improved readability

- Reformatted the parameterization of the test_from_sparse_different_dtypes method for better clarity.
- Enhanced the readability of nested loops in the test cases by adjusting line breaks and indentation.
- Ensured that the tests maintain their functionality while improving overall code structure.

* Remove dependency on pre-commit in CI workflow for Python checks

* Refactor test cases in test_meyerwavelet.py and test_udct_from_sparse.py for improved clarity

- Simplified the assertion for dtype preservation in the MeyerWavelet tests by consolidating the absolute tolerance logic.
- Removed unnecessary parameters from test methods in both test suites, enhancing readability and reducing complexity.
- Ensured that the functionality of the tests remains intact while improving overall code structure.

Author: cako

.github/workflows/ci.yml

@@ -36,7 +36,6 @@ jobs:
   checks:
     name: Check Python ${{ matrix.python-version }} on ${{ matrix.runs-on }}
     runs-on: ${{ matrix.runs-on }}
-    needs: [pre-commit]
     strategy:
       fail-fast: false
       matrix:

docs/conf.py

@@ -54,12 +54,14 @@
     ("py:class", "U"),
     ("py:class", "optional"),
     ("py:class", '"curvelet"'),
-    ("py:class", '"meyer"'),
     ("py:class", '"wavelet"'),
     ("py:class", '{"curvelet"'),
     ("py:class", '"wavelet"}'),
     ("py:class", "ParamUDCT"),
     ("py:class", "UDCTWindows"),
+    ("py:class", "MUDCTCoefficients"),
+    ("py:obj", "MUDCTCoefficients"),
+    ("py:obj", "curvelets.numpy.MUDCTCoefficients"),
 ]
 
 always_document_param_types = True

docs/index.rst

@@ -63,14 +63,15 @@ Features
 +--------------------------+----------------------------+
 | Wavelet at highest scale |   ✅  [#wavelet-ftn]_      |
 +--------------------------+----------------------------+
-| Monogenic coefficients   |          ❌                |
+| Monogenic coefficients   |       ✅ [#monogenic-ftn]_ |
 +--------------------------+----------------------------+
 | PyTorch bindings         |          ❌                |
 +--------------------------+----------------------------+
 
 .. [#real-ftn] Supports real inputs with reduced storage requirements which exploit the symmetry of the real-valued Fourier transform.
 .. [#asymmetric-ftn] The directional resolution is asymmetric in the sense that the number of wedges per direction is different for each direction. See :ref:`sphx_glr_auto_examples_plot_02_direction_resolution.py` for an example.
-.. [#wavelet-ftn] Meyer and isotropic wavelets are supported. See :ref:`sphx_glr_auto_examples_plot_05_curvelet_vs_wavelet.py` for an example.
+.. [#wavelet-ftn] Isotropic wavelets are supported. See :ref:`sphx_glr_auto_examples_plot_05_curvelet_vs_wavelet.py` for an example.
+.. [#monogenic-ftn] The monogenic curvelet transform was originally defined for 2D signals by Storath 2010, but this implementation extends it to arbitrary N-D signals by using all Riesz transform components (one per dimension).
 
 
 Credits

examples/plot_05_curvelet_vs_wavelet.py

@@ -1,10 +1,10 @@
 """
 Wavelet Transforms at the Highest Scale
 =============================================
-This example compares curvelet, Meyer, and wavelet transforms at the highest scale
+This example compares curvelet and wavelet transforms at the highest scale
 (using a simple nscale=2) on a zone plate test image. It demonstrates the difference
-between directional curvelet windows, Meyer wavelet windows, and a
-ring-shaped highpass window, showing both frequency-domain windows and spatial coefficients.
+between directional curvelet windows and a ring-shaped highpass window,
+showing both frequency-domain windows and spatial coefficients.
 """
 
 from __future__ import annotations
@@ -26,20 +26,14 @@
 shape = (256, 256)
 zone_plate = make_zone_plate(shape)
 
-# Create three UDCT transforms with num_scales=2
-# One in curvelet mode, one in Meyer mode, one in wavelet mode
+# Create two UDCT transforms with num_scales=2
+# One in curvelet mode, one in wavelet mode
 C_curvelet = UDCT(
     shape=shape,
     num_scales=2,
     wedges_per_direction=3,
     high_frequency_mode="curvelet",
 )
-C_meyer = UDCT(
-    shape=shape,
-    num_scales=2,
-    wedges_per_direction=3,
-    high_frequency_mode="meyer",
-)
 C_wavelet = UDCT(
     shape=shape,
     num_scales=2,
@@ -52,15 +46,12 @@
 # ##################
 
 coeffs_curvelet = C_curvelet.forward(zone_plate)
-coeffs_meyer = C_meyer.forward(zone_plate)
 coeffs_wavelet = C_wavelet.forward(zone_plate)
 
 print(f"Input shape: {zone_plate.shape}")  # noqa: T201
 print(f"Curvelet scales: {len(coeffs_curvelet)}")  # noqa: T201
-print(f"Meyer scales: {len(coeffs_meyer)}")  # noqa: T201
 print(f"Wavelet scales: {len(coeffs_wavelet)}")  # noqa: T201
 print(f"Curvelet scale 1 directions: {len(coeffs_curvelet[1])}")  # noqa: T201
-print(f"Meyer scale 1 bands: {len(coeffs_meyer[1])}")  # noqa: T201
 print(f"Wavelet scale 1 directions: {len(coeffs_wavelet[1])}")  # noqa: T201
 
 # %%
@@ -82,10 +73,10 @@
 # %%
 # Frequency Domain Windows at Highest Scale (Scale 1)
 # ###################################################
-# Visualize the frequency domain windows for curvelet vs Meyer vs wavelet at scale 1.
-# curvelet windows are directional (wedges), Meyer windows are separable
-# (product of 1D filters), and wavelet window is a single ring-shaped window
-# that encompasses the entire high-frequency ring (complement of lowpass filter).
+# Visualize the frequency domain windows for curvelet vs wavelet at scale 1.
+# Curvelet windows are directional (wedges), and wavelet window is a single
+# ring-shaped window that encompasses the entire high-frequency ring
+# (complement of lowpass filter).
 
 # Frequency coordinates for axis labels
 nx, ny = shape
@@ -103,34 +94,14 @@
         curvelet_windows_scale1.append(window_shifted)
         curvelet_window_info.append((idir, iwedge))
 
-# Extract Meyer windows for scale 1 (highpass bands)
-# Access the Meyer wavelet filters from the UDCT object
-if C_meyer._meyer_wavelet is None:
-    msg = "Meyer wavelet not initialized"
-    raise RuntimeError(msg)
-lowpass_1d, highpass_1d = C_meyer._meyer_wavelet._filters[shape[0]]
-
-# Construct 2D frequency domain windows using outer products
-# For 2D Meyer wavelet, the windows are separable (product of 1D filters)
-meyer_window_0_2d = np.outer(lowpass_1d, highpass_1d)  # Low-High
-meyer_window_1_2d = np.outer(highpass_1d, lowpass_1d)  # High-Low
-meyer_window_2_2d = np.outer(highpass_1d, highpass_1d)  # High-High
-
-# Apply fftshift to center the frequency domain for visualization
-meyer_windows_scale1 = [
-    np.real(fftshift(meyer_window_0_2d)).astype(np.float64),
-    np.real(fftshift(meyer_window_1_2d)).astype(np.float64),
-    np.real(fftshift(meyer_window_2_2d)).astype(np.float64),
-]
-
 # Extract wavelet window for scale 1 (single ring-shaped window, complement of lowpass)
 wavelet_window_sparse = C_wavelet.windows[1][0][0]
 wavelet_window_dense = C_wavelet.from_sparse(wavelet_window_sparse)
 wavelet_window_shifted = fftshift(wavelet_window_dense)
 wavelet_windows_scale1 = [wavelet_window_shifted]
 
 # Find common vmax for all windows
-all_windows = curvelet_windows_scale1 + meyer_windows_scale1 + wavelet_windows_scale1
+all_windows = curvelet_windows_scale1 + wavelet_windows_scale1
 vmax_windows = max(np.abs(w).max() for w in all_windows)
 window_opts = {
     "aspect": "equal",
@@ -140,13 +111,12 @@
     "extent": [kx[0], kx[-1], ky[-1], ky[0]],
 }
 
-# Plot curvelet, meyer, and wavelet windows
+# Plot curvelet and wavelet windows
 # sphinx_gallery_thumbnail_number = 2
 n_curvelet_windows = len(curvelet_windows_scale1)
-n_meyer_windows = len(meyer_windows_scale1)
 n_wavelet_windows = len(wavelet_windows_scale1)
-n_cols = max(n_curvelet_windows, n_meyer_windows, n_wavelet_windows)
-fig, axs = plt.subplots(3, n_cols, figsize=(4 * n_cols, 12))
+n_cols = max(n_curvelet_windows, n_wavelet_windows)
+fig, axs = plt.subplots(2, n_cols, figsize=(4 * n_cols, 8))
 fig.suptitle("Frequency Domain Windows at Highest Scale (Scale 1)", fontsize=14)
 
 # Top row: curvelet windows
@@ -174,38 +144,10 @@
 for i in range(n_curvelet_windows, n_cols):
     axs[0, i].axis("off")
 
-# Middle row: Meyer windows
-meyer_names = [
-    "Meyer Highpass 0\n(Low-High)",
-    "Meyer Highpass 1\n(High-Low)",
-    "Meyer Highpass 2\n(High-High)",
-]
-for i, (window, name) in enumerate(zip(meyer_windows_scale1, meyer_names)):
-    ax = axs[1, i]
-    window_real = np.real(window).astype(np.float64)
-    im = ax.imshow(window_real.T, **window_opts)
-    _, cb = create_colorbar(im=im, ax=ax)
-    fmt = ticker.FuncFormatter(lambda x, _: f"{x:.2f}")
-    cb.ax.yaxis.set_major_formatter(fmt)
-    ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.1))
-    ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.1))
-    ax.set(
-        xlim=[kx[0], -kx[0]],
-        ylim=[-ky[0], ky[0]],
-        xlabel="Normalized $k_x$",
-        ylabel="Normalized $k_y$",
-        title=name,
-    )
-    despine(ax)
-
-# Hide unused subplots in middle row
-for i in range(n_meyer_windows, n_cols):
-    axs[1, i].axis("off")
-
 # Bottom row: wavelet window (single ring-shaped window)
 wavelet_names = ["Wavelet\n(Ring)"]
 for i, (window, name) in enumerate(zip(wavelet_windows_scale1, wavelet_names)):
-    ax = axs[2, i]
+    ax = axs[1, i]
     window_real = np.real(window).astype(np.float64)
     im = ax.imshow(window_real.T, **window_opts)
     _, cb = create_colorbar(im=im, ax=ax)
@@ -224,14 +166,14 @@
 
 # Hide unused subplots in bottom row
 for i in range(n_wavelet_windows, n_cols):
-    axs[2, i].axis("off")
+    axs[1, i].axis("off")
 
 fig.tight_layout()
 
 # %%
 # Coefficients at Highest Scale (Scale 1)
 # #######################################
-# Visualize the spatial coefficients for curvelet vs Meyer vs wavelet at scale 1.
+# Visualize the spatial coefficients for curvelet vs wavelet at scale 1.
 # These show how the transforms capture different features of the input.
 
 # Extract curvelet coefficients for scale 1
@@ -243,21 +185,11 @@
         curvelet_coeffs_scale1.append(coeff)
         curvelet_coeff_info.append((idir, iwedge))
 
-# Extract meyer coefficients for scale 1 (highpass bands)
-# Structure: coeffs_meyer[1][direction][wedge]
-# Extract 3 coefficients matching the 3 windows: Low-High, High-Low, High-High
-# Use direction 0, wedges 0, 1, 2 correspond to the 3 highpass bands
-meyer_coeffs_scale1 = [
-    coeffs_meyer[1][0][0],  # Low-High
-    coeffs_meyer[1][0][1],  # High-Low
-    coeffs_meyer[1][0][2],  # High-High
-]
-
 # Extract wavelet coefficient for scale 1 (single coefficient)
 wavelet_coeffs_scale1 = [coeffs_wavelet[1][0][0]]
 
 # Find common vmax for amplitude visualization
-all_coeffs = curvelet_coeffs_scale1 + meyer_coeffs_scale1 + wavelet_coeffs_scale1
+all_coeffs = curvelet_coeffs_scale1 + wavelet_coeffs_scale1
 vmax_coeffs = max(np.abs(c).max() for c in all_coeffs)
 coeff_opts = {
     "aspect": "equal",
@@ -268,10 +200,9 @@
 
 # Plot coefficients
 n_curvelet_coeffs = len(curvelet_coeffs_scale1)
-n_meyer_coeffs = len(meyer_coeffs_scale1)
 n_wavelet_coeffs = len(wavelet_coeffs_scale1)
-n_cols = max(n_curvelet_coeffs, n_meyer_coeffs, n_wavelet_coeffs)
-fig, axs = plt.subplots(3, n_cols, figsize=(4 * n_cols, 12))
+n_cols = max(n_curvelet_coeffs, n_wavelet_coeffs)
+fig, axs = plt.subplots(2, n_cols, figsize=(4 * n_cols, 8))
 fig.suptitle("Coefficients at Highest Scale (Scale 1)", fontsize=14)
 
 # Top row: curvelet coefficients
@@ -294,31 +225,10 @@
 for i in range(n_curvelet_coeffs, n_cols):
     axs[0, i].axis("off")
 
-# Middle row: Meyer coefficients
-meyer_coeff_names = [
-    "Meyer Highpass 0\n(Low-High)",
-    "Meyer Highpass 1\n(High-Low)",
-    "Meyer Highpass 2\n(High-High)",
-]
-for i, (coeff, name) in enumerate(zip(meyer_coeffs_scale1, meyer_coeff_names)):
-    ax = axs[1, i]
-    # Meyer coefficients are real-valued
-    coeff_real = np.real(coeff).astype(np.float64)
-    im = ax.imshow(coeff_real.T, **coeff_opts)
-    _, cb = create_colorbar(im=im, ax=ax)
-    fmt = ticker.FuncFormatter(lambda x, _: f"{x:.2e}")
-    cb.ax.yaxis.set_major_formatter(fmt)
-    ax.set(title=name)
-    despine(ax)
-
-# Hide unused subplots in middle row
-for i in range(n_meyer_coeffs, n_cols):
-    axs[1, i].axis("off")
-
 # Bottom row: wavelet coefficient (single coefficient)
 wavelet_coeff_names = ["Wavelet\n(Ring)"]
 for i, (coeff, name) in enumerate(zip(wavelet_coeffs_scale1, wavelet_coeff_names)):
-    ax = axs[2, i]
+    ax = axs[1, i]
     # wavelet coefficient is real-valued
     coeff_real = np.real(coeff).astype(np.float64)
     im = ax.imshow(coeff_real.T, **coeff_opts)
@@ -330,7 +240,7 @@
 
 # Hide unused subplots in bottom row
 for i in range(n_wavelet_coeffs, n_cols):
-    axs[2, i].axis("off")
+    axs[1, i].axis("off")
 
 fig.tight_layout()
 
@@ -346,12 +256,6 @@
         f"  Dir {idir} Wedge {iwedge}: Energy={energy:.2e}, Max={max_val:.2e}"
     )
 
-print("\nMeyer Scale 1 Statistics:")  # noqa: T201
-for i, coeff in enumerate(meyer_coeffs_scale1):
-    energy = np.sum(np.abs(coeff) ** 2)
-    max_val = np.abs(coeff).max()
-    print(f"  Highpass {i}: Energy={energy:.2e}, Max={max_val:.2e}")  # noqa: T201
-
 print("\nWavelet Scale 1 Statistics:")  # noqa: T201
 for i, coeff in enumerate(wavelet_coeffs_scale1):
     energy = np.sum(np.abs(coeff) ** 2)