From 09f7da11b6ef7a34918954b7dc70237144541315 Mon Sep 17 00:00:00 2001 From: johnnynunez Date: Sun, 28 Jun 2026 18:35:05 +0200 Subject: [PATCH] feat(image-processing): add composable neighborhood DPPs --- .../image_processing/convolution_fast.h | 300 ++++++++++++++ .../image_processing/image_processing.h | 3 + .../algorithms/image_processing/median_fast.h | 388 ++++++++++++++++++ .../image_processing/morphology_fast.h | 278 +++++++++++++ .../test_convolution_fast_dpp.h | 175 ++++++++ tests/image_processing/test_median_fast_dpp.h | 160 ++++++++ .../test_morphology_fast_dpp.h | 170 ++++++++ 7 files changed, 1474 insertions(+) create mode 100644 include/fused_kernel/algorithms/image_processing/convolution_fast.h create mode 100644 include/fused_kernel/algorithms/image_processing/median_fast.h create mode 100644 include/fused_kernel/algorithms/image_processing/morphology_fast.h create mode 100644 tests/image_processing/test_convolution_fast_dpp.h create mode 100644 tests/image_processing/test_median_fast_dpp.h create mode 100644 tests/image_processing/test_morphology_fast_dpp.h diff --git a/include/fused_kernel/algorithms/image_processing/convolution_fast.h b/include/fused_kernel/algorithms/image_processing/convolution_fast.h new file mode 100644 index 00000000..f4928494 --- /dev/null +++ b/include/fused_kernel/algorithms/image_processing/convolution_fast.h @@ -0,0 +1,300 @@ +/* Copyright 2026 Oscar Amoros Huguet, Johnny Nunez + + 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. */ + +#ifndef FK_CONVOLUTION_FAST_H +#define FK_CONVOLUTION_FAST_H + +#include +#include +#include + +#include + +namespace fk { + +constexpr int FK_CONV_MAX_KERNEL_SIDE = 15; +constexpr int FK_CONV_MAX_COEFFICIENTS = + FK_CONV_MAX_KERNEL_SIDE * FK_CONV_MAX_KERNEL_SIDE; + +struct ConvQuadDetails { + int width; + int height; + int kernelWidth; + int kernelHeight; + int anchorX; + int anchorY; + float coefficients[FK_CONV_MAX_COEFFICIENTS]; + + template + FK_HOST_DEVICE_CNST bool validFor() const { + const int effectiveWidth = + STATIC_KW > 0 ? STATIC_KW : kernelWidth; + const int effectiveHeight = + STATIC_KH > 0 ? STATIC_KH : kernelHeight; + return width > 0 && height > 0 && + effectiveWidth > 0 && effectiveHeight > 0 && + effectiveWidth <= FK_CONV_MAX_KERNEL_SIDE && + effectiveHeight <= FK_CONV_MAX_KERNEL_SIDE && + anchorX >= 0 && anchorX < effectiveWidth && + anchorY >= 0 && anchorY < effectiveHeight; + } +}; + +template +struct ConvQuadDPP; + +template +struct ConvQuadDPP { +private: + using SelfType = ConvQuadDPP; + + FK_HOST_FUSE int clamp(const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + } + +public: + FK_STATIC_STRUCT(ConvQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::CPU; + + FK_HOST_DEVICE_FUSE bool accepts(const ConvQuadDetails& details) { + return details.template validFor(); + } + + template + FK_HOST_FUSE void exec(const ConvQuadDetails& details, + const InIOp& input, + const ComputeIOps& compute, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "ConvQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "ConvQuadDPP requires a Write IOp"); + static_assert(std::is_same_v, + "ConvQuadDPP currently supports float pixels"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "ConvQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "ConvQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "ConvQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "ConvQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "ConvQuadDPP sizes are both static or both runtime"); + if (!accepts(details)) return; + + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + const auto& multiply = get<0>(compute); + const auto& combine = get<1>(compute); + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + for (int y = 0; y < details.height; ++y) { + for (int x = 0; x < details.width; ++x) { + float value = 0.f; + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + const float product = make_tuple( + details.coefficients[ky * kernelWidth + kx], + source(x + kx - details.anchorX, + y + ky - details.anchorY)) | multiply; + value = make_tuple(value, product) | combine; + } + } + OutIOp::Operation::exec(Point{x, y, 0}, value, output); + } + } + } +}; + +#if defined(__NVCC__) +template +struct ConvQuadDPP { +private: + using SelfType = ConvQuadDPP< + ParArch::GPU_NVIDIA, T, EX, EY, KW, KH>; + static constexpr int MAX_KERNEL_WIDTH = + KW > 0 ? KW : FK_CONV_MAX_KERNEL_SIDE; + static constexpr int MAX_SPAN = EX + MAX_KERNEL_WIDTH - 1; + +public: + FK_STATIC_STRUCT(ConvQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::GPU_NVIDIA; + static constexpr int BLOCK_THREADS = 256; + + struct LaunchConfig { + unsigned int blocks; + unsigned int threads; + }; + + FK_HOST_DEVICE_FUSE bool accepts(const ConvQuadDetails& details) { + return details.template validFor(); + } + + FK_HOST_FUSE LaunchConfig launchConfig(const ConvQuadDetails& details) { + const int columns = (details.width + EX - 1) / EX; + const int rows = (details.height + EY - 1) / EY; + const int workItems = columns * rows; + return {static_cast( + (workItems + BLOCK_THREADS - 1) / BLOCK_THREADS), + static_cast(BLOCK_THREADS)}; + } + + template + FK_DEVICE_FUSE void exec(const ConvQuadDetails& details, + const InIOp& input, + const ComputeIOps& compute, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "ConvQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "ConvQuadDPP requires a Write IOp"); + static_assert(std::is_same_v, + "ConvQuadDPP currently supports float pixels"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "ConvQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "ConvQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "ConvQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "ConvQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "ConvQuadDPP sizes are both static or both runtime"); +#if defined(__CUDA_ARCH__) + if (!accepts(details)) return; + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + const auto& multiply = get<0>(compute); + const auto& combine = get<1>(compute); + const int columns = (details.width + EX - 1) / EX; + const int workItem = static_cast( + blockIdx.x * blockDim.x + threadIdx.x); + const int x0 = (workItem % columns) * EX; + const int y0 = (workItem / columns) * EY; + if (x0 >= details.width || y0 >= details.height) return; + + const int span = EX + kernelWidth - 1; + if (span > MAX_SPAN) return; + const int rowsNeeded = EY + kernelHeight - 1; + const int firstColumn = x0 - details.anchorX; + auto clamp = [](const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + }; + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + float accumulated[EY][EX]; + #pragma unroll + for (int localY = 0; localY < EY; ++localY) { + #pragma unroll + for (int localX = 0; localX < EX; ++localX) { + accumulated[localY][localX] = 0.f; + } + } + + #pragma unroll + for (int row = 0; row < rowsNeeded; ++row) { + float sourceRow[MAX_SPAN]; + #pragma unroll + for (int index = 0; index < span; ++index) { + sourceRow[index] = source( + firstColumn + index, + y0 - details.anchorY + row); + } + #pragma unroll + for (int localY = 0; localY < EY; ++localY) { + const int ky = row - localY; + if (ky < 0 || ky >= kernelHeight || + y0 + localY >= details.height) continue; + #pragma unroll + for (int kx = 0; kx < kernelWidth; ++kx) { + const float coefficient = + details.coefficients[ky * kernelWidth + kx]; + #pragma unroll + for (int localX = 0; localX < EX; ++localX) { + const float product = make_tuple( + coefficient, + sourceRow[localX + kx]) | multiply; + accumulated[localY][localX] = make_tuple( + accumulated[localY][localX], product) | combine; + } + } + } + } + + #pragma unroll + for (int localY = 0; localY < EY; ++localY) { + const int y = y0 + localY; + if (y >= details.height) break; + #pragma unroll + for (int localX = 0; localX < EX; ++localX) { + const int x = x0 + localX; + if (x >= details.width) break; + OutIOp::Operation::exec( + Point{x, y, 0}, accumulated[localY][localX], output); + } + } +#endif // defined(__CUDA_ARCH__) + } +}; + +template +__global__ void launchConvQuadDPP_Kernel( + const __grid_constant__ ConvQuadDetails details, + const __grid_constant__ IOps... iOps) { + DPP::exec(details, iOps...); +} + +template +FK_HOST_FUSE void executeConvQuad( + Stream_& stream, + const ConvQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::GPU_NVIDIA, + "GPU stream requires the NVIDIA ConvQuadDPP specialization"); + if (!DPP::accepts(details)) return; + const auto launch = DPP::launchConfig(details); + launchConvQuadDPP_Kernel + <<>>( + details, iOps...); + gpuErrchk(cudaGetLastError()); +} +#endif // defined(__NVCC__) + +template +FK_HOST_FUSE void executeConvQuad( + Stream_&, + const ConvQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::CPU, + "CPU stream requires the CPU ConvQuadDPP specialization"); + DPP::exec(details, iOps...); +} + +} // namespace fk + +#endif // FK_CONVOLUTION_FAST_H diff --git a/include/fused_kernel/algorithms/image_processing/image_processing.h b/include/fused_kernel/algorithms/image_processing/image_processing.h index 42d1a4d7..f4915305 100644 --- a/include/fused_kernel/algorithms/image_processing/image_processing.h +++ b/include/fused_kernel/algorithms/image_processing/image_processing.h @@ -24,5 +24,8 @@ #include #include #include +#include +#include +#include #endif // FK_IMAGE_PROCESSING diff --git a/include/fused_kernel/algorithms/image_processing/median_fast.h b/include/fused_kernel/algorithms/image_processing/median_fast.h new file mode 100644 index 00000000..9a70a47c --- /dev/null +++ b/include/fused_kernel/algorithms/image_processing/median_fast.h @@ -0,0 +1,388 @@ +/* Copyright 2026 Oscar Amoros Huguet, Johnny Nunez + + 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. */ + +#ifndef FK_MEDIAN_FAST_H +#define FK_MEDIAN_FAST_H + +#include +#include +#include +#include +#include + +#include + +namespace fk { + +constexpr int FK_MEDIAN_MAX_KERNEL_SIDE = 7; +constexpr int FK_MEDIAN_SORT_SIZE = 64; + +struct MedianQuadDetails { + int width; + int height; + int kernelWidth; + int kernelHeight; + int anchorX; + int anchorY; + + template + FK_HOST_DEVICE_CNST bool validFor() const { + const int effectiveWidth = + STATIC_KW > 0 ? STATIC_KW : kernelWidth; + const int effectiveHeight = + STATIC_KH > 0 ? STATIC_KH : kernelHeight; + return width > 0 && height > 0 && + effectiveWidth > 0 && effectiveHeight > 0 && + effectiveWidth <= FK_MEDIAN_MAX_KERNEL_SIDE && + effectiveHeight <= FK_MEDIAN_MAX_KERNEL_SIDE && + effectiveWidth * effectiveHeight <= FK_MEDIAN_SORT_SIZE && + anchorX >= 0 && anchorX < effectiveWidth && + anchorY >= 0 && anchorY < effectiveHeight; + } +}; + +namespace median_detail { +template +FK_HOST_DEVICE_FUSE void compareSwap( + T* values, const int first, const int second, + const MinIOp& minimum, const MaxIOp& maximum, + const bool ascending = true) { + const T left = values[first]; + const T right = values[second]; + const T low = make_tuple(left, right) | minimum; + const T high = make_tuple(left, right) | maximum; + values[first] = ascending ? low : high; + values[second] = ascending ? high : low; +} + +template +FK_HOST_DEVICE_FUSE T medianNet9( + T* values, const MinIOp& minimum, const MaxIOp& maximum) { +#define FK_MEDIAN_SWAP(a, b) \ + compareSwap(values, a, b, minimum, maximum) + FK_MEDIAN_SWAP(0,1); FK_MEDIAN_SWAP(2,3); FK_MEDIAN_SWAP(4,5); + FK_MEDIAN_SWAP(6,7); FK_MEDIAN_SWAP(0,2); FK_MEDIAN_SWAP(1,3); + FK_MEDIAN_SWAP(4,6); FK_MEDIAN_SWAP(5,7); FK_MEDIAN_SWAP(1,2); + FK_MEDIAN_SWAP(5,6); FK_MEDIAN_SWAP(0,4); FK_MEDIAN_SWAP(1,5); + FK_MEDIAN_SWAP(2,6); FK_MEDIAN_SWAP(3,7); FK_MEDIAN_SWAP(2,4); + FK_MEDIAN_SWAP(3,5); FK_MEDIAN_SWAP(1,2); FK_MEDIAN_SWAP(3,4); + FK_MEDIAN_SWAP(5,6); FK_MEDIAN_SWAP(0,8); FK_MEDIAN_SWAP(4,8); + FK_MEDIAN_SWAP(2,4); FK_MEDIAN_SWAP(3,5); FK_MEDIAN_SWAP(6,8); + FK_MEDIAN_SWAP(1,2); FK_MEDIAN_SWAP(3,4); FK_MEDIAN_SWAP(5,6); + FK_MEDIAN_SWAP(7,8); +#undef FK_MEDIAN_SWAP + return values[4]; +} + +template +FK_HOST_DEVICE_FUSE T medianNet25( + T* values, const MinIOp& minimum, const MaxIOp& maximum) { +#define FK_MEDIAN_SWAP(a, b) \ + compareSwap(values, a, b, minimum, maximum) + FK_MEDIAN_SWAP(0,1); FK_MEDIAN_SWAP(2,3); FK_MEDIAN_SWAP(4,5); + FK_MEDIAN_SWAP(6,7); FK_MEDIAN_SWAP(8,9); FK_MEDIAN_SWAP(10,11); + FK_MEDIAN_SWAP(12,13); FK_MEDIAN_SWAP(14,15); FK_MEDIAN_SWAP(16,17); + FK_MEDIAN_SWAP(18,19); FK_MEDIAN_SWAP(20,21); FK_MEDIAN_SWAP(22,23); + FK_MEDIAN_SWAP(0,2); FK_MEDIAN_SWAP(1,3); FK_MEDIAN_SWAP(4,6); + FK_MEDIAN_SWAP(5,7); FK_MEDIAN_SWAP(8,10); FK_MEDIAN_SWAP(9,11); + FK_MEDIAN_SWAP(12,14); FK_MEDIAN_SWAP(13,15); FK_MEDIAN_SWAP(16,18); + FK_MEDIAN_SWAP(17,19); FK_MEDIAN_SWAP(20,22); FK_MEDIAN_SWAP(21,23); + FK_MEDIAN_SWAP(1,2); FK_MEDIAN_SWAP(5,6); FK_MEDIAN_SWAP(9,10); + FK_MEDIAN_SWAP(13,14); FK_MEDIAN_SWAP(17,18); FK_MEDIAN_SWAP(21,22); + FK_MEDIAN_SWAP(0,4); FK_MEDIAN_SWAP(1,5); FK_MEDIAN_SWAP(2,6); + FK_MEDIAN_SWAP(3,7); FK_MEDIAN_SWAP(8,12); FK_MEDIAN_SWAP(9,13); + FK_MEDIAN_SWAP(10,14); FK_MEDIAN_SWAP(11,15); FK_MEDIAN_SWAP(16,20); + FK_MEDIAN_SWAP(17,21); FK_MEDIAN_SWAP(18,22); FK_MEDIAN_SWAP(19,23); + FK_MEDIAN_SWAP(2,4); FK_MEDIAN_SWAP(3,5); FK_MEDIAN_SWAP(10,12); + FK_MEDIAN_SWAP(11,13); FK_MEDIAN_SWAP(18,20); FK_MEDIAN_SWAP(19,21); + FK_MEDIAN_SWAP(1,2); FK_MEDIAN_SWAP(3,4); FK_MEDIAN_SWAP(5,6); + FK_MEDIAN_SWAP(9,10); FK_MEDIAN_SWAP(11,12); FK_MEDIAN_SWAP(13,14); + FK_MEDIAN_SWAP(17,18); FK_MEDIAN_SWAP(19,20); FK_MEDIAN_SWAP(21,22); + FK_MEDIAN_SWAP(0,8); FK_MEDIAN_SWAP(1,9); FK_MEDIAN_SWAP(2,10); + FK_MEDIAN_SWAP(3,11); FK_MEDIAN_SWAP(4,12); FK_MEDIAN_SWAP(5,13); + FK_MEDIAN_SWAP(6,14); FK_MEDIAN_SWAP(7,15); FK_MEDIAN_SWAP(16,24); + FK_MEDIAN_SWAP(4,8); FK_MEDIAN_SWAP(5,9); FK_MEDIAN_SWAP(6,10); + FK_MEDIAN_SWAP(7,11); FK_MEDIAN_SWAP(20,24); FK_MEDIAN_SWAP(2,4); + FK_MEDIAN_SWAP(3,5); FK_MEDIAN_SWAP(6,8); FK_MEDIAN_SWAP(7,9); + FK_MEDIAN_SWAP(10,12); FK_MEDIAN_SWAP(11,13); FK_MEDIAN_SWAP(18,20); + FK_MEDIAN_SWAP(19,21); FK_MEDIAN_SWAP(22,24); FK_MEDIAN_SWAP(1,2); + FK_MEDIAN_SWAP(3,4); FK_MEDIAN_SWAP(5,6); FK_MEDIAN_SWAP(7,8); + FK_MEDIAN_SWAP(9,10); FK_MEDIAN_SWAP(11,12); FK_MEDIAN_SWAP(13,14); + FK_MEDIAN_SWAP(17,18); FK_MEDIAN_SWAP(19,20); FK_MEDIAN_SWAP(21,22); + FK_MEDIAN_SWAP(23,24); FK_MEDIAN_SWAP(0,16); FK_MEDIAN_SWAP(1,17); + FK_MEDIAN_SWAP(2,18); FK_MEDIAN_SWAP(3,19); FK_MEDIAN_SWAP(4,20); + FK_MEDIAN_SWAP(5,21); FK_MEDIAN_SWAP(6,22); FK_MEDIAN_SWAP(7,23); + FK_MEDIAN_SWAP(8,24); FK_MEDIAN_SWAP(8,16); FK_MEDIAN_SWAP(9,17); + FK_MEDIAN_SWAP(10,18); FK_MEDIAN_SWAP(11,19); FK_MEDIAN_SWAP(12,20); + FK_MEDIAN_SWAP(13,21); FK_MEDIAN_SWAP(14,22); FK_MEDIAN_SWAP(15,23); + FK_MEDIAN_SWAP(4,8); FK_MEDIAN_SWAP(5,9); FK_MEDIAN_SWAP(6,10); + FK_MEDIAN_SWAP(7,11); FK_MEDIAN_SWAP(12,16); FK_MEDIAN_SWAP(13,17); + FK_MEDIAN_SWAP(14,18); FK_MEDIAN_SWAP(15,19); FK_MEDIAN_SWAP(20,24); + FK_MEDIAN_SWAP(2,4); FK_MEDIAN_SWAP(3,5); FK_MEDIAN_SWAP(6,8); + FK_MEDIAN_SWAP(7,9); FK_MEDIAN_SWAP(10,12); FK_MEDIAN_SWAP(11,13); + FK_MEDIAN_SWAP(14,16); FK_MEDIAN_SWAP(15,17); FK_MEDIAN_SWAP(18,20); + FK_MEDIAN_SWAP(19,21); FK_MEDIAN_SWAP(22,24); FK_MEDIAN_SWAP(1,2); + FK_MEDIAN_SWAP(3,4); FK_MEDIAN_SWAP(5,6); FK_MEDIAN_SWAP(7,8); + FK_MEDIAN_SWAP(9,10); FK_MEDIAN_SWAP(11,12); FK_MEDIAN_SWAP(13,14); + FK_MEDIAN_SWAP(15,16); FK_MEDIAN_SWAP(17,18); FK_MEDIAN_SWAP(19,20); + FK_MEDIAN_SWAP(21,22); FK_MEDIAN_SWAP(23,24); +#undef FK_MEDIAN_SWAP + return values[12]; +} + +template +FK_HOST_DEVICE_FUSE T medianBitonic64( + T* values, const int count, + const MinIOp& minimum, const MaxIOp& maximum) { + for (int index = count; index < FK_MEDIAN_SORT_SIZE; ++index) { + values[index] = maxValue; + } + for (int size = 2; size <= FK_MEDIAN_SORT_SIZE; size <<= 1) { + for (int stride = size >> 1; stride > 0; stride >>= 1) { + for (int index = 0; index < FK_MEDIAN_SORT_SIZE; ++index) { + const int partner = index ^ stride; + if (partner > index) { + compareSwap(values, index, partner, minimum, maximum, + (index & size) == 0); + } + } + } + } + return values[count / 2]; +} + +template +FK_HOST_DEVICE_FUSE T median( + T* values, const int count, const int kernelWidth, + const int kernelHeight, const ComputeIOps& compute) { + const auto& minimum = get<0>(compute); + const auto& maximum = get<1>(compute); + if (kernelWidth == 3 && kernelHeight == 3) { + return medianNet9(values, minimum, maximum); + } + if (kernelWidth == 5 && kernelHeight == 5) { + return medianNet25(values, minimum, maximum); + } + return medianBitonic64(values, count, minimum, maximum); +} +} // namespace median_detail + +template +struct MedianQuadDPP; + +template +struct MedianQuadDPP { +private: + using SelfType = MedianQuadDPP; + + FK_HOST_FUSE int clamp(const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + } + +public: + FK_STATIC_STRUCT(MedianQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::CPU; + + FK_HOST_DEVICE_FUSE bool accepts(const MedianQuadDetails& details) { + return details.template validFor(); + } + + template + FK_HOST_FUSE void exec(const MedianQuadDetails& details, + const InIOp& input, + const ComputeIOps& compute, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "MedianQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "MedianQuadDPP requires a Write IOp"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "MedianQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "MedianQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "MedianQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "MedianQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "MedianQuadDPP sizes are both static or both runtime"); + static_assert(cn == 1, + "MedianQuadDPP currently supports scalar pixel types"); + if (!accepts(details)) return; + + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + for (int y = 0; y < details.height; ++y) { + for (int x = 0; x < details.width; ++x) { + T values[FK_MEDIAN_SORT_SIZE]; + int count = 0; + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + values[count++] = source( + x + kx - details.anchorX, + y + ky - details.anchorY); + } + } + const T value = median_detail::median( + values, count, kernelWidth, kernelHeight, compute); + OutIOp::Operation::exec(Point{x, y, 0}, value, output); + } + } + } +}; + +#if defined(__NVCC__) +template +struct MedianQuadDPP { +private: + using SelfType = MedianQuadDPP< + ParArch::GPU_NVIDIA, T, EX, EY, KW, KH>; + +public: + FK_STATIC_STRUCT(MedianQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::GPU_NVIDIA; + static constexpr int BLOCK_THREADS = 256; + + struct LaunchConfig { + unsigned int blocks; + unsigned int threads; + }; + + FK_HOST_DEVICE_FUSE bool accepts(const MedianQuadDetails& details) { + return details.template validFor(); + } + + FK_HOST_FUSE LaunchConfig launchConfig(const MedianQuadDetails& details) { + const int columns = (details.width + EX - 1) / EX; + const int rows = (details.height + EY - 1) / EY; + const int workItems = columns * rows; + return {static_cast( + (workItems + BLOCK_THREADS - 1) / BLOCK_THREADS), + static_cast(BLOCK_THREADS)}; + } + + template + FK_DEVICE_FUSE void exec(const MedianQuadDetails& details, + const InIOp& input, + const ComputeIOps& compute, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "MedianQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "MedianQuadDPP requires a Write IOp"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "MedianQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "MedianQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "MedianQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "MedianQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "MedianQuadDPP sizes are both static or both runtime"); + static_assert(cn == 1, + "MedianQuadDPP currently supports scalar pixel types"); +#if defined(__CUDA_ARCH__) + if (!accepts(details)) return; + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + const int columns = (details.width + EX - 1) / EX; + const int workItem = static_cast( + blockIdx.x * blockDim.x + threadIdx.x); + const int x0 = (workItem % columns) * EX; + const int y0 = (workItem / columns) * EY; + if (x0 >= details.width || y0 >= details.height) return; + + auto clamp = [](const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + }; + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + #pragma unroll + for (int localY = 0; localY < EY; ++localY) { + const int y = y0 + localY; + if (y >= details.height) break; + #pragma unroll + for (int localX = 0; localX < EX; ++localX) { + const int x = x0 + localX; + if (x >= details.width) break; + T values[FK_MEDIAN_SORT_SIZE]; + int count = 0; + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + values[count++] = source( + x + kx - details.anchorX, + y + ky - details.anchorY); + } + } + const T value = median_detail::median( + values, count, kernelWidth, kernelHeight, compute); + OutIOp::Operation::exec(Point{x, y, 0}, value, output); + } + } +#endif // defined(__CUDA_ARCH__) + } +}; + +template +__global__ void launchMedianQuadDPP_Kernel( + const __grid_constant__ MedianQuadDetails details, + const __grid_constant__ IOps... iOps) { + DPP::exec(details, iOps...); +} + +template +FK_HOST_FUSE void executeMedianQuad( + Stream_& stream, + const MedianQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::GPU_NVIDIA, + "GPU stream requires the NVIDIA MedianQuadDPP specialization"); + if (!DPP::accepts(details)) return; + const auto launch = DPP::launchConfig(details); + launchMedianQuadDPP_Kernel + <<>>( + details, iOps...); + gpuErrchk(cudaGetLastError()); +} +#endif // defined(__NVCC__) + +template +FK_HOST_FUSE void executeMedianQuad( + Stream_&, + const MedianQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::CPU, + "CPU stream requires the CPU MedianQuadDPP specialization"); + DPP::exec(details, iOps...); +} + +} // namespace fk + +#endif // FK_MEDIAN_FAST_H diff --git a/include/fused_kernel/algorithms/image_processing/morphology_fast.h b/include/fused_kernel/algorithms/image_processing/morphology_fast.h new file mode 100644 index 00000000..7c072dad --- /dev/null +++ b/include/fused_kernel/algorithms/image_processing/morphology_fast.h @@ -0,0 +1,278 @@ +/* Copyright 2026 Oscar Amoros Huguet, Johnny Nunez + + 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. */ + +#ifndef FK_MORPHOLOGY_FAST_H +#define FK_MORPHOLOGY_FAST_H + +#include +#include +#include +#include + +#include + +namespace fk { + +struct MorphQuadDetails { + int width; + int height; + int kernelWidth; + int kernelHeight; + int anchorX; + int anchorY; + + template + FK_HOST_DEVICE_CNST bool validFor( + const int maxRuntimeKernelWidth) const { + const int effectiveWidth = + STATIC_KW > 0 ? STATIC_KW : kernelWidth; + const int effectiveHeight = + STATIC_KH > 0 ? STATIC_KH : kernelHeight; + const bool runtimeWidthSupported = + STATIC_KW > 0 || effectiveWidth <= maxRuntimeKernelWidth; + return width > 0 && height > 0 && + effectiveWidth > 0 && effectiveHeight > 0 && + runtimeWidthSupported && + anchorX >= 0 && anchorX < effectiveWidth && + anchorY >= 0 && anchorY < effectiveHeight; + } +}; + +template +struct MorphQuadDPP; + +template +struct MorphQuadDPP { +private: + using SelfType = MorphQuadDPP; + + FK_HOST_FUSE int clamp(const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + } + +public: + FK_STATIC_STRUCT(MorphQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::CPU; + static constexpr int MAX_RUNTIME_KERNEL_WIDTH = 31; + + FK_HOST_DEVICE_FUSE bool accepts(const MorphQuadDetails& details) { + return details.template validFor( + MAX_RUNTIME_KERNEL_WIDTH); + } + + template + FK_HOST_FUSE void exec(const MorphQuadDetails& details, + const InIOp& input, + const ReduceIOp& reduce, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "MorphQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "MorphQuadDPP requires a Write IOp"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "MorphQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "MorphQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "MorphQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "MorphQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "MorphQuadDPP sizes are both static or both runtime"); + static_assert(cn == 1, + "MorphQuadDPP currently supports scalar pixel types"); + if (!accepts(details)) return; + + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + for (int y = 0; y < details.height; ++y) { + for (int x = 0; x < details.width; ++x) { + T value = source(x - details.anchorX, + y - details.anchorY); + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + if (kx == 0 && ky == 0) continue; + value = make_tuple( + value, + source(x + kx - details.anchorX, + y + ky - details.anchorY)) | reduce; + } + } + OutIOp::Operation::exec(Point{x, y, 0}, value, output); + } + } + } +}; + +#if defined(__NVCC__) +template +struct MorphQuadDPP { +private: + using SelfType = MorphQuadDPP< + ParArch::GPU_NVIDIA, T, EX, EY, KW, KH>; + static constexpr int MAX_RUNTIME_KERNEL_WIDTH = 31; + static constexpr int MAX_KERNEL_WIDTH = + KW > 0 ? KW : MAX_RUNTIME_KERNEL_WIDTH; + static constexpr int MAX_SPAN = EX + MAX_KERNEL_WIDTH - 1; + +public: + FK_STATIC_STRUCT(MorphQuadDPP, SelfType) + static constexpr ParArch PAR_ARCH = ParArch::GPU_NVIDIA; + static constexpr int BLOCK_THREADS = 256; + + struct LaunchConfig { + unsigned int blocks; + unsigned int threads; + }; + + FK_HOST_DEVICE_FUSE bool accepts(const MorphQuadDetails& details) { + return details.template validFor( + MAX_RUNTIME_KERNEL_WIDTH); + } + + FK_HOST_FUSE LaunchConfig launchConfig( + const MorphQuadDetails& details) { + const int columns = (details.width + EX - 1) / EX; + const int rows = (details.height + EY - 1) / EY; + const int workItems = columns * rows; + return {static_cast( + (workItems + BLOCK_THREADS - 1) / BLOCK_THREADS), + static_cast(BLOCK_THREADS)}; + } + + template + FK_DEVICE_FUSE void exec(const MorphQuadDetails& details, + const InIOp& input, + const ReduceIOp& reduce, + const OutIOp& output) { + static_assert(isAnyCompleteReadType, + "MorphQuadDPP requires a complete Read IOp"); + static_assert(isAnyWriteType, + "MorphQuadDPP requires a Write IOp"); + static_assert(std::is_same_v< + typename InIOp::Operation::OutputType, T>, + "MorphQuadDPP Read IOp must produce T"); + static_assert(std::is_same_v< + typename OutIOp::Operation::InputType, T>, + "MorphQuadDPP Write IOp must consume T"); + static_assert(EX > 0 && EY > 0, + "MorphQuadDPP output tile must be positive"); + static_assert(KW >= 0 && KH >= 0, + "MorphQuadDPP static kernel sizes cannot be negative"); + static_assert((KW == 0) == (KH == 0), + "MorphQuadDPP sizes are both static or both runtime"); + static_assert(cn == 1, + "MorphQuadDPP currently supports scalar pixel types"); +#if defined(__CUDA_ARCH__) + if (!accepts(details)) return; + const int kernelWidth = KW > 0 ? KW : details.kernelWidth; + const int kernelHeight = KH > 0 ? KH : details.kernelHeight; + const int columns = (details.width + EX - 1) / EX; + const int workItem = static_cast( + blockIdx.x * blockDim.x + threadIdx.x); + const int x0 = (workItem % columns) * EX; + const int y0 = (workItem / columns) * EY; + if (x0 >= details.width || y0 >= details.height) return; + + const int span = EX + kernelWidth - 1; + if (span > MAX_SPAN) return; + const int firstColumn = x0 - details.anchorX; + auto clamp = [](const int value, const int upper) { + return value < 0 ? 0 : (value >= upper ? upper - 1 : value); + }; + auto source = [&](const int x, const int y) { + return InIOp::Operation::exec( + Point{clamp(x, details.width), + clamp(y, details.height), 0}, input); + }; + + #pragma unroll + for (int localY = 0; localY < EY; ++localY) { + const int y = y0 + localY; + if (y >= details.height) break; + T columnsReduced[MAX_SPAN]; + #pragma unroll + for (int index = 0; index < span; ++index) { + T value = source(firstColumn + index, + y - details.anchorY); + for (int ky = 1; ky < kernelHeight; ++ky) { + value = make_tuple( + value, + source(firstColumn + index, + y + ky - details.anchorY)) | reduce; + } + columnsReduced[index] = value; + } + + #pragma unroll + for (int localX = 0; localX < EX; ++localX) { + const int x = x0 + localX; + if (x >= details.width) break; + T value = columnsReduced[localX]; + for (int kx = 1; kx < kernelWidth; ++kx) { + value = make_tuple( + value, columnsReduced[localX + kx]) | reduce; + } + OutIOp::Operation::exec(Point{x, y, 0}, value, output); + } + } +#endif // defined(__CUDA_ARCH__) + } +}; + +template +__global__ void launchMorphQuadDPP_Kernel( + const __grid_constant__ MorphQuadDetails details, + const __grid_constant__ IOps... iOps) { + DPP::exec(details, iOps...); +} + +template +FK_HOST_FUSE void executeMorphQuad( + Stream_& stream, + const MorphQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::GPU_NVIDIA, + "GPU stream requires the NVIDIA MorphQuadDPP specialization"); + if (!DPP::accepts(details)) return; + const auto launch = DPP::launchConfig(details); + launchMorphQuadDPP_Kernel + <<>>( + details, iOps...); + gpuErrchk(cudaGetLastError()); +} +#endif // defined(__NVCC__) + +template +FK_HOST_FUSE void executeMorphQuad( + Stream_&, + const MorphQuadDetails& details, + const IOps&... iOps) { + static_assert(DPP::PAR_ARCH == ParArch::CPU, + "CPU stream requires the CPU MorphQuadDPP specialization"); + DPP::exec(details, iOps...); +} + +} // namespace fk + +#endif // FK_MORPHOLOGY_FAST_H diff --git a/tests/image_processing/test_convolution_fast_dpp.h b/tests/image_processing/test_convolution_fast_dpp.h new file mode 100644 index 00000000..f975cbe7 --- /dev/null +++ b/tests/image_processing/test_convolution_fast_dpp.h @@ -0,0 +1,175 @@ +/* Copyright 2026 Johnny Nunez + + 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. */ + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace fk; + +namespace { +constexpr float READ_BIAS = 0.25f; +constexpr float WRITE_BIAS = -0.5f; + +float sourceValue(const int x, const int y) { + return static_cast(((x * 19 + y * 23) % 97) - 48) / 17.f; +} + +std::vector makeCoefficients(const int width, const int height) { + std::vector values(static_cast(width * height)); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + values[static_cast(y * width + x)] = + static_cast(((x + 1) * 7 + (y + 2) * 11) % 29 - 14) / + static_cast(width * height * 5); + } + } + return values; +} + +float oracleAt(const int width, const int height, + const int kernelWidth, const int kernelHeight, + const int anchorX, const int anchorY, + const std::vector& coefficients, + const int x, const int y) { + float result = 0.f; + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + const int sx = std::max(0, std::min( + width - 1, x + kx - anchorX)); + const int sy = std::max(0, std::min( + height - 1, y + ky - anchorY)); + result += coefficients[static_cast( + ky * kernelWidth + kx)] * + (sourceValue(sx, sy) + READ_BIAS); + } + } + return result + WRITE_BIAS; +} + +struct Result { + bool passed; + std::vector output; +}; + +template +Result runCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY) { + constexpr bool GPU = PA == ParArch::GPU_NVIDIA; + const auto memoryType = GPU ? MemType::DeviceAndPinned : MemType::Host; + Ptr2D input(width, height, 0, memoryType); + Ptr2D output(width, height, 0, memoryType); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + input.at(Point{x, y, 0}) = sourceValue(x, y); + output.at(Point{x, y, 0}) = 1234.f; + } + } + + const int kernelWidth = KW > 0 ? KW : runtimeKW; + const int kernelHeight = KH > 0 ? KH : runtimeKH; + const auto coefficients = makeCoefficients(kernelWidth, kernelHeight); + ConvQuadDetails details{ + width, height, runtimeKW, runtimeKH, anchorX, anchorY, {}}; + std::copy(coefficients.begin(), coefficients.end(), details.coefficients); + + Stream_ stream; +#if defined(__NVCC__) + if constexpr (GPU) { + input.upload(stream); + output.upload(stream); + } +#endif + const auto read = PerThreadRead::build(input) + .then(Add::build(READ_BIAS)); + const auto compute = make_tuple( + Mul::build(), + Add::build()); + const auto write = Add::build(WRITE_BIAS) + .then(PerThreadWrite::build(output)); + using DPP = ConvQuadDPP; + executeConvQuad(stream, details, read, compute, write); +#if defined(__NVCC__) + if constexpr (GPU) output.download(stream); +#endif + stream.sync(); + + bool passed = true; + std::vector actual(static_cast(width * height)); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + const float value = output.at(Point{x, y, 0}); + actual[static_cast(y * width + x)] = value; + passed = std::fabs(value - oracleAt( + width, height, kernelWidth, kernelHeight, + anchorX, anchorY, coefficients, x, y)) <= 2e-5f && passed; + } + } + return {passed, std::move(actual)}; +} + +template +bool verifyCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY, + const char* name) { + const auto cpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY); + bool ok = cpu.passed; +#if defined(__NVCC__) + const auto gpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY); + bool parity = gpu.output.size() == cpu.output.size(); + for (std::size_t i = 0; i < gpu.output.size() && parity; ++i) { + parity = std::fabs(gpu.output[i] - cpu.output[i]) <= 2e-5f; + } + ok = ok && gpu.passed && parity; + std::printf("ConvQuad %-18s %dx%d k%dx%d CPU/GPU %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#else + std::printf("ConvQuad CPU %-14s %dx%d k%dx%d %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#endif + return ok; +} +} // namespace + +int launch() { + bool ok = true; + ok = verifyCase<4, 4, 3, 3>( + 257, 129, 3, 3, 1, 1, "3x3-odd") && ok; + ok = verifyCase<4, 4, 5, 5>( + 192, 108, 5, 5, 2, 2, "5x5") && ok; + ok = verifyCase<2, 3, 7, 7>( + 65, 37, 7, 7, 3, 3, "7x7-ragged") && ok; + ok = verifyCase<4, 4, 7, 7>( + 39, 27, 1, 1, 3, 3, "static-details") && ok; + ok = verifyCase<4, 4, 0, 0>( + 73, 41, 3, 5, 0, 3, "runtime-3x5") && ok; + return ok ? 0 : -1; +} diff --git a/tests/image_processing/test_median_fast_dpp.h b/tests/image_processing/test_median_fast_dpp.h new file mode 100644 index 00000000..56e23208 --- /dev/null +++ b/tests/image_processing/test_median_fast_dpp.h @@ -0,0 +1,160 @@ +/* Copyright 2026 Johnny Nunez + + 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. */ + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace fk; + +namespace { +constexpr unsigned char READ_BIAS = 2; +constexpr unsigned char WRITE_BIAS = 1; + +unsigned char sourceValue(const int x, const int y) { + return static_cast((x * 37 + y * 13 + (x ^ (y * 3))) % 241); +} + +unsigned char oracleAt(const int width, const int height, + const int kernelWidth, const int kernelHeight, + const int anchorX, const int anchorY, + const int x, const int y) { + std::vector values; + values.reserve(static_cast(kernelWidth * kernelHeight)); + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + const int sx = std::max(0, std::min( + width - 1, x + kx - anchorX)); + const int sy = std::max(0, std::min( + height - 1, y + ky - anchorY)); + values.push_back(static_cast( + sourceValue(sx, sy) + READ_BIAS)); + } + } + std::sort(values.begin(), values.end()); + return static_cast( + values[values.size() / 2] + WRITE_BIAS); +} + +struct Result { + bool passed; + std::vector output; +}; + +template +Result runCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY) { + constexpr bool GPU = PA == ParArch::GPU_NVIDIA; + const auto memoryType = GPU ? MemType::DeviceAndPinned : MemType::Host; + Ptr2D input(width, height, 0, memoryType); + Ptr2D output(width, height, 0, memoryType); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + input.at(Point{x, y, 0}) = sourceValue(x, y); + output.at(Point{x, y, 0}) = 0xA5; + } + } + + Stream_ stream; +#if defined(__NVCC__) + if constexpr (GPU) { + input.upload(stream); + output.upload(stream); + } +#endif + const auto read = PerThreadRead::build(input) + .then(Add::build(READ_BIAS)); + const auto compare = make_tuple( + Min::build(), + Max::build()); + const auto write = Add::build(WRITE_BIAS) + .then(PerThreadWrite::build(output)); + using DPP = MedianQuadDPP; + const MedianQuadDetails details{ + width, height, runtimeKW, runtimeKH, anchorX, anchorY}; + executeMedianQuad(stream, details, read, compare, write); +#if defined(__NVCC__) + if constexpr (GPU) output.download(stream); +#endif + stream.sync(); + + const int kernelWidth = KW > 0 ? KW : runtimeKW; + const int kernelHeight = KH > 0 ? KH : runtimeKH; + bool passed = true; + std::vector actual( + static_cast(width * height)); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + const auto value = output.at(Point{x, y, 0}); + actual[static_cast(y * width + x)] = value; + passed = value == oracleAt( + width, height, kernelWidth, kernelHeight, + anchorX, anchorY, x, y) && passed; + } + } + return {passed, std::move(actual)}; +} + +template +bool verifyCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY, + const char* name) { + const auto cpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY); + bool ok = cpu.passed; +#if defined(__NVCC__) + const auto gpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY); + ok = ok && gpu.passed && gpu.output == cpu.output; + std::printf("MedianQuad %-16s %dx%d k%dx%d CPU/GPU %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#else + std::printf("MedianQuad CPU %-12s %dx%d k%dx%d %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#endif + return ok; +} +} // namespace + +int launch() { + bool ok = true; + ok = verifyCase<4, 4, 3, 3>( + 257, 129, 3, 3, 1, 1, "3x3-odd") && ok; + ok = verifyCase<4, 4, 5, 5>( + 192, 108, 5, 5, 2, 2, "5x5") && ok; + ok = verifyCase<2, 3, 7, 7>( + 65, 37, 7, 7, 3, 3, "7x7-ragged") && ok; + ok = verifyCase<4, 4, 7, 7>( + 39, 27, 1, 1, 3, 3, "static-details") && ok; + ok = verifyCase<4, 4, 0, 0>( + 73, 41, 3, 5, 0, 3, "runtime-3x5") && ok; + return ok ? 0 : -1; +} diff --git a/tests/image_processing/test_morphology_fast_dpp.h b/tests/image_processing/test_morphology_fast_dpp.h new file mode 100644 index 00000000..335fd38b --- /dev/null +++ b/tests/image_processing/test_morphology_fast_dpp.h @@ -0,0 +1,170 @@ +/* Copyright 2026 Johnny Nunez + + 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. */ + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace fk; + +namespace { +constexpr unsigned char READ_BIAS = 3; +constexpr unsigned char WRITE_BIAS = 1; + +enum class MorphologyKind { Erode, Dilate }; + +unsigned char sourceValue(const int x, const int y) { + return static_cast((x * 31 + y * 17 + (x ^ y)) % 241); +} + +unsigned char oracleAt(const int width, const int height, + const int kernelWidth, const int kernelHeight, + const int anchorX, const int anchorY, + const int x, const int y, + const MorphologyKind kind) { + unsigned char result = kind == MorphologyKind::Erode ? 255 : 0; + for (int ky = 0; ky < kernelHeight; ++ky) { + for (int kx = 0; kx < kernelWidth; ++kx) { + const int sx = std::max(0, std::min( + width - 1, x + kx - anchorX)); + const int sy = std::max(0, std::min( + height - 1, y + ky - anchorY)); + const auto value = static_cast( + sourceValue(sx, sy) + READ_BIAS); + result = kind == MorphologyKind::Erode + ? std::min(result, value) : std::max(result, value); + } + } + return static_cast(result + WRITE_BIAS); +} + +struct Result { + bool passed; + std::vector output; +}; + +template +Result runCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY, + const MorphologyKind kind) { + constexpr bool GPU = PA == ParArch::GPU_NVIDIA; + const auto memoryType = GPU ? MemType::DeviceAndPinned : MemType::Host; + Ptr2D input(width, height, 0, memoryType); + Ptr2D output(width, height, 0, memoryType); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + input.at(Point{x, y, 0}) = sourceValue(x, y); + output.at(Point{x, y, 0}) = 0xA5; + } + } + + Stream_ stream; +#if defined(__NVCC__) + if constexpr (GPU) { + input.upload(stream); + output.upload(stream); + } +#endif + const auto read = PerThreadRead::build(input) + .then(Add::build(READ_BIAS)); + const auto reduce = ReduceIOp::build(); + const auto write = Add::build(WRITE_BIAS) + .then(PerThreadWrite::build(output)); + using DPP = MorphQuadDPP; + const MorphQuadDetails details{ + width, height, runtimeKW, runtimeKH, anchorX, anchorY}; + executeMorphQuad(stream, details, read, reduce, write); +#if defined(__NVCC__) + if constexpr (GPU) output.download(stream); +#endif + stream.sync(); + + const int kernelWidth = KW > 0 ? KW : runtimeKW; + const int kernelHeight = KH > 0 ? KH : runtimeKH; + bool passed = true; + std::vector actual( + static_cast(width * height)); + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + const auto value = output.at(Point{x, y, 0}); + actual[static_cast(y * width + x)] = value; + passed = value == oracleAt( + width, height, kernelWidth, kernelHeight, + anchorX, anchorY, x, y, kind) && passed; + } + } + return {passed, std::move(actual)}; +} + +template +bool verifyCase(const int width, const int height, + const int runtimeKW, const int runtimeKH, + const int anchorX, const int anchorY, + const MorphologyKind kind, const char* name) { + const auto cpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY, kind); + bool ok = cpu.passed; +#if defined(__NVCC__) + const auto gpu = runCase( + width, height, runtimeKW, runtimeKH, anchorX, anchorY, kind); + ok = ok && gpu.passed && gpu.output == cpu.output; + std::printf("MorphQuad %-20s %dx%d k%dx%d CPU/GPU %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#else + std::printf("MorphQuad CPU %-16s %dx%d k%dx%d %s\n", + name, width, height, + KW > 0 ? KW : runtimeKW, KH > 0 ? KH : runtimeKH, + ok ? "PASS" : "FAIL"); +#endif + return ok; +} +} // namespace + +int launch() { + using MinIOp = Min; + using MaxIOp = Max; + bool ok = true; + ok = verifyCase( + 257, 129, 3, 3, 1, 1, MorphologyKind::Erode, + "erode-3x3-odd") && ok; + ok = verifyCase( + 192, 108, 5, 5, 2, 2, MorphologyKind::Dilate, + "dilate-5x5") && ok; + ok = verifyCase( + 65, 37, 7, 7, 3, 3, MorphologyKind::Erode, + "erode-7x7-ragged") && ok; + ok = verifyCase( + 39, 27, 1, 1, 3, 3, MorphologyKind::Dilate, + "static-details") && ok; + ok = verifyCase( + 73, 41, 5, 3, 1, 0, MorphologyKind::Erode, + "runtime-5x3") && ok; + return ok ? 0 : -1; +}