Add linear and median filter operations

cmake/discover_tests.cmake

@@ -42,7 +42,7 @@ function (discover_tests DIR)
 
         string(FIND ${cuda_source} "npp" is_npp)    	
 
-        target_link_libraries(${cuda_target} PRIVATE CUDA::nppc CUDA::nppial CUDA::nppidei CUDA::nppig) 								              
+        target_link_libraries(${cuda_target} PRIVATE CUDA::nppc CUDA::nppial CUDA::nppidei CUDA::nppig CUDA::nppim CUDA::nppif) 	             
 
     endforeach()
 endfunction()

include/fast_npp.h

@@ -23,6 +23,8 @@
 #include <fused_kernel/algorithms/image_processing/resize.h>
 #include <fused_kernel/algorithms/basic_ops/vector_ops.h>
 #include <fused_kernel/algorithms/basic_ops/arithmetic.h>
+#include <fused_kernel/algorithms/image_processing/linear_filter.h>
+#include <fused_kernel/algorithms/image_processing/median_filter.h>
 #include <fused_kernel/algorithms/basic_ops/bitwise.h>
 #include <fused_kernel/algorithms/basic_ops/math.h>
 #include <fused_kernel/algorithms/basic_ops/memory_operations.h>
@@ -35,6 +37,45 @@
 
 namespace fastNPP {
 
+    // ===== Linear filters =====
+    // FilterBoxBorder: mean filter over a mask. Matches nppiFilterBoxBorder
+    // (sum-then-divide, truncated) with NPP_BORDER_REPLICATE.
+#define FASTNPP_DEFINE_BOX(NPPNAME, T)                                              \
+    inline auto NPPNAME(const fk::Ptr2D<T>& pSrc, int nMaskWidth, int nMaskHeight,  \
+                        int nAnchorX, int nAnchorY) {                               \
+        return fk::BoxFilter<fk::ND::_2D, T, fk::FilterBorder::REPLICATE>::build(    \
+            pSrc, nMaskWidth, nMaskHeight, nAnchorX, nAnchorY);                     \
+    }
+    FASTNPP_DEFINE_BOX(FilterBoxBorder_8u_C1R_Ctx,  uchar)
+    FASTNPP_DEFINE_BOX(FilterBoxBorder_8u_C3R_Ctx,  uchar3)
+    FASTNPP_DEFINE_BOX(FilterBoxBorder_16u_C1R_Ctx, ushort)
+    FASTNPP_DEFINE_BOX(FilterBoxBorder_32f_C1R_Ctx, float)
+
+    // FilterBorder: general convolution with a float coefficient kernel.
+    // Matches nppiFilterBorder_32f with NPP_BORDER_REPLICATE.
+#define FASTNPP_DEFINE_CONV(NPPNAME, T)                                             \
+    inline auto NPPNAME(const fk::Ptr2D<T>& pSrc, const fk::Ptr2D<float>& pKernel,  \
+                        int nKernelWidth, int nKernelHeight, int nAnchorX, int nAnchorY) { \
+        return fk::LinearFilter<fk::ND::_2D, T, fk::FilterBorder::REPLICATE,         \
+                                fk::FilterRounding::NEAREST>::build(                 \
+            pSrc, pKernel, nKernelWidth, nKernelHeight, nAnchorX, nAnchorY);        \
+    }
+    FASTNPP_DEFINE_CONV(FilterBorder_32f_C1R_Ctx, float)
+    FASTNPP_DEFINE_CONV(FilterBorder_32f_C3R_Ctx, float3)
+
+    // FilterMedianBorder: per-channel window median. Matches nppiFilterMedianBorder
+    // with NPP_BORDER_REPLICATE. (Note: FKL computes the median directly, so the
+    // NPP scratch-buffer argument is not needed.)
+#define FASTNPP_DEFINE_MEDIAN(NPPNAME, T)                                          \
+    inline auto NPPNAME(const fk::Ptr2D<T>& pSrc, int nMaskWidth, int nMaskHeight, \
+                        int nAnchorX, int nAnchorY) {                              \
+        return fk::MedianFilter<fk::ND::_2D, T, 49, fk::FilterBorder::REPLICATE>::build( \
+            pSrc, nMaskWidth, nMaskHeight, nAnchorX, nAnchorY);                    \
+    }
+    FASTNPP_DEFINE_MEDIAN(FilterMedianBorder_8u_C1R_Ctx,  uchar)
+    FASTNPP_DEFINE_MEDIAN(FilterMedianBorder_8u_C3R_Ctx,  uchar3)
+    FASTNPP_DEFINE_MEDIAN(FilterMedianBorder_16u_C1R_Ctx, ushort)
+    FASTNPP_DEFINE_MEDIAN(FilterMedianBorder_32f_C1R_Ctx, float)
     // ===== AbsDiff with constant: |src - C| =====
     constexpr inline auto AbsDiffC_8u_C1R_Ctx(const uchar& nConstant) {
         return fk::AbsDiff<uchar>::build(nConstant);

tests/filtering/fastNPP_filter_test.cu

@@ -0,0 +1,154 @@
+/* Copyright 2025 Oscar Amoros Huguet
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License. */
+
+// Validates FastNPP linear filters (FilterBoxBorder, FilterBorder convolution)
+// against NVIDIA NPP with NPP_BORDER_REPLICATE.
+
+#ifdef WIN32
+#include <tests/main.h>
+#endif
+
+#include <fast_npp.h>
+#include <cuda_runtime.h>
+#include <vector>
+#include <cstdio>
+#include <random>
+
+namespace {
+
+NppStreamContext makeCtx() {
+    NppStreamContext c{};
+    c.hStream = 0;
+    cudaGetDevice(&c.nCudaDeviceId);
+    cudaDeviceProp p{};
+    cudaGetDeviceProperties(&p, c.nCudaDeviceId);
+    c.nMultiProcessorCount = p.multiProcessorCount;
+    c.nMaxThreadsPerMultiProcessor = p.maxThreadsPerMultiProcessor;
+    c.nMaxThreadsPerBlock = p.maxThreadsPerBlock;
+    c.nSharedMemPerBlock = p.sharedMemPerBlock;
+    c.nCudaDevAttrComputeCapabilityMajor = p.major;
+    c.nCudaDevAttrComputeCapabilityMinor = p.minor;
+    cudaStreamGetFlags(c.hStream, &c.nStreamFlags);
+    return c;
+}
+
+int box8u(int mW, int mH, int aX, int aY) {
+    const int W = 128, H = 96;
+    const size_t N = (size_t)W * H;
+    std::vector<Npp8u> h(N), ref(N), fkl(N);
+    std::mt19937 rng(7); std::uniform_int_distribution<int> di(0, 255);
+    for (size_t i = 0; i < N; ++i) h[i] = (Npp8u)di(rng);
+    fk::Ptr2D<uchar> s(W, H), out(W, H);
+    const int pitch = (int)s.ptr().dims.pitch, rb = W;
+    Npp8u *ds, *dd;
+    cudaMalloc(&ds, (size_t)pitch * H); cudaMalloc(&dd, (size_t)pitch * H);
+    cudaMemcpy2D(ds, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    cudaMemset(dd, 0, (size_t)pitch * H);
+    NppiSize ms{mW, mH}; NppiPoint an{aX, aY}; NppiSize ss{W, H}; NppiPoint so{0, 0};
+    nppiFilterBoxBorder_8u_C1R_Ctx(ds, pitch, ss, so, dd, pitch, {W, H}, ms, an, NPP_BORDER_REPLICATE, makeCtx());
+    cudaDeviceSynchronize();
+    cudaMemcpy2D(ref.data(), rb, dd, pitch, rb, H, cudaMemcpyDeviceToHost);
+    cudaMemcpy2D(s.ptr().data, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    fk::Stream st;
+    fk::executeOperations<fk::TransformDPP<>>(st,
+        fastNPP::FilterBoxBorder_8u_C1R_Ctx(s, mW, mH, aX, aY),
+        fk::PerThreadWrite<fk::ND::_2D, uchar>::build(out));
+    st.sync();
+    cudaMemcpy2D(fkl.data(), rb, out.ptr().data, pitch, rb, H, cudaMemcpyDeviceToHost);
+    int bad = 0; for (size_t i = 0; i < N; ++i) if (ref[i] != fkl[i]) ++bad;
+    printf("[%s] FilterBoxBorder_8u_C1R %dx%d a(%d,%d) mismatches=%d/%zu\n",
+           bad ? "FAIL" : "PASS", mW, mH, aX, aY, bad, N);
+    cudaFree(ds); cudaFree(dd);
+    return bad;
+}
+
+int conv32f() {
+    const int W = 128, H = 96, kW = 3, kH = 3, aX = 1, aY = 1;
+    const size_t N = (size_t)W * H;
+    std::vector<float> h(N), ref(N), fkl(N);
+    std::mt19937 rng(9); std::uniform_real_distribution<float> dr(-10, 10);
+    for (size_t i = 0; i < N; ++i) h[i] = dr(rng);
+    std::vector<float> hk = {0, -1, 0, -1, 5, -1, 0, -1, 0};
+    fk::Ptr2D<float> s(W, H), out(W, H), ker(kW, kH);
+    const int pitch = (int)s.ptr().dims.pitch, rb = W * sizeof(float), kp = (int)ker.ptr().dims.pitch;
+    Npp32f *ds, *dd, *dk;
+    cudaMalloc(&ds, (size_t)pitch * H); cudaMalloc(&dd, (size_t)pitch * H); cudaMalloc(&dk, kW * kH * sizeof(float));
+    cudaMemcpy2D(ds, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    cudaMemcpy(dk, hk.data(), kW * kH * sizeof(float), cudaMemcpyHostToDevice);
+    cudaMemset(dd, 0, (size_t)pitch * H);
+    NppiSize ks{kW, kH}; NppiPoint an{aX, aY}; NppiSize ss{W, H}; NppiPoint so{0, 0};
+    nppiFilterBorder_32f_C1R_Ctx(ds, pitch, ss, so, dd, pitch, {W, H}, dk, ks, an, NPP_BORDER_REPLICATE, makeCtx());
+    cudaDeviceSynchronize();
+    cudaMemcpy2D(ref.data(), rb, dd, pitch, rb, H, cudaMemcpyDeviceToHost);
+    cudaMemcpy2D(s.ptr().data, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    cudaMemcpy2D(ker.ptr().data, kp, hk.data(), kW * sizeof(float), kW * sizeof(float), kH, cudaMemcpyHostToDevice);
+    fk::Stream st;
+    fk::executeOperations<fk::TransformDPP<>>(st,
+        fastNPP::FilterBorder_32f_C1R_Ctx(s, ker, kW, kH, aX, aY),
+        fk::PerThreadWrite<fk::ND::_2D, float>::build(out));
+    st.sync();
+    cudaMemcpy2D(fkl.data(), rb, out.ptr().data, pitch, rb, H, cudaMemcpyDeviceToHost);
+    int bad = 0; double maxd = 0;
+    for (size_t i = 0; i < N; ++i) { double d = fabs((double)ref[i] - (double)fkl[i]); if (d > maxd) maxd = d; if (d > 1e-2) ++bad; }
+    printf("[%s] FilterBorder_32f_C1R 3x3 conv mismatches=%d/%zu maxdiff=%.3e\n", bad ? "FAIL" : "PASS", bad, N, maxd);
+    cudaFree(ds); cudaFree(dd); cudaFree(dk);
+    return bad;
+}
+
+int medianTest(int mW, int mH, int aX, int aY) {
+    const int W = 128, H = 96;
+    const size_t N = (size_t)W * H;
+    std::vector<Npp8u> h(N), ref(N), fkl(N);
+    std::mt19937 rng(5); std::uniform_int_distribution<int> di(0, 255);
+    for (size_t i = 0; i < N; ++i) h[i] = (Npp8u)di(rng);
+    fk::Ptr2D<uchar> s(W, H), out(W, H);
+    const int pitch = (int)s.ptr().dims.pitch, rb = W;
+    Npp8u *ds, *dd;
+    cudaMalloc(&ds, (size_t)pitch * H); cudaMalloc(&dd, (size_t)pitch * H);
+    cudaMemcpy2D(ds, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    cudaMemset(dd, 0, (size_t)pitch * H);
+    NppiSize ms{mW, mH}; NppiPoint an{aX, aY}; NppiSize ss{W, H}; NppiPoint so{0, 0};
+    Npp32u bufSz = 0;
+    nppiFilterMedianBorderGetBufferSize_8u_C1R_Ctx({W, H}, ms, &bufSz, NPP_BORDER_REPLICATE, makeCtx());
+    Npp8u* dbuf = nullptr; if (bufSz) cudaMalloc(&dbuf, bufSz);
+    nppiFilterMedianBorder_8u_C1R_Ctx(ds, pitch, ss, so, dd, pitch, {W, H}, ms, an, dbuf, NPP_BORDER_REPLICATE, makeCtx());
+    cudaDeviceSynchronize();
+    cudaMemcpy2D(ref.data(), rb, dd, pitch, rb, H, cudaMemcpyDeviceToHost);
+    cudaMemcpy2D(s.ptr().data, pitch, h.data(), rb, rb, H, cudaMemcpyHostToDevice);
+    fk::Stream st;
+    fk::executeOperations<fk::TransformDPP<>>(st,
+        fastNPP::FilterMedianBorder_8u_C1R_Ctx(s, mW, mH, aX, aY),
+        fk::PerThreadWrite<fk::ND::_2D, uchar>::build(out));
+    st.sync();
+    cudaMemcpy2D(fkl.data(), rb, out.ptr().data, pitch, rb, H, cudaMemcpyDeviceToHost);
+    int bad = 0; for (size_t i = 0; i < N; ++i) if (ref[i] != fkl[i]) ++bad;
+    printf("[%s] FilterMedianBorder_8u_C1R %dx%d a(%d,%d) mismatches=%d/%zu\n",
+           bad ? "FAIL" : "PASS", mW, mH, aX, aY, bad, N);
+    if (dbuf) cudaFree(dbuf);
+    cudaFree(ds); cudaFree(dd);
+    return bad;
+}
+
+} // namespace
+
+int launch() {
+    int bad = 0;
+    bad += box8u(3, 3, 1, 1);
+    bad += box8u(5, 5, 2, 2);
+    bad += conv32f();
+    bad += medianTest(3, 3, 1, 1);
+    bad += medianTest(5, 5, 2, 2);
+    printf("%s\n", bad == 0 ? "ALL PASS" : "FAILURES DETECTED");
+    return bad == 0 ? 0 : 1;
+}