Fuse and vectorize MergeDebevec for 1- and 3-channel HDR.

Fold LUT/split/mul into one pass and SIMD the accumulate/finalize loops to cut full-image temporaries.
This commit is contained in:
Madan mohan Manokar
2026-09-17 09:21:54 +00:00
parent a31b3d5e85
commit 46b9862506
3 changed files with 421 additions and 26 deletions
+16
View File
@@ -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;
};
+309 -26
View File
@@ -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<float>& hat)
{
const int n = weights_lut.rows;
hat.resize((size_t)n);
const float inv_cn = 1.0f / static_cast<float>(channels);
for (int z = 0; z < n; z++)
hat[z] = weights_lut.at<float>(z) * inv_cn;
}
static void extractLogResponseC1(const Mat& log_response, std::vector<float>& g0)
{
const int n = log_response.rows;
g0.resize((size_t)n);
for (int z = 0; z < n; z++)
g0[z] = log_response.at<float>(z);
}
static void extractLogResponseC3(const Mat& log_response,
std::vector<float>& g0, std::vector<float>& g1, std::vector<float>& 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<Vec3f>(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<typename IdxT>
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<IdxT>(y);
float* p0 = s0.ptr<float>(y);
float* pw = wsum.ptr<float>(y);
int x = 0;
#if CV_SIMD || CV_SIMD_SCALABLE
const v_float32 ln = vx_setall_f32(ln_t);
const int n32 = VTraits<v_float32>::vlanes();
if (sizeof(IdxT) == 1)
{
for (; x <= cols - n32; x += n32)
{
const v_uint32 q = vx_load_expand_q(reinterpret_cast<const uchar*>(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<const ushort*>(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<typename IdxT>
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<IdxT>(y);
float* p0 = s0.ptr<float>(y);
float* p1 = s1.ptr<float>(y);
float* p2 = s2.ptr<float>(y);
float* pw = wsum.ptr<float>(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<v_uint8>::vlanes();
for (; x <= cols - n8; x += n8, src += n8 * 3)
{
v_uint8 vb, vg, vr;
v_load_deinterleave(reinterpret_cast<const uchar*>(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<v_float32>::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<v_uint16>::vlanes();
for (; x <= cols - n16; x += n16, src += n16 * 3)
{
v_uint16 vb, vg, vr;
v_load_deinterleave(reinterpret_cast<const ushort*>(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<v_float32>::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<IdxT>(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<float>(y);
const float* pw = wsum.ptr<float>(y);
float* d = dst.ptr<float>(y);
int x = 0;
#if CV_SIMD || CV_SIMD_SCALABLE
const int n32 = VTraits<v_float32>::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<float>(y);
const float* p1 = s1.ptr<float>(y);
const float* p2 = s2.ptr<float>(y);
const float* pw = wsum.ptr<float>(y);
float* d = dst.ptr<float>(y);
int x = 0;
#if CV_SIMD || CV_SIMD_SCALABLE
const int n32 = VTraits<v_float32>::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<Vec3f>(0) = response.at<Vec3f>(1);
if (channels == 3)
response.at<Vec3f>(0) = response.at<Vec3f>(1);
else if (channels == 1)
response.at<float>(0) = response.at<float>(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<Mat> 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<Mat> 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<float> 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<float>((int)i));
Mat weight_sum = Mat::zeros(size, CV_32F);
if (channels == 3)
{
std::vector<float> 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<float>((int)i);
if (idxDepth == CV_8U)
debevecAccumulateC3<uchar>(lutImages[i], hat.data(), g0.data(), g1.data(), g2.data(),
ln_t, s0, s1, s2, weight_sum);
else
debevecAccumulateC3<ushort>(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<float> 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<float>((int)i);
if (idxDepth == CV_8U)
debevecAccumulateC1<uchar>(lutImages[i], hat.data(), g0.data(), ln_t, s0, weight_sum);
else
debevecAccumulateC1<ushort>(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
+96
View File
@@ -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<Mat>& bgr, vector<Mat>& 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<float>(i) = static_cast<float>(i);
response.at<float>(0) = response.at<float>(1);
return response;
}
TEST(Photo_MergeDebevec, regression_1channel)
{
string test_path = string(cvtest::TS::ptr()->get_data_path()) + "hdr/";
vector<Mat> images_bgr;
vector<float> times;
loadExposureSeq(test_path + "exposures/", images_bgr, times);
vector<Mat> images;
toGraySeq(images_bgr, images);
Ptr<MergeDebevec> 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<float>(i) = response3.at<Vec3f>(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<Mat> images_bgr;
vector<float> times;
loadExposureSeq(test_path + "exposures/", images_bgr, times);
vector<Mat> images8;
toGraySeq(images_bgr, images8);
vector<Mat> 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<MergeDebevec> debevec = createMergeDebevec();
Ptr<Tonemap> 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<Mat> 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/";