diff --git a/Interpolation.hpp b/Interpolation.hpp index f8ca311..c3e3fad 100644 --- a/Interpolation.hpp +++ b/Interpolation.hpp @@ -1,26 +1,120 @@ +/** + * @file Interpolation.hpp + * @brief Provides various interpolation methods, including linear, cubic Hermite, cubic Lagrange, and cubic B-spline interpolations. + * + * This file defines template structures and methods for performing different types of interpolation. + * Each interpolation method allows for smooth estimation of values between known data points + * using specific mathematical techniques: + * - Linear interpolation computes values by connecting points with straight lines. + * - Cubic Hermite interpolation provides smooth results with continuous first derivatives. + * - Cubic Lagrange interpolation uses Lagrange polynomials to approximate values between points. + * - Cubic B-spline interpolation uses B-spline basis functions to ensure smooth curves with continuity in function values and derivatives. + * + * The templates are designed to be flexible and accept various data types for interpolation, such as float or double. + */ + #ifndef INTERPOLATION_HPP #define INTERPOLATION_HPP // Enumeration of interpolation types -enum class InterpType { None, Linear, CubicHermite, CubicLagrange, CubicBSpline }; +/** + * @enum InterpType + * @brief Defines different types of interpolation methods. + * + * This enumeration specifies various interpolation techniques that can be + * used to estimate values between known data points. + */ +enum class InterpType { + None, /**< No interpolation applied. */ + Linear, /**< Linear interpolation. Computes the value by connecting data points with a straight line. */ + CubicHermite, /**< Cubic Hermite interpolation. Provides smooth results with continuity in both the function and its first derivative. */ + CubicLagrange, /**< Cubic Lagrange interpolation. Uses the Lagrange polynomial for approximating values. */ + CubicBSpline /**< Cubic B-Spline interpolation. Uses B-spline curves for smooth interpolation. */ +}; // Linear +/** + * @struct linear_interp + * @brief A template structure for performing linear interpolation. + * + * This structure provides a mechanism for calculating the linear interpolation + * of a given value between two points. The template parameter @p T represents + * the data type of the values being interpolated. + * + * @tparam T The type of the data points to be interpolated (e.g., float, double). + */ + template struct linear_interp { + + /** + * @brief Performs linear interpolation between two values. + * + * This operator overload computes the linear interpolation between two values @p y0 and @p y1 + * at a point @p x, where @p x is typically a fractional value between 0 and 1. + * + * @param x The interpolation factor. Typically a value between 0 and 1, representing + * the relative position between @p y0 and @p y1. + * @param y0 The starting value (interpolated at @p x = 0). + * @param y1 The ending value (interpolated at @p x = 1). + * @return The interpolated value between @p y0 and @p y1. + */ + T operator()(const T& x, const T& y0, const T& y1) { return (y0 + x * ((y1 - y0))); } }; // Cubic Hermite +/** + * @struct cubic_hermite_interp + * @brief A template structure for performing cubic Hermite interpolation. + * + * This structure provides a mechanism for calculating cubic Hermite interpolation, + * which is a method of interpolation that provides smooth curves with continuity + * in both the function values and their first derivatives. + * + * The template parameter @p T represents the data type of the values being interpolated. + * + * @tparam T The type of the data points to be interpolated (e.g., float, double). + */ + template struct cubic_hermite_interp { + + /** + * @brief Default constructor for the cubic Hermite interpolation structure. + * + * Initializes constants used in the cubic Hermite interpolation process. + * These constants represent fractional coefficients that are used in the interpolation formula: + * - @p _5div2 is initialized to 2.5. + * - @p _3div2 is initialized to 1.5. + * - @p _1div2 is initialized to 0.5. + */ + cubic_hermite_interp() : _5div2(2.5), _3div2(1.5), _1div2(0.5) {} + /** + * @brief Performs cubic Hermite interpolation between four values. + * + * This operator overload computes the cubic Hermite interpolation between four points + * @p y0, @p y1, @p y2, and @p y3 at a given point @p x. The values @p y1 and @p y2 + * represent the interval over which interpolation occurs, while @p y0 and @p y3 are + * used to compute the tangent at the interval boundaries, ensuring smooth transitions. + * + * @param x The interpolation factor, typically between 0 and 1, representing the relative + * position between @p y1 and @p y2. + * @param y0 The value before the starting point of interpolation, used to calculate the slope at @p y1. + * @param y1 The starting value of the interpolation interval (interpolated at @p x = 0). + * @param y2 The ending value of the interpolation interval (interpolated at @p x = 1). + * @param y3 The value after the ending point of interpolation, used to calculate the slope at @p y2. + * @return The interpolated value between @p y1 and @p y2 using cubic Hermite interpolation. + */ + T operator()(const T& x, const T& y0, const T& y1, const T& y2, const T& y3) { const T c0 = y1; @@ -33,18 +127,81 @@ struct cubic_hermite_interp private: + /** + * @brief A constant representing the value 5/2 (or 2.5), used in cubic Hermite interpolation calculations. + * + * This constant is used internally in the cubic Hermite interpolation formula to scale the results + * based on specific mathematical operations. + */ + const T _5div2; + + /** + * @brief A constant representing the value 3/2 (or 1.5), used in cubic Hermite interpolation calculations. + * + * This constant is utilized in the cubic Hermite interpolation formula to help scale and calculate + * the interpolated value between points. + */ + const T _3div2; + + /** + * @brief A constant representing the value 1/2 (or 0.5), used in cubic Hermite interpolation calculations. + * + * This constant is used as part of the cubic Hermite interpolation formula to help adjust the scaling + * of the interpolated value between points. + */ + const T _1div2; }; // Cubic Lagrange +/** + * @struct cubic_lagrange_interp + * @brief A template structure for performing cubic Lagrange interpolation. + * + * This structure provides a mechanism for calculating cubic Lagrange interpolation, + * which approximates a smooth curve through four data points using Lagrange polynomials. + * The template parameter @p T represents the data type of the values being interpolated. + * + * @tparam T The type of the data points to be interpolated (e.g., float, double). + */ + template struct cubic_lagrange_interp { + + /** + * @brief Default constructor for the cubic Lagrange interpolation structure. + * + * Initializes constants used in the cubic Lagrange interpolation process. + * These constants represent fractional coefficients that are used in the interpolation formula: + * - @p _1div3 is initialized to 1/3. + * - @p _1div6 is initialized to 1/6. + * - @p _1div2 is initialized to 0.5. + * + * The values are computed based on the template parameter type @p T. + */ + cubic_lagrange_interp() : _1div3(T(1)/T(3)), _1div6(T(1)/T(6)), _1div2(0.5) {} + /** + * @brief Performs cubic Lagrange interpolation between four values. + * + * This operator overload computes the cubic Lagrange interpolation between four points + * @p y0, @p y1, @p y2, and @p y3 at a given point @p x. Lagrange interpolation uses + * polynomials to estimate the value at @p x based on these four surrounding points. + * + * @param x The interpolation factor, typically between 0 and 1, representing the relative position + * within the interval between @p y1 and @p y2. + * @param y0 The value before the starting point of interpolation, used to calculate the Lagrange polynomial. + * @param y1 The starting value of the interpolation interval (interpolated at @p x = 0). + * @param y2 The ending value of the interpolation interval (interpolated at @p x = 1). + * @param y3 The value after the ending point of interpolation, used to calculate the Lagrange polynomial. + * @return The interpolated value between @p y1 and @p y2 using cubic Lagrange interpolation. + */ + T operator()(const T& x, const T& y0, const T& y1, const T& y2, const T& y3) { const T c0 = y1; @@ -57,18 +214,83 @@ struct cubic_lagrange_interp private: + /** + * @brief A constant representing the value 1/3, used in cubic Lagrange interpolation calculations. + * + * This constant is utilized as part of the cubic Lagrange interpolation formula to scale + * and compute the interpolated value between points. + */ + const T _1div3; + + /** + * @brief A constant representing the value 1/6, used in cubic Lagrange interpolation calculations. + * + * This constant is part of the cubic Lagrange interpolation formula, helping to scale and + * compute the interpolated value between data points. + */ + const T _1div6; + + /** + * @brief A constant representing the value 1/2 (or 0.5), used in cubic Lagrange interpolation calculations. + * + * This constant is utilized in the cubic Lagrange interpolation formula to help adjust + * and scale the interpolated value between data points. + */ + const T _1div2; }; // Cubic B-spline +/** + * @struct cubic_bspline_interp + * @brief A template structure for performing cubic B-spline interpolation. + * + * This structure provides a mechanism for calculating cubic B-spline interpolation, + * which produces smooth curves through control points using B-spline basis functions. + * The cubic B-spline method ensures continuity in both the function and its first two derivatives. + * The template parameter @p T represents the data type of the values being interpolated. + * + * @tparam T The type of the data points to be interpolated (e.g., float, double). + */ + template struct cubic_bspline_interp { + + /** + * @brief Default constructor for the cubic B-spline interpolation structure. + * + * Initializes constants used in the cubic B-spline interpolation process. + * These constants represent fractional coefficients that are used in the B-spline interpolation formula: + * - @p _2div3 is initialized to 2/3. + * - @p _1div6 is initialized to 1/6. + * - @p _1div2 is initialized to 0.5. + * + * The values are computed based on the template parameter type @p T. + */ + cubic_bspline_interp() : _2div3(T(2)/T(3)), _1div6(T(1)/T(6)), _1div2(0.5) {} + /** + * @brief Performs cubic B-spline interpolation between four values. + * + * This operator overload computes the cubic B-spline interpolation between four control points + * @p y0, @p y1, @p y2, and @p y3 at a given point @p x. B-spline interpolation uses basis functions + * to create a smooth curve through the control points, ensuring continuity in both the function + * and its first two derivatives. + * + * @param x The interpolation factor, typically between 0 and 1, representing the relative position + * within the interval between @p y1 and @p y2. + * @param y0 The value before the starting point of interpolation, used to define the curve. + * @param y1 The starting value of the interpolation interval (interpolated at @p x = 0). + * @param y2 The ending value of the interpolation interval (interpolated at @p x = 1). + * @param y3 The value after the ending point of interpolation, used to define the curve. + * @return The interpolated value between @p y1 and @p y2 using cubic B-spline interpolation. + */ + T operator()(const T& x, const T& y0, const T& y1, const T& y2, const T& y3) { const T y0py2 = y0 + y2; @@ -82,8 +304,31 @@ struct cubic_bspline_interp private: + /** + * @brief A constant representing the value 2/3, used in cubic B-spline interpolation calculations. + * + * This constant is part of the cubic B-spline interpolation formula, used to scale and compute + * the interpolated value between control points. + */ + const T _2div3; + + /** + * @brief A constant representing the value 1/6, used in cubic B-spline interpolation calculations. + * + * This constant is utilized as part of the cubic B-spline interpolation formula to help scale + * and compute the interpolated value between control points. + */ + const T _1div6; + + /** + * @brief A constant representing the value 1/2 (or 0.5), used in cubic B-spline interpolation calculations. + * + * This constant is part of the cubic B-spline interpolation formula and is used to help scale + * the interpolated value between control points. + */ + const T _1div2; }; diff --git a/KernelSmoother.hpp b/KernelSmoother.hpp index 4bc24af..ed31833 100644 --- a/KernelSmoother.hpp +++ b/KernelSmoother.hpp @@ -1,4 +1,24 @@ +/** + * @file KernelSmoother.hpp + * + * @brief This file defines the `kernel_smoother` class and associated methods for applying kernel smoothing + * and spectral processing on input data. + * + * The `kernel_smoother` class provides functionality for smoothing data using various kernel-based methods, + * including FFT-based approaches for efficient convolution. The class supports different edge-handling strategies + * and allows users to specify filter parameters such as width, gain, and symmetry. The file also includes utility + * classes and methods for managing input data, applying filters, and handling FFT operations. + * + * Key Components: + * - `kernel_smoother`: The main class that provides kernel smoothing capabilities. + * - Edge handling strategies through `Ends` and `EdgeMode` enumerations. + * - FFT-based filtering methods for efficient convolution. + * - Utility methods for creating and applying filters to data. + * + * @note This file contains template-based functions and classes to support a wide range of data types and customizable behavior. + */ + #ifndef KERNELSMOOTHER_HPP #define KERNELSMOOTHER_HPP @@ -14,44 +34,253 @@ #include "SpectralProcessor.hpp" #include "TableReader.hpp" +/** + * @brief A class that performs kernel smoothing on spectral data. + * + * This template class `kernel_smoother` is designed to smooth spectral data using a specified kernel. + * It inherits from `spectral_processor` to leverage spectral processing capabilities, + * and operates on data types defined by the template parameters. + * + * @tparam T The type of data being processed (e.g., float, double). + * @tparam Allocator The allocator type for memory management, defaults to `aligned_allocator` for optimized memory alignment. + * @tparam auto_resize_fft A boolean flag indicating whether the FFT size should automatically resize to accommodate different data sizes. + * Defaults to `false`, meaning the FFT size remains fixed unless manually changed. + */ + template class kernel_smoother : private spectral_processor { + + /** + * @typedef processor + * + * An alias for the `spectral_processor` class template specialized with the types `T` and `Allocator`. + * + * This typedef simplifies access to the `spectral_processor` base class, allowing the `kernel_smoother` + * class to refer to it more conveniently. + * + * @tparam T The data type being processed by the spectral processor (e.g., float, double). + * @tparam Allocator The allocator type used for memory management in the spectral processor. + */ + using processor = spectral_processor; + + /** + * @typedef op_sizes + * + * An alias for the `op_sizes` type defined in the `processor` class (i.e., `spectral_processor`). + * + * This typedef provides access to the internal structure or type that handles operational sizes + * within the `spectral_processor` class, which could involve FFT sizes or buffer lengths used + * during spectral processing. + * + * @tparam T The data type being processed by the spectral processor. + * @tparam Allocator The allocator type used in the spectral processor. + */ + using op_sizes = typename processor::op_sizes; + + /** + * @typedef zipped_pointer + * + * An alias for the `zipped_pointer` type defined in the `processor` class (i.e., `spectral_processor`). + * + * This typedef provides access to a specialized pointer type used in the `spectral_processor` class, + * likely designed to point to multiple data arrays (zipped together) for efficient processing in + * spectral operations. + * + * @tparam T The data type being processed by the spectral processor. + * @tparam Allocator The allocator type used in the spectral processor. + */ + using zipped_pointer = typename processor::zipped_pointer; + + /** + * @typedef in_ptr + * + * An alias for the `in_ptr` type defined in the `processor` class (i.e., `spectral_processor`). + * + * This typedef represents a pointer type used to access input data in the `spectral_processor` class, + * allowing for efficient manipulation and processing of the input data during spectral operations. + * + * @tparam T The data type being processed by the spectral processor. + * @tparam Allocator The allocator type used in the spectral processor. + */ + using in_ptr = typename processor::in_ptr; + /** + * @typedef Split + * + * An alias for the `Split` type defined within the `FFTTypes` class. + * + * This typedef simplifies the usage of the `Split` type for the current + * template type `T`, allowing easier access to a specialized FFT split data structure. + * + * @tparam T The data type for which the FFT types and associated split types are defined. + */ + using Split = typename FFTTypes::Split; + /** + * @typedef enable_if_t + * + * A template alias for `std::enable_if`, which is used for SFINAE (Substitution Failure Is Not An Error) + * to conditionally enable or disable template specializations or functions based on a boolean constant. + * + * If the boolean value `B` is true, `enable_if_t` resolves to `int`. Otherwise, this alias is invalid, + * and the function or specialization is excluded from the overload set. + * + * @tparam B A boolean constant used to control whether the alias is enabled or disabled. + */ + template using enable_if_t = typename std::enable_if::type; + /** + * @enum Ends + * + * @brief An enumeration that defines different types of endpoint behaviors for kernel smoothing or spectral processing. + * + * This enum is used to specify how the ends or boundaries of the data should be handled during processing. + * Different options allow for various treatments, such as zero-padding or maintaining non-zero values. + * + * - `Zero`: The endpoints are set to zero. + * - `NonZero`: The endpoints maintain their non-zero values. + * - `SymZero`: The endpoints are symmetrically zero-padded. + * - `SymNonZero`: The endpoints are symmetrically handled with non-zero values. + */ + enum class Ends { Zero, NonZero, SymZero, SymNonZero }; public: + /** + * @enum EdgeMode + * + * @brief An enumeration that defines various strategies for handling edge conditions in kernel smoothing or spectral processing. + * + * This enum specifies how the edges of the data should be treated when applying filters or transformations, particularly when data at the boundaries is required for computations. + * + * - `ZeroPad`: The edges are padded with zeros. + * - `Extend`: The edge values are extended by repeating the boundary value. + * - `Wrap`: The data wraps around, treating the end as if it continues from the beginning. + * - `Fold`: The data is folded over, mirroring at the edges. + * - `Mirror`: The data is reflected at the boundaries. + */ + enum class EdgeMode { ZeroPad, Extend, Wrap, Fold, Mirror }; + /** + * @brief Constructs a `kernel_smoother` object with an optional maximum FFT size. + * + * This constructor initializes the `kernel_smoother` object, and it is only enabled if the `Allocator` + * type is default constructible. The constructor also initializes the base class `spectral_processor` + * with the given maximum FFT size. + * + * @tparam U The allocator type, which defaults to the `Allocator` type defined for the class. + * This template parameter is only enabled if `U` is default constructible. + * + * @param max_fft_size The maximum size for the FFT, defaulting to \(2^{18}\) (262144). This size determines the largest FFT that the kernel smoother can handle. + * + * @note This constructor is conditionally enabled using SFINAE, meaning it is only available if the `Allocator` type + * is default constructible. + */ + template ::value> = 0> kernel_smoother(uintptr_t max_fft_size = 1 << 18) : spectral_processor(max_fft_size) {} + /** + * @brief Constructs a `kernel_smoother` object with a custom allocator and an optional maximum FFT size. + * + * This constructor initializes the `kernel_smoother` object using a provided allocator and optionally a maximum FFT size. + * It is only enabled if the `Allocator` type is copy constructible. The constructor also passes the allocator + * and FFT size to the base class `spectral_processor`. + * + * @tparam U The allocator type, which defaults to the `Allocator` type defined for the class. + * This template parameter is only enabled if `U` is copy constructible. + * + * @param allocator A reference to the allocator used for managing memory within the `kernel_smoother`. + * @param max_fft_size The maximum size for the FFT, defaulting to \(2^{18}\) (262144). This size determines the largest FFT that the kernel smoother can handle. + * + * @note This constructor is conditionally enabled using SFINAE, meaning it is only available if the `Allocator` type + * is copy constructible. + */ + template ::value> = 0> kernel_smoother(const Allocator& allocator, uintptr_t max_fft_size = 1 << 18) : spectral_processor(allocator, max_fft_size) {} + /** + * @brief Constructs a `kernel_smoother` object using a move-constructed allocator and an optional maximum FFT size. + * + * This constructor initializes the `kernel_smoother` object using a move-constructed allocator and an optional + * maximum FFT size. It is only enabled if the `Allocator` type is move constructible. The allocator and FFT size + * are passed to the base class `spectral_processor`. + * + * @tparam U The allocator type, which defaults to the `Allocator` type defined for the class. + * This template parameter is only enabled if `U` is move constructible. + * + * @param allocator An rvalue reference to the allocator, which is move-constructed to manage memory within the `kernel_smoother`. + * @param max_fft_size The maximum size for the FFT, defaulting to \(2^{18}\) (262144). This size determines the largest FFT that the kernel smoother can handle. + * + * @note This constructor is conditionally enabled using SFINAE, meaning it is only available if the `Allocator` type + * is move constructible. + */ + template ::value> = 0> kernel_smoother(const Allocator&& allocator, uintptr_t max_fft_size = 1 << 18) : spectral_processor(allocator, max_fft_size) {} + /** + * @brief Sets the maximum size for the FFT operations in the kernel smoother. + * + * This method updates the maximum FFT size used by the `kernel_smoother` for spectral processing. It forwards + * the size to the base class `spectral_processor` to configure the internal FFT size limit. + * + * @param size The maximum FFT size to set. This value defines the upper limit for the size of FFT that + * the smoother can handle. + */ + void set_max_fft_size(uintptr_t size) { processor::set_max_fft_size(size); } + /** + * @brief Retrieves the maximum size for FFT operations in the kernel smoother. + * + * This method returns the current maximum FFT size used by the `kernel_smoother`, + * as defined in the base class `spectral_processor`. + * + * @return The maximum FFT size that the kernel smoother can handle. + */ + uintptr_t max_fft_size() { return processor::max_fft_size(); } + /** + * @brief Performs kernel smoothing on the input data. + * + * This method applies kernel smoothing to the input data `in` and stores the result in the output array `out`. + * The smoothing is done using a specified kernel and various parameters that control the width, symmetry, and + * edge behavior. + * + * @param out A pointer to the output array where the smoothed data will be stored. + * @param in A pointer to the input data array that will be smoothed. + * @param kernel A pointer to the kernel array used for smoothing the data. + * @param length The length of the input data array. + * @param kernel_length The length of the kernel array. + * @param width_lo The lower bound width scaling factor for the kernel smoothing. + * @param width_hi The upper bound width scaling factor for the kernel smoothing. + * @param symmetric A boolean flag indicating whether the smoothing should be symmetric. + * If `true`, the kernel is applied symmetrically around the data points. + * @param edges Specifies how the edges of the input data should be handled, using the `EdgeMode` enumeration. + * + * @note The behavior at the edges is controlled by the `EdgeMode` parameter, which determines how data is treated + * at the boundaries during the smoothing process. + */ + void smooth(T *out, const T *in, const T *kernel, uintptr_t length, uintptr_t kernel_length, double width_lo, double width_hi, bool symmetric, EdgeMode edges) { if (!length || !kernel_length) @@ -210,16 +439,73 @@ class kernel_smoother : private spectral_processor private: + /** + * @brief A struct that specializes `table_fetcher` for `double` precision values. + * + * The `fetcher` struct inherits from `table_fetcher` and provides functionality for fetching + * or processing table-based data with double precision. It may include additional functionality + * or specializations specific to the `kernel_smoother`. + * + * @see table_fetcher + */ + struct fetcher : table_fetcher { + + /** + * @brief Constructs a `fetcher` object for fetching data from a given input array. + * + * This constructor initializes the `fetcher` by passing the size of the input data and a scaling factor + * to the base class `table_fetcher`. It also stores a pointer to the input data array. + * + * @param in A pointer to the input data array that the `fetcher` will operate on. + * @param size The size (number of elements) of the input data array. + * + * @note The constructor passes a scaling factor of 1.0 to the `table_fetcher` base class. + */ + fetcher(const T *in, intptr_t size) : table_fetcher(size, 1.0), data(in) {} + /** + * @brief Retrieves the value at the specified index from the input data array. + * + * This function call operator allows `fetcher` objects to be used like functions to access + * the data at a given index. It returns the value from the input data array at the specified index. + * + * @param idx The index of the data element to retrieve. + * @return The value of type `T` located at the given index in the input data array. + */ + T operator()(intptr_t idx) { return data[idx]; } + /** + * @brief A pointer to the input data array. + * + * This member stores the address of the input data that the `fetcher` operates on. + * It is a constant pointer to ensure that the data cannot be modified through this member. + * + * @note The data is of type `T` and is used for fetching or processing elements by the `fetcher`. + */ + const T *data; }; + /** + * @brief Copies the edges of the input data to the output array with a specified filter size. + * + * This method is responsible for copying the edges of the input data to the output data, handling the boundaries according to a given filter size. The method template allows for a customizable edge-handling strategy, determined by the template template parameter `U`. + * + * @tparam U A template template parameter representing a class template that defines the edge handling strategy. + * It takes a single type `V` as its template argument. + * @param in A pointer to the input data array. + * @param out A pointer to the output data array where the edge-handled values will be copied. + * @param length The length of the input data array. + * @param filter_size The size of the filter to be applied to the edges. + * + * @note This function handles edge conditions using the specified template template parameter `U` to define how the edges are copied or processed. + */ + template