diff --git a/modules/photo/include/opencv2/photo.hpp b/modules/photo/include/opencv2/photo.hpp index a8b7c6049a..a4a6e6214a 100644 --- a/modules/photo/include/opencv2/photo.hpp +++ b/modules/photo/include/opencv2/photo.hpp @@ -612,13 +612,29 @@ public: /** @brief The resulting HDR image is calculated as weighted average of the exposures considering exposure values and camera response. +Input images must be 1-channel or 3-channel (CV_8U, CV_16U, or CV_32F). + For more information see @cite DM97 . */ class CV_EXPORTS_W MergeDebevec : public MergeExposures { public: + /** @brief Merges images. + + @param src vector of input images, all 1-channel or all 3-channel, CV_8U, CV_16U, or CV_32F + @param dst result image + @param times vector of exposure time values for each image + @param response 256x1 or 65536x1 matrix with inverse camera response function for each pixel value; it must + have the same number of channels as images (empty means a linear response) + */ CV_WRAP virtual void process(InputArrayOfArrays src, OutputArray dst, InputArray times, InputArray response) CV_OVERRIDE = 0; + /** @brief Short version of process that uses a linear camera response. + + @param src vector of input images, all 1-channel or all 3-channel, CV_8U, CV_16U, or CV_32F + @param dst result image + @param times vector of exposure time values for each image + */ CV_WRAP virtual void process(InputArrayOfArrays src, OutputArray dst, InputArray times) = 0; }; diff --git a/modules/photo/src/merge.cpp b/modules/photo/src/merge.cpp index 2f7868b4bd..cbb28be886 100644 --- a/modules/photo/src/merge.cpp +++ b/modules/photo/src/merge.cpp @@ -11,6 +11,7 @@ // For Open Source Computer Vision Library // // Copyright (C) 2013, OpenCV Foundation, all rights reserved. +// Copyright (C) 2026, Advanced Micro Devices, Inc., all rights reserved. // Third party copyrights are property of their respective owners. // // Redistribution and use in source and binary forms, with or without modification, @@ -42,10 +43,277 @@ #include "precomp.hpp" #include "opencv2/photo.hpp" #include "opencv2/imgproc.hpp" +#include "opencv2/core/hal/intrin.hpp" #include "hdr_common.hpp" namespace cv { +namespace +{ + +static void scaleHatWeights(const Mat& weights_lut, int channels, std::vector& hat) +{ + const int n = weights_lut.rows; + hat.resize((size_t)n); + const float inv_cn = 1.0f / static_cast(channels); + for (int z = 0; z < n; z++) + hat[z] = weights_lut.at(z) * inv_cn; +} + +static void extractLogResponseC1(const Mat& log_response, std::vector& g0) +{ + const int n = log_response.rows; + g0.resize((size_t)n); + for (int z = 0; z < n; z++) + g0[z] = log_response.at(z); +} + +static void extractLogResponseC3(const Mat& log_response, + std::vector& g0, std::vector& g1, std::vector& g2) +{ + const int n = log_response.rows; + g0.resize((size_t)n); + g1.resize((size_t)n); + g2.resize((size_t)n); + for (int z = 0; z < n; z++) + { + const Vec3f v = log_response.at(z); + g0[z] = v[0]; + g1[z] = v[1]; + g2[z] = v[2]; + } +} + +#if CV_SIMD || CV_SIMD_SCALABLE +static inline v_float32 debevecTriangleHat(const v_int32& idx, float maxVal, float inv_cn) +{ + const v_float32 z = v_cvt_f32(idx); + v_float32 h = v_min(z, v_sub(vx_setall_f32(maxVal), z)); + h = v_max(h, vx_setall_f32(1e-6f)); + return v_mul(h, vx_setall_f32(inv_cn)); +} + +static inline void debevecAccQuad(const v_uint32& ib, const v_uint32& ig, const v_uint32& ir, + const float* g0, const float* g1, const float* g2, + const v_float32& ln, float maxVal, float inv_cn, + float* p0, float* p1, float* p2, float* pw) +{ + const v_int32 sb = v_reinterpret_as_s32(ib); + const v_int32 sg = v_reinterpret_as_s32(ig); + const v_int32 sr = v_reinterpret_as_s32(ir); + const v_float32 w = v_add(v_add(debevecTriangleHat(sb, maxVal, inv_cn), + debevecTriangleHat(sg, maxVal, inv_cn)), + debevecTriangleHat(sr, maxVal, inv_cn)); + v_store(p0, v_fma(w, v_sub(v_lut(g0, sb), ln), vx_load(p0))); + v_store(p1, v_fma(w, v_sub(v_lut(g1, sg), ln), vx_load(p1))); + v_store(p2, v_fma(w, v_sub(v_lut(g2, sr), ln), vx_load(p2))); + v_store(pw, v_add(vx_load(pw), w)); +} + +static inline void debevecAccQuadC1(const v_uint32& iz, const float* g0, + const v_float32& ln, float maxVal, float inv_cn, + float* p0, float* pw) +{ + const v_int32 sz = v_reinterpret_as_s32(iz); + const v_float32 w = debevecTriangleHat(sz, maxVal, inv_cn); + v_store(p0, v_fma(w, v_sub(v_lut(g0, sz), ln), vx_load(p0))); + v_store(pw, v_add(vx_load(pw), w)); +} +#endif + +template +static void debevecAccumulateC1(const Mat& img, const float* hat, const float* g0, float ln_t, + Mat& s0, Mat& wsum) +{ + const int cols = img.cols; +#if CV_SIMD || CV_SIMD_SCALABLE + const float maxVal = (sizeof(IdxT) == 1) ? 255.f : 65535.f; + const float inv_cn = 1.0f; +#endif + parallel_for_(Range(0, img.rows), [&](const Range& range) + { + for (int y = range.start; y < range.end; y++) + { + const IdxT* src = img.ptr(y); + float* p0 = s0.ptr(y); + float* pw = wsum.ptr(y); + int x = 0; +#if CV_SIMD || CV_SIMD_SCALABLE + const v_float32 ln = vx_setall_f32(ln_t); + const int n32 = VTraits::vlanes(); + if (sizeof(IdxT) == 1) + { + for (; x <= cols - n32; x += n32) + { + const v_uint32 q = vx_load_expand_q(reinterpret_cast(src + x)); + debevecAccQuadC1(q, g0, ln, maxVal, inv_cn, p0 + x, pw + x); + } + } + else + { + for (; x <= cols - n32; x += n32) + { + const v_uint32 q = vx_load_expand(reinterpret_cast(src + x)); + debevecAccQuadC1(q, g0, ln, maxVal, inv_cn, p0 + x, pw + x); + } + } +#endif + for (; x < cols; x++) + { + const IdxT z = src[x]; + const float w = hat[z]; + p0[x] += w * (g0[z] - ln_t); + pw[x] += w; + } + } + }); +} + +template +static void debevecAccumulateC3(const Mat& img, const float* hat, + const float* g0, const float* g1, const float* g2, + float ln_t, Mat& s0, Mat& s1, Mat& s2, Mat& wsum) +{ + const int cols = img.cols; +#if CV_SIMD || CV_SIMD_SCALABLE + const float maxVal = (sizeof(IdxT) == 1) ? 255.f : 65535.f; + const float inv_cn = 1.0f / 3.0f; +#endif + parallel_for_(Range(0, img.rows), [&](const Range& range) + { + for (int y = range.start; y < range.end; y++) + { + const IdxT* src = img.ptr(y); + float* p0 = s0.ptr(y); + float* p1 = s1.ptr(y); + float* p2 = s2.ptr(y); + float* pw = wsum.ptr(y); + int x = 0; +#if CV_SIMD || CV_SIMD_SCALABLE + const v_float32 ln = vx_setall_f32(ln_t); + if (sizeof(IdxT) == 1) + { + const int n8 = VTraits::vlanes(); + for (; x <= cols - n8; x += n8, src += n8 * 3) + { + v_uint8 vb, vg, vr; + v_load_deinterleave(reinterpret_cast(src), vb, vg, vr); + v_uint16 b0, b1, g0u, g1u, r0u, r1u; + v_expand(vb, b0, b1); + v_expand(vg, g0u, g1u); + v_expand(vr, r0u, r1u); + v_uint32 qb0, qb1, qb2, qb3, qg0, qg1, qg2, qg3, qr0, qr1, qr2, qr3; + v_expand(b0, qb0, qb1); + v_expand(b1, qb2, qb3); + v_expand(g0u, qg0, qg1); + v_expand(g1u, qg2, qg3); + v_expand(r0u, qr0, qr1); + v_expand(r1u, qr2, qr3); + const int n32 = VTraits::vlanes(); + debevecAccQuad(qb0, qg0, qr0, g0, g1, g2, ln, maxVal, inv_cn, p0 + x, p1 + x, p2 + x, pw + x); + debevecAccQuad(qb1, qg1, qr1, g0, g1, g2, ln, maxVal, inv_cn, + p0 + x + n32, p1 + x + n32, p2 + x + n32, pw + x + n32); + debevecAccQuad(qb2, qg2, qr2, g0, g1, g2, ln, maxVal, inv_cn, + p0 + x + 2 * n32, p1 + x + 2 * n32, p2 + x + 2 * n32, pw + x + 2 * n32); + debevecAccQuad(qb3, qg3, qr3, g0, g1, g2, ln, maxVal, inv_cn, + p0 + x + 3 * n32, p1 + x + 3 * n32, p2 + x + 3 * n32, pw + x + 3 * n32); + } + } + else + { + const int n16 = VTraits::vlanes(); + for (; x <= cols - n16; x += n16, src += n16 * 3) + { + v_uint16 vb, vg, vr; + v_load_deinterleave(reinterpret_cast(src), vb, vg, vr); + v_uint32 qb0, qb1, qg0, qg1, qr0, qr1; + v_expand(vb, qb0, qb1); + v_expand(vg, qg0, qg1); + v_expand(vr, qr0, qr1); + const int n32 = VTraits::vlanes(); + debevecAccQuad(qb0, qg0, qr0, g0, g1, g2, ln, maxVal, inv_cn, p0 + x, p1 + x, p2 + x, pw + x); + debevecAccQuad(qb1, qg1, qr1, g0, g1, g2, ln, maxVal, inv_cn, + p0 + x + n32, p1 + x + n32, p2 + x + n32, pw + x + n32); + } + } + src = img.ptr(y) + x * 3; +#endif + for (; x < cols; x++, src += 3) + { + const IdxT b = src[0], g = src[1], r = src[2]; + const float w = hat[b] + hat[g] + hat[r]; + p0[x] += w * (g0[b] - ln_t); + p1[x] += w * (g1[g] - ln_t); + p2[x] += w * (g2[r] - ln_t); + pw[x] += w; + } + } + }); +} + +static void debevecFinalizeC1(const Mat& s0, const Mat& wsum, Mat& dst) +{ + const int cols = dst.cols; + parallel_for_(Range(0, dst.rows), [&](const Range& range) + { + for (int y = range.start; y < range.end; y++) + { + const float* p0 = s0.ptr(y); + const float* pw = wsum.ptr(y); + float* d = dst.ptr(y); + int x = 0; +#if CV_SIMD || CV_SIMD_SCALABLE + const int n32 = VTraits::vlanes(); + const v_float32 one = vx_setall_f32(1.0f); + for (; x <= cols - n32; x += n32) + { + const v_float32 inv = v_div(one, vx_load(pw + x)); + v_store(d + x, v_exp(v_mul(vx_load(p0 + x), inv))); + } +#endif + for (; x < cols; x++) + d[x] = std::exp(p0[x] * (1.0f / pw[x])); + } + }); +} + +static void debevecFinalizeC3(const Mat& s0, const Mat& s1, const Mat& s2, const Mat& wsum, Mat& dst) +{ + const int cols = dst.cols; + parallel_for_(Range(0, dst.rows), [&](const Range& range) + { + for (int y = range.start; y < range.end; y++) + { + const float* p0 = s0.ptr(y); + const float* p1 = s1.ptr(y); + const float* p2 = s2.ptr(y); + const float* pw = wsum.ptr(y); + float* d = dst.ptr(y); + int x = 0; +#if CV_SIMD || CV_SIMD_SCALABLE + const int n32 = VTraits::vlanes(); + const v_float32 one = vx_setall_f32(1.0f); + for (; x <= cols - n32; x += n32) + { + const v_float32 inv = v_div(one, vx_load(pw + x)); + v_store_interleave(d + x * 3, + v_exp(v_mul(vx_load(p0 + x), inv)), + v_exp(v_mul(vx_load(p1 + x), inv)), + v_exp(v_mul(vx_load(p2 + x), inv))); + } +#endif + for (; x < cols; x++) + { + const float inv = 1.0f / pw[x]; + d[x * 3 + 0] = std::exp(p0[x] * inv); + d[x * 3 + 1] = std::exp(p1[x] * inv); + d[x * 3 + 2] = std::exp(p2[x] * inv); + } + } + }); +} + +} // namespace class MergeDebevecImpl CV_FINAL : public MergeDebevec { @@ -70,6 +338,7 @@ public: CV_Assert(depth == CV_8U || depth == CV_16U || depth == CV_32F); int channels = images[0].channels(); + CV_Assert(channels == 1 || channels == 3); Size size = images[0].size(); int CV_32FCC = CV_MAKETYPE(CV_32F, channels); @@ -100,7 +369,10 @@ public: if(response.empty()) { response = linearResponse(channels, lutLength); - response.at(0) = response.at(1); + if (channels == 3) + response.at(0) = response.at(1); + else if (channels == 1) + response.at(0) = response.at(1); } Mat log_response; @@ -111,40 +383,51 @@ public: Mat exp_values(times.clone()); log(exp_values, exp_values); - result = Mat::zeros(size, CV_32FCC); - std::vector result_split; - split(result, result_split); - Mat weight_sum = Mat::zeros(size, CV_32F); - Mat weights_lut = use16bitLUT ? triangleWeights(lutLength) : weights; - for(size_t i = 0; i < images.size(); i++) { - std::vector splitted; - split(lutImages[i], splitted); + const int idxDepth = lutImages[0].depth(); + CV_Assert(idxDepth == CV_8U || idxDepth == CV_16U); - Mat w = Mat::zeros(size, CV_32F); - for(int c = 0; c < channels; c++) { - LUT(splitted[c], weights_lut, splitted[c]); - w += splitted[c]; - } - w /= channels; + std::vector hat; + scaleHatWeights(weights_lut, channels, hat); - Mat response_img; - LUT(lutImages[i], log_response, response_img); - split(response_img, splitted); - for(int c = 0; c < channels; c++) { - result_split[c] += w.mul(splitted[c] - exp_values.at((int)i)); + Mat weight_sum = Mat::zeros(size, CV_32F); + if (channels == 3) + { + std::vector g0, g1, g2; + extractLogResponseC3(log_response, g0, g1, g2); + Mat s0 = Mat::zeros(size, CV_32F); + Mat s1 = Mat::zeros(size, CV_32F); + Mat s2 = Mat::zeros(size, CV_32F); + for (size_t i = 0; i < lutImages.size(); i++) + { + const float ln_t = exp_values.at((int)i); + if (idxDepth == CV_8U) + debevecAccumulateC3(lutImages[i], hat.data(), g0.data(), g1.data(), g2.data(), + ln_t, s0, s1, s2, weight_sum); + else + debevecAccumulateC3(lutImages[i], hat.data(), g0.data(), g1.data(), g2.data(), + ln_t, s0, s1, s2, weight_sum); } - weight_sum += w; + debevecFinalizeC3(s0, s1, s2, weight_sum, result); } - weight_sum = 1.0f / weight_sum; - for(int c = 0; c < channels; c++) { - result_split[c] = result_split[c].mul(weight_sum); + else + { + std::vector g0; + extractLogResponseC1(log_response, g0); + Mat s0 = Mat::zeros(size, CV_32F); + for (size_t i = 0; i < lutImages.size(); i++) + { + const float ln_t = exp_values.at((int)i); + if (idxDepth == CV_8U) + debevecAccumulateC1(lutImages[i], hat.data(), g0.data(), ln_t, s0, weight_sum); + else + debevecAccumulateC1(lutImages[i], hat.data(), g0.data(), ln_t, s0, weight_sum); + } + debevecFinalizeC1(s0, weight_sum, result); } - merge(result_split, result); - exp(result, result); } void process(InputArrayOfArrays src, OutputArray dst, InputArray times) CV_OVERRIDE diff --git a/modules/photo/test/test_hdr.cpp b/modules/photo/test/test_hdr.cpp index 9922a01fff..cd4dc8373f 100644 --- a/modules/photo/test/test_hdr.cpp +++ b/modules/photo/test/test_hdr.cpp @@ -240,6 +240,102 @@ TEST(Photo_MergeDebevec, regression_depth_consistency) checkEqual(hdr8, hdr32, 2e-2f, "Debevec realdata 32F vs 8U"); } +static void toGraySeq(const vector& bgr, vector& gray) +{ + gray.resize(bgr.size()); + for (size_t i = 0; i < bgr.size(); ++i) + cvtColor(bgr[i], gray[i], COLOR_BGR2GRAY); +} + +static Mat linearResponseC1(int length) +{ + Mat response(length, 1, CV_32FC1); + for (int i = 0; i < length; i++) + response.at(i) = static_cast(i); + response.at(0) = response.at(1); + return response; +} + +TEST(Photo_MergeDebevec, regression_1channel) +{ + string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/"; + + vector images_bgr; + vector times; + loadExposureSeq(test_path + "exposures/", images_bgr, times); + + vector images; + toGraySeq(images_bgr, images); + + Ptr merge = createMergeDebevec(); + + Mat result_default; + merge->process(images, result_default, times); + ASSERT_EQ(result_default.type(), CV_32FC1); + ASSERT_EQ(result_default.size(), images[0].size()); + ASSERT_TRUE(cv::checkRange(result_default)) << "Debevec 1ch default produced non-finite values"; + + Mat result_linear; + merge->process(images, result_linear, times, linearResponseC1(256)); + checkEqual(result_default, result_linear, 1e-5f, "Debevec 1ch default vs linear CRF"); + + Mat response3; + loadResponseCSV(test_path + "exposures/response.csv", response3); + Mat response1(256, 1, CV_32FC1); + for (int i = 0; i < 256; i++) + response1.at(i) = response3.at(i)[0]; + Mat result_custom; + merge->process(images, result_custom, times, response1); + ASSERT_EQ(result_custom.type(), CV_32FC1); + ASSERT_TRUE(cv::checkRange(result_custom)) << "Debevec 1ch custom CRF produced non-finite values"; +} + +TEST(Photo_MergeDebevec, regression_1channel_depth_consistency) +{ + string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/"; + + vector images_bgr; + vector times; + loadExposureSeq(test_path + "exposures/", images_bgr, times); + + vector images8; + toGraySeq(images_bgr, images8); + + vector images16(images8.size()), images32(images8.size()); + for (size_t i = 0; i < images8.size(); ++i) + { + images8[i].convertTo(images16[i], CV_16UC1, 257.0); + images8[i].convertTo(images32[i], CV_32FC1, 1.0 / 255.0); + } + + Ptr debevec = createMergeDebevec(); + Ptr map = createTonemap(); + + Mat hdr8, hdr16, hdr32; + debevec->process(images8, hdr8, times); + debevec->process(images16, hdr16, times); + debevec->process(images32, hdr32, times); + ASSERT_TRUE(cv::checkRange(hdr8) && cv::checkRange(hdr16) && cv::checkRange(hdr32)); + + auto to3 = [](const Mat& c1) { + Mat c3; + std::vector ch{c1, c1, c1}; + cv::merge(ch, c3); + return c3; + }; + Mat t8, t16, t32; + map->process(to3(hdr8), t8); + map->process(to3(hdr16), t16); + map->process(to3(hdr32), t32); + + checkEqual(t8, t16, 2e-2f, "Debevec 1ch realdata 16U vs 8U"); + checkEqual(t8, t32, 2e-2f, "Debevec 1ch realdata 32F vs 8U"); + + Mat hdr16_linear; + debevec->process(images16, hdr16_linear, times, linearResponseC1(65536)); + checkEqual(hdr16, hdr16_linear, 1e-5f, "Debevec 1ch 16U default vs linear CRF"); +} + TEST(Photo_MergeRobertson, regression_depth_consistency) { string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/";