N-dimensional bilateral_filter template function implementation

//  N-dimensional bilateral_filter template function implementation
//  bilateral_filter template function implementation
template<typename ElementT, class ExecutionPolicy, class Fr, class Gs, std::integral SizeT = std::size_t>
requires((std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) and
         (std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    const Fr fr,
    const Gs gs,
    const BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    const ElementT value_for_constant_padding = ElementT{}
    )
{
    if (input.getDimensionality() != 2)
    {
        throw std::runtime_error("Unsupported dimension!");
    }
    return windowed_filter(
        std::forward<ExecutionPolicy>(execution_policy),
        input,
        window_size,
        [&](auto&& window) {
            return bilateral_filter_detail(window, fr, gs);
        },
        boundaryCondition,
        value_for_constant_padding
    );
}

//  bilateral_filter template function implementation
template<class ExecutionPolicy, typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
requires ((std::same_as<ElementT, RGB>) ||
          (std::same_as<ElementT, RGB_DOUBLE>) ||
          (std::same_as<ElementT, HSV>) ||
          (is_MultiChannel<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return apply_each(
        input,
        [&](auto&& each_plane)
        {
            return bilateral_filter(
                std::forward<ExecutionPolicy>(execution_policy),
                each_plane,
                window_size,
                fr,
                gs,
                boundaryCondition
                );
        });
}

//  bilateral_filter template function implementation (without Execution Policy)
template<typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
constexpr static auto bilateral_filter(
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return bilateral_filter(
        std::execution::seq,
        input,
        window_size,
        fr,
        gs,
        boundaryCondition
    );
}

bilateral_filter_detail template function implementation

//  bilateral_filter_detail template function implementation
template<typename ElementT, class Fr, class Gs, class FloatingType = double>
requires(std::invocable<Fr, ElementT> && std::invocable<Gs, std::size_t>)
constexpr static auto bilateral_filter_detail(
    const Image<ElementT>& input,
    Fr fr,
    Gs gs
)
{
    const ElementT center_pixel = get_center_pixel(std::execution::seq, input);
    auto center_location = get_center_location(std::execution::seq, input);
    FloatingType sum{}, weight_sum{};
    const std::size_t total_elements = input.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());
        auto fr_result = std::invoke(
            fr,
            std::abs(input.at_without_boundary_check(indices) - center_pixel)
        );
        auto location_difference = difference(std::execution::seq, indices, center_location);
        auto gs_result = std::invoke(
            gs,
            std::reduce(std::ranges::cbegin(location_difference), std::ranges::cend(location_difference))
        );
        sum += input.at_without_boundary_check(indices) * fr_result * gs_result;
        weight_sum += fr_result * gs_result;
    }

    if (weight_sum == 0)
    {
        return static_cast<ElementT>(sum);
    }
    return static_cast<ElementT>(sum / weight_sum);
}

windowed_filter template function implementation

//  windowed_filter template function implementation
template<class ElementT, class ExecutionPolicy, class Filter, class SizeT = std::size_t>
requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
constexpr static auto windowed_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    SizeT window_size,
    Filter filter,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    ElementT value_for_constant_padding = ElementT{},
    std::function<void(double)> progress_callback = nullptr
)
{
    const std::size_t dim = input.getDimensionality();

    auto padded_image = input;
    if (dim == 2)                   //  padding algorithm supported
    {
        padded_image = generate_padded_image(
            std::forward<ExecutionPolicy>(execution_policy),
            input,
            window_size,
            window_size,
            boundaryCondition,
            value_for_constant_padding
        );
    }
    // Iterate over all elements in the original image
    const std::size_t total_elements = input.count();
    #ifndef USE_OPENMP
    auto indices_view = std::views::iota(std::size_t{ 0 }, total_elements);
    // Progress tracking variables
    std::atomic<std::size_t> processed_count = 0;
    std::atomic<std::size_t> next_report = 0;
    std::mutex report_mutex;
    const std::size_t report_step = std::max<std::size_t>(1, total_elements / 100);
    auto process_element = [&](auto&& idx) {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);

        // Progress reporting
        if (progress_callback) {
            auto count = ++processed_count;
            if (count >= next_report.load(std::memory_order_relaxed)) {
                std::lock_guard lock(report_mutex);
                if (count >= next_report.load(std::memory_order_relaxed)) {
                    double progress = static_cast<double>(count) / total_elements;
                    progress_callback(progress);
                    next_report.store(count + report_step, std::memory_order_relaxed);
                }
            }
        }
        };
    std::for_each(
        std::forward<ExecutionPolicy>(execution_policy),
        indices_view.begin(),
        indices_view.end(),
        process_element
    );
    #else
    #pragma omp parallel for
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);
    }
    #endif
    // Crop padded image to original size
    std::vector<std::size_t> original_sizes = input.getSize();
    std::vector<std::size_t> centers;
    for (auto& s : padded_image.getSize())
    {
        centers.emplace_back(s / 2);  // Center of padded image
    }
    auto output = subimage(std::forward<ExecutionPolicy>(execution_policy), padded_image, original_sizes, centers);
    return output;
}

subimage template function implementation for N dimensional image

//  subimage template function implementation for N dimensional image
template<class ExecutionPolicy, typename ElementT, std::ranges::input_range Sizes, std::ranges::input_range Centers>
requires (
    (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) &&
    (std::same_as<std::ranges::range_value_t<Sizes>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Sizes>, int>) &&
    (std::same_as<std::ranges::range_value_t<Centers>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Centers>, int>)
)
constexpr static auto subimage(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const Sizes& new_sizes,
    const Centers& centers,
    const ElementT default_element = ElementT{})
{
    const std::size_t dim = input.getDimensionality();
    if (new_sizes.size() != dim || centers.size() != dim)
    {
        throw std::runtime_error("subimage: new_sizes and centers must match input dimensionality");
    }

    // Compute start indices for each dimension
    std::vector<std::size_t> start(dim);
    for (std::size_t i = 0; i < dim; ++i)
    {
        start[i] = centers[i] - static_cast<std::size_t>(std::floor(static_cast<double>(new_sizes[i]) / 2.0));
    }

    Image<ElementT> output(new_sizes);

    std::vector<std::size_t> new_sizes_vec;
    new_sizes_vec.resize(std::ranges::size(new_sizes));

    std::transform(
        std::forward<ExecutionPolicy>(execution_policy),
        std::ranges::cbegin(new_sizes),
        std::ranges::cend(new_sizes),
        std::ranges::begin(new_sizes_vec),
        [](auto&& element) { return static_cast<std::size_t>(element); }
    );

    // Precompute strides using new_sizes_vec
    std::vector<std::size_t> strides(dim, 1);
    for (std::size_t i = 1; i < dim; ++i)
    {
        strides[i] = strides[i - 1] * new_sizes_vec[i - 1];
    }

    // Iterate through all output elements
    const std::size_t total_elements = output.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Use helper function to get indices
        const auto output_indices = linear_index_to_indices(idx, strides, new_sizes_vec);

        // Compute input indices and validate
        bool valid = true;
        std::vector<std::size_t> input_indices(dim);
        for (std::size_t i = 0; i < dim; ++i)
        {
            input_indices[i] = start[i] + output_indices[i];
            if (input_indices[i] >= input.getSize(i))
            {
                valid = false;
                break;
            }
        }

        // Assign value
        output.set(idx) = valid ? 
            input.at_without_boundary_check(input_indices) : 
            default_element;
    }

    return output;
}

linear_index_to_indices template function implementation

//  linear_index_to_indices template function implementation
//  Helper function to convert linear index to N-dimensional indices
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    const std::size_t linearIndex,
    const Sizes& strides,
    const Sizes& sizes)
{
    const std::size_t dim = std::ranges::size(sizes);
    std::vector<std::size_t> indices(dim);
    for (std::size_t i = 0; i < dim; ++i) {
        indices[i] = (linearIndex / strides[i]) % sizes[i];
    }
    return indices;
}

//  linear_index_to_indices template function implementation
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    std::size_t linear_idx,
    const Sizes& sizes
)
{
    std::vector<std::size_t> indices;
    std::size_t stride = 1;
    for (auto& s : sizes) {
        indices.emplace_back((linear_idx / stride) % s);
        stride *= s;
    }
    return indices;
}

N-dimensional bilateral_filter template function implementation

//  N-dimensional bilateral_filter template function implementation
//  bilateral_filter template function implementation
template<typename ElementT, class ExecutionPolicy, class Fr, class Gs, std::integral SizeT = std::size_t>
requires((std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) and
         (std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    const Fr fr,
    const Gs gs,
    const BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    const ElementT value_for_constant_padding = ElementT{}
    )
{
    return windowed_filter(
        std::forward<ExecutionPolicy>(execution_policy),
        input,
        window_size,
        [&](auto&& window) {
            return bilateral_filter_detail(window, fr, gs);
        },
        boundaryCondition,
        value_for_constant_padding
    );
}

//  bilateral_filter template function implementation
template<class ExecutionPolicy, typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
requires ((std::same_as<ElementT, RGB>) ||
          (std::same_as<ElementT, RGB_DOUBLE>) ||
          (std::same_as<ElementT, HSV>) ||
          (is_MultiChannel<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return apply_each(
        input,
        [&](auto&& each_plane)
        {
            return bilateral_filter(
                std::forward<ExecutionPolicy>(execution_policy),
                each_plane,
                window_size,
                fr,
                gs,
                boundaryCondition
                );
        });
}

//  bilateral_filter template function implementation (without Execution Policy)
template<typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
constexpr static auto bilateral_filter(
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return bilateral_filter(
        std::execution::seq,
        input,
        window_size,
        fr,
        gs,
        boundaryCondition
    );
}

bilateral_filter_detail template function implementation

//  bilateral_filter_detail template function implementation
template<typename ElementT, class Fr, class Gs, class FloatingType = double>
requires(std::invocable<Fr, ElementT> && std::invocable<Gs, std::size_t>)
constexpr static auto bilateral_filter_detail(
    const Image<ElementT>& input,
    Fr fr,
    Gs gs
)
{
    const ElementT center_pixel = get_center_pixel(std::execution::seq, input);
    auto center_location = get_center_location(std::execution::seq, input);
    FloatingType sum{}, weight_sum{};
    const std::size_t total_elements = input.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());
        auto fr_result = std::invoke(
            fr,
            std::abs(input.at_without_boundary_check(indices) - center_pixel)
        );
        auto location_difference = difference(std::execution::seq, indices, center_location);
        auto gs_result = std::invoke(
            gs,
            std::reduce(std::ranges::cbegin(location_difference), std::ranges::cend(location_difference))
        );
        sum += input.at_without_boundary_check(indices) * fr_result * gs_result;
        weight_sum += fr_result * gs_result;
    }

    if (weight_sum == 0)
    {
        return static_cast<ElementT>(sum);
    }
    return static_cast<ElementT>(sum / weight_sum);
}

windowed_filter template function implementation

//  windowed_filter template function implementation
template<class ElementT, class ExecutionPolicy, class Filter, class SizeT = std::size_t>
requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
constexpr static auto windowed_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    SizeT window_size,
    Filter filter,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    ElementT value_for_constant_padding = ElementT{},
    std::function<void(double)> progress_callback = nullptr
)
{
    const std::size_t dim = input.getDimensionality();

    auto padded_image = input;
    if (dim == 2)                   //  padding algorithm supported
    {
        padded_image = generate_padded_image(
            std::forward<ExecutionPolicy>(execution_policy),
            input,
            window_size,
            window_size,
            boundaryCondition,
            value_for_constant_padding
        );
    }
    // Iterate over all elements in the original image
    const std::size_t total_elements = input.count();
    #ifndef USE_OPENMP
    auto indices_view = std::views::iota(std::size_t{ 0 }, total_elements);
    // Progress tracking variables
    std::atomic<std::size_t> processed_count = 0;
    std::atomic<std::size_t> next_report = 0;
    std::mutex report_mutex;
    const std::size_t report_step = std::max<std::size_t>(1, total_elements / 100);
    auto process_element = [&](auto&& idx) {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);

        // Progress reporting
        if (progress_callback) {
            auto count = ++processed_count;
            if (count >= next_report.load(std::memory_order_relaxed)) {
                std::lock_guard lock(report_mutex);
                if (count >= next_report.load(std::memory_order_relaxed)) {
                    double progress = static_cast<double>(count) / total_elements;
                    progress_callback(progress);
                    next_report.store(count + report_step, std::memory_order_relaxed);
                }
            }
        }
        };
    std::for_each(
        std::forward<ExecutionPolicy>(execution_policy),
        indices_view.begin(),
        indices_view.end(),
        process_element
    );
    #else
    #pragma omp parallel for
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);
    }
    #endif
    // Crop padded image to original size
    std::vector<std::size_t> original_sizes = input.getSize();
    std::vector<std::size_t> centers;
    for (auto& s : padded_image.getSize())
    {
        centers.emplace_back(s / 2);  // Center of padded image
    }
    auto output = subimage(std::forward<ExecutionPolicy>(execution_policy), padded_image, original_sizes, centers);
    return output;
}

subimage template function implementation for N dimensional image

//  subimage template function implementation for N dimensional image
template<class ExecutionPolicy, typename ElementT, std::ranges::input_range Sizes, std::ranges::input_range Centers>
requires (
    (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) &&
    (std::same_as<std::ranges::range_value_t<Sizes>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Sizes>, int>) &&
    (std::same_as<std::ranges::range_value_t<Centers>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Centers>, int>)
)
constexpr static auto subimage(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const Sizes& new_sizes,
    const Centers& centers,
    const ElementT default_element = ElementT{})
{
    const std::size_t dim = input.getDimensionality();
    if (new_sizes.size() != dim || centers.size() != dim)
    {
        throw std::runtime_error("subimage: new_sizes and centers must match input dimensionality");
    }

    // Compute start indices for each dimension
    std::vector<std::size_t> start(dim);
    for (std::size_t i = 0; i < dim; ++i)
    {
        start[i] = centers[i] - static_cast<std::size_t>(std::floor(static_cast<double>(new_sizes[i]) / 2.0));
    }

    Image<ElementT> output(new_sizes);

    std::vector<std::size_t> new_sizes_vec;
    new_sizes_vec.resize(std::ranges::size(new_sizes));

    std::transform(
        std::forward<ExecutionPolicy>(execution_policy),
        std::ranges::cbegin(new_sizes),
        std::ranges::cend(new_sizes),
        std::ranges::begin(new_sizes_vec),
        [](auto&& element) { return static_cast<std::size_t>(element); }
    );

    // Precompute strides using new_sizes_vec
    std::vector<std::size_t> strides(dim, 1);
    for (std::size_t i = 1; i < dim; ++i)
    {
        strides[i] = strides[i - 1] * new_sizes_vec[i - 1];
    }

    // Iterate through all output elements
    const std::size_t total_elements = output.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Use helper function to get indices
        const auto output_indices = linear_index_to_indices(idx, strides, new_sizes_vec);

        // Compute input indices and validate
        bool valid = true;
        std::vector<std::size_t> input_indices(dim);
        for (std::size_t i = 0; i < dim; ++i)
        {
            input_indices[i] = start[i] + output_indices[i];
            if (input_indices[i] >= input.getSize(i))
            {
                valid = false;
                break;
            }
        }

        // Assign value
        output.set(idx) = valid ? 
            input.at_without_boundary_check(input_indices) : 
            default_element;
    }

    return output;
}

linear_index_to_indices template function implementation

//  linear_index_to_indices template function implementation
//  Helper function to convert linear index to N-dimensional indices
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    const std::size_t linearIndex,
    const Sizes& strides,
    const Sizes& sizes)
{
    const std::size_t dim = std::ranges::size(sizes);
    std::vector<std::size_t> indices(dim);
    for (std::size_t i = 0; i < dim; ++i) {
        indices[i] = (linearIndex / strides[i]) % sizes[i];
    }
    return indices;
}

//  linear_index_to_indices template function implementation
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    std::size_t linear_idx,
    const Sizes& sizes
)
{
    std::vector<std::size_t> indices;
    std::size_t stride = 1;
    for (auto& s : sizes) {
        indices.emplace_back((linear_idx / stride) % s);
        stride *= s;
    }
    return indices;
}

N-dimensional bilateral_filter template function implementation

//  N-dimensional bilateral_filter template function implementation
//  bilateral_filter template function implementation
template<typename ElementT, class ExecutionPolicy, class Fr, class Gs, std::integral SizeT = std::size_t>
requires((std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) and
         (std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    const Fr fr,
    const Gs gs,
    const BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    const ElementT value_for_constant_padding = ElementT{}
    )
{
    if (input.getDimensionality() != 2)
    {
        throw std::runtime_error("Unsupported dimension!");
    }
    return windowed_filter(
        std::forward<ExecutionPolicy>(execution_policy),
        input,
        window_size,
        [&](auto&& window) {
            return bilateral_filter_detail(window, fr, gs);
        },
        boundaryCondition,
        value_for_constant_padding
    );
}

//  bilateral_filter template function implementation
template<class ExecutionPolicy, typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
requires ((std::same_as<ElementT, RGB>) ||
          (std::same_as<ElementT, RGB_DOUBLE>) ||
          (std::same_as<ElementT, HSV>) ||
          (is_MultiChannel<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return apply_each(
        input,
        [&](auto&& each_plane)
        {
            return bilateral_filter(
                std::forward<ExecutionPolicy>(execution_policy),
                each_plane,
                window_size,
                fr,
                gs,
                boundaryCondition
                );
        });
}

//  bilateral_filter template function implementation (without Execution Policy)
template<typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
constexpr static auto bilateral_filter(
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return bilateral_filter(
        std::execution::seq,
        input,
        window_size,
        fr,
        gs,
        boundaryCondition
    );
}

windowed_filter template function implementation

//  windowed_filter template function implementation
template<class ElementT, class ExecutionPolicy, class Filter, class SizeT = std::size_t>
requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
constexpr static auto windowed_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    SizeT window_size,
    Filter filter,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    ElementT value_for_constant_padding = ElementT{},
    std::function<void(double)> progress_callback = nullptr
)
{
    const std::size_t dim = input.getDimensionality();

    auto padded_image = input;
    if (dim == 2)                   //  padding algorithm supported
    {
        padded_image = generate_padded_image(
            std::forward<ExecutionPolicy>(execution_policy),
            input,
            window_size,
            window_size,
            boundaryCondition,
            value_for_constant_padding
        );
    }
    // Iterate over all elements in the original image
    const std::size_t total_elements = input.count();
    #ifndef USE_OPENMP
    auto indices_view = std::views::iota(std::size_t{ 0 }, total_elements);
    // Progress tracking variables
    std::atomic<std::size_t> processed_count = 0;
    std::atomic<std::size_t> next_report = 0;
    std::mutex report_mutex;
    const std::size_t report_step = std::max<std::size_t>(1, total_elements / 100);
    auto process_element = [&](auto&& idx) {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);

        // Progress reporting
        if (progress_callback) {
            auto count = ++processed_count;
            if (count >= next_report.load(std::memory_order_relaxed)) {
                std::lock_guard lock(report_mutex);
                if (count >= next_report.load(std::memory_order_relaxed)) {
                    double progress = static_cast<double>(count) / total_elements;
                    progress_callback(progress);
                    next_report.store(count + report_step, std::memory_order_relaxed);
                }
            }
        }
        };
    std::for_each(
        std::forward<ExecutionPolicy>(execution_policy),
        indices_view.begin(),
        indices_view.end(),
        process_element
    );
    #else
    #pragma omp parallel for
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);
    }
    #endif
    // Crop padded image to original size
    std::vector<std::size_t> original_sizes = input.getSize();
    std::vector<std::size_t> centers;
    for (auto& s : padded_image.getSize())
    {
        centers.emplace_back(s / 2);  // Center of padded image
    }
    auto output = subimage(std::forward<ExecutionPolicy>(execution_policy), padded_image, original_sizes, centers);
    return output;
}

subimage template function implementation for N dimensional image

//  subimage template function implementation for N dimensional image
template<class ExecutionPolicy, typename ElementT, std::ranges::input_range Sizes, std::ranges::input_range Centers>
requires (
    (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) &&
    (std::same_as<std::ranges::range_value_t<Sizes>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Sizes>, int>) &&
    (std::same_as<std::ranges::range_value_t<Centers>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Centers>, int>)
)
constexpr static auto subimage(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const Sizes& new_sizes,
    const Centers& centers,
    const ElementT default_element = ElementT{})
{
    const std::size_t dim = input.getDimensionality();
    if (new_sizes.size() != dim || centers.size() != dim)
    {
        throw std::runtime_error("subimage: new_sizes and centers must match input dimensionality");
    }

    // Compute start indices for each dimension
    std::vector<std::size_t> start(dim);
    for (std::size_t i = 0; i < dim; ++i)
    {
        start[i] = centers[i] - static_cast<std::size_t>(std::floor(static_cast<double>(new_sizes[i]) / 2.0));
    }

    Image<ElementT> output(new_sizes);

    std::vector<std::size_t> new_sizes_vec;
    new_sizes_vec.resize(std::ranges::size(new_sizes));

    std::transform(
        std::forward<ExecutionPolicy>(execution_policy),
        std::ranges::cbegin(new_sizes),
        std::ranges::cend(new_sizes),
        std::ranges::begin(new_sizes_vec),
        [](auto&& element) { return static_cast<std::size_t>(element); }
    );

    // Precompute strides using new_sizes_vec
    std::vector<std::size_t> strides(dim, 1);
    for (std::size_t i = 1; i < dim; ++i)
    {
        strides[i] = strides[i - 1] * new_sizes_vec[i - 1];
    }

    // Iterate through all output elements
    const std::size_t total_elements = output.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Use helper function to get indices
        const auto output_indices = linear_index_to_indices(idx, strides, new_sizes_vec);

        // Compute input indices and validate
        bool valid = true;
        std::vector<std::size_t> input_indices(dim);
        for (std::size_t i = 0; i < dim; ++i)
        {
            input_indices[i] = start[i] + output_indices[i];
            if (input_indices[i] >= input.getSize(i))
            {
                valid = false;
                break;
            }
        }

        // Assign value
        output.set(idx) = valid ? 
            input.at_without_boundary_check(input_indices) : 
            default_element;
    }

    return output;
}

linear_index_to_indices template function implementation

//  linear_index_to_indices template function implementation
//  Helper function to convert linear index to N-dimensional indices
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    const std::size_t linearIndex,
    const Sizes& strides,
    const Sizes& sizes)
{
    const std::size_t dim = std::ranges::size(sizes);
    std::vector<std::size_t> indices(dim);
    for (std::size_t i = 0; i < dim; ++i) {
        indices[i] = (linearIndex / strides[i]) % sizes[i];
    }
    return indices;
}

//  linear_index_to_indices template function implementation
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    std::size_t linear_idx,
    const Sizes& sizes
)
{
    std::vector<std::size_t> indices;
    std::size_t stride = 1;
    for (auto& s : sizes) {
        indices.emplace_back((linear_idx / stride) % s);
        stride *= s;
    }
    return indices;
}

N-dimensional bilateral_filter template function implementation

//  N-dimensional bilateral_filter template function implementation
//  bilateral_filter template function implementation
template<typename ElementT, class ExecutionPolicy, class Fr, class Gs, std::integral SizeT = std::size_t>
requires((std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) and
         (std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    const Fr fr,
    const Gs gs,
    const BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    const ElementT value_for_constant_padding = ElementT{}
    )
{
    if (input.getDimensionality() != 2)
    {
        throw std::runtime_error("Unsupported dimension!");
    }
    return windowed_filter(
        std::forward<ExecutionPolicy>(execution_policy),
        input,
        window_size,
        [&](auto&& window) {
            return bilateral_filter_detail(window, fr, gs);
        },
        boundaryCondition,
        value_for_constant_padding
    );
}

//  bilateral_filter template function implementation
template<class ExecutionPolicy, typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
requires ((std::same_as<ElementT, RGB>) ||
          (std::same_as<ElementT, RGB_DOUBLE>) ||
          (std::same_as<ElementT, HSV>) ||
          (is_MultiChannel<ElementT>::value))
constexpr static auto bilateral_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return apply_each(
        input,
        [&](auto&& each_plane)
        {
            return bilateral_filter(
                std::forward<ExecutionPolicy>(execution_policy),
                each_plane,
                window_size,
                fr,
                gs,
                boundaryCondition
                );
        });
}

//  bilateral_filter template function implementation (without Execution Policy)
template<typename ElementT, class Fr, class Gs, std::integral SizeT = std::size_t>
constexpr static auto bilateral_filter(
    const Image<ElementT>& input,
    const SizeT window_size,
    Fr fr,
    Gs gs,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror
)
{
    return bilateral_filter(
        std::execution::seq,
        input,
        window_size,
        fr,
        gs,
        boundaryCondition
    );
}

bilateral_filter_detail template function implementation

//  bilateral_filter_detail template function implementation
template<typename ElementT, class Fr, class Gs, class FloatingType = double>
requires(std::invocable<Fr, ElementT> && std::invocable<Gs, std::size_t>)
constexpr static auto bilateral_filter_detail(
    const Image<ElementT>& input,
    Fr fr,
    Gs gs
)
{
    const ElementT center_pixel = get_center_pixel(std::execution::seq, input);
    auto center_location = get_center_location(std::execution::seq, input);
    FloatingType sum{}, weight_sum{};
    const std::size_t total_elements = input.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());
        auto fr_result = std::invoke(
            fr,
            std::abs(input.at_without_boundary_check(indices) - center_pixel)
        );
        auto location_difference = difference(std::execution::seq, indices, center_location);
        auto gs_result = std::invoke(
            gs,
            std::reduce(std::ranges::cbegin(location_difference), std::ranges::cend(location_difference))
        );
        sum += input.at_without_boundary_check(indices) * fr_result * gs_result;
        weight_sum += fr_result * gs_result;
    }

    if (weight_sum == 0)
    {
        return static_cast<ElementT>(sum);
    }
    return static_cast<ElementT>(sum / weight_sum);
}

windowed_filter template function implementation

//  windowed_filter template function implementation
template<class ElementT, class ExecutionPolicy, class Filter, class SizeT = std::size_t>
requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
constexpr static auto windowed_filter(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    SizeT window_size,
    Filter filter,
    BoundaryCondition boundaryCondition = BoundaryCondition::mirror,
    ElementT value_for_constant_padding = ElementT{},
    std::function<void(double)> progress_callback = nullptr
)
{
    const std::size_t dim = input.getDimensionality();

    auto padded_image = input;
    if (dim == 2)                   //  padding algorithm supported
    {
        padded_image = generate_padded_image(
            std::forward<ExecutionPolicy>(execution_policy),
            input,
            window_size,
            window_size,
            boundaryCondition,
            value_for_constant_padding
        );
    }
    // Iterate over all elements in the original image
    const std::size_t total_elements = input.count();
    #ifndef USE_OPENMP
    auto indices_view = std::views::iota(std::size_t{ 0 }, total_elements);
    // Progress tracking variables
    std::atomic<std::size_t> processed_count = 0;
    std::atomic<std::size_t> next_report = 0;
    std::mutex report_mutex;
    const std::size_t report_step = std::max<std::size_t>(1, total_elements / 100);
    auto process_element = [&](auto&& idx) {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);

        // Progress reporting
        if (progress_callback) {
            auto count = ++processed_count;
            if (count >= next_report.load(std::memory_order_relaxed)) {
                std::lock_guard lock(report_mutex);
                if (count >= next_report.load(std::memory_order_relaxed)) {
                    double progress = static_cast<double>(count) / total_elements;
                    progress_callback(progress);
                    next_report.store(count + report_step, std::memory_order_relaxed);
                }
            }
        }
        };
    std::for_each(
        std::forward<ExecutionPolicy>(execution_policy),
        indices_view.begin(),
        indices_view.end(),
        process_element
    );
    #else
    #pragma omp parallel for
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Convert linear index to N-D indices (original image)
        auto indices = linear_index_to_indices(idx, input.getSize());

        // Convert to padded indices
        std::vector<std::size_t> padded_indices;
        for (auto& i : indices) {
            padded_indices.emplace_back(i + window_size);
        }

        // Extract window
        std::vector<std::size_t> window_sizes(input.getDimensionality(), window_size);
        auto window = subimage(
            std::forward<ExecutionPolicy>(execution_policy),
            padded_image,
            window_sizes,
            padded_indices
        );

        // Apply filter and store result
        padded_image.at_without_boundary_check(padded_indices) =
            std::invoke(filter, window);
    }
    #endif
    // Crop padded image to original size
    std::vector<std::size_t> original_sizes = input.getSize();
    std::vector<std::size_t> centers;
    for (auto& s : padded_image.getSize())
    {
        centers.emplace_back(s / 2);  // Center of padded image
    }
    auto output = subimage(std::forward<ExecutionPolicy>(execution_policy), padded_image, original_sizes, centers);
    return output;
}

subimage template function implementation for N dimensional image

//  subimage template function implementation for N dimensional image
template<class ExecutionPolicy, typename ElementT, std::ranges::input_range Sizes, std::ranges::input_range Centers>
requires (
    (std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) &&
    (std::same_as<std::ranges::range_value_t<Sizes>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Sizes>, int>) &&
    (std::same_as<std::ranges::range_value_t<Centers>, std::size_t> ||
     std::same_as<std::ranges::range_value_t<Centers>, int>)
)
constexpr static auto subimage(
    ExecutionPolicy&& execution_policy,
    const Image<ElementT>& input,
    const Sizes& new_sizes,
    const Centers& centers,
    const ElementT default_element = ElementT{})
{
    const std::size_t dim = input.getDimensionality();
    if (new_sizes.size() != dim || centers.size() != dim)
    {
        throw std::runtime_error("subimage: new_sizes and centers must match input dimensionality");
    }

    // Compute start indices for each dimension
    std::vector<std::size_t> start(dim);
    for (std::size_t i = 0; i < dim; ++i)
    {
        start[i] = centers[i] - static_cast<std::size_t>(std::floor(static_cast<double>(new_sizes[i]) / 2.0));
    }

    Image<ElementT> output(new_sizes);

    std::vector<std::size_t> new_sizes_vec;
    new_sizes_vec.resize(std::ranges::size(new_sizes));

    std::transform(
        std::forward<ExecutionPolicy>(execution_policy),
        std::ranges::cbegin(new_sizes),
        std::ranges::cend(new_sizes),
        std::ranges::begin(new_sizes_vec),
        [](auto&& element) { return static_cast<std::size_t>(element); }
    );

    // Precompute strides using new_sizes_vec
    std::vector<std::size_t> strides(dim, 1);
    for (std::size_t i = 1; i < dim; ++i)
    {
        strides[i] = strides[i - 1] * new_sizes_vec[i - 1];
    }

    // Iterate through all output elements
    const std::size_t total_elements = output.count();
    for (std::size_t idx = 0; idx < total_elements; ++idx)
    {
        // Use helper function to get indices
        const auto output_indices = linear_index_to_indices(idx, strides, new_sizes_vec);

        // Compute input indices and validate
        bool valid = true;
        std::vector<std::size_t> input_indices(dim);
        for (std::size_t i = 0; i < dim; ++i)
        {
            input_indices[i] = start[i] + output_indices[i];
            if (input_indices[i] >= input.getSize(i))
            {
                valid = false;
                break;
            }
        }

        // Assign value
        output.set(idx) = valid ? 
            input.at_without_boundary_check(input_indices) : 
            default_element;
    }

    return output;
}

linear_index_to_indices template function implementation

//  linear_index_to_indices template function implementation
//  Helper function to convert linear index to N-dimensional indices
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    const std::size_t linearIndex,
    const Sizes& strides,
    const Sizes& sizes)
{
    const std::size_t dim = std::ranges::size(sizes);
    std::vector<std::size_t> indices(dim);
    for (std::size_t i = 0; i < dim; ++i) {
        indices[i] = (linearIndex / strides[i]) % sizes[i];
    }
    return indices;
}

//  linear_index_to_indices template function implementation
template<std::ranges::input_range Sizes>
requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
constexpr static std::vector<std::size_t> linear_index_to_indices(
    std::size_t linear_idx,
    const Sizes& sizes
)
{
    std::vector<std::size_t> indices;
    std::size_t stride = 1;
    for (auto& s : sizes) {
        indices.emplace_back((linear_idx / stride) % s);
        stride *= s;
    }
    return indices;
}