Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "blend" step of "raster pipeline"
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : * Copyright (c) 2009, Frank Warmerdam
10 :
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdalalg_raster_blend.h"
15 :
16 : #include "cpl_conv.h"
17 : #include "gdal_priv.h"
18 :
19 : #include <algorithm>
20 : #include <array>
21 : #include <cassert>
22 : #include <limits>
23 : #include <utility>
24 :
25 : #if defined(__x86_64) || defined(_M_X64)
26 : #define HAVE_SSE2
27 : #include <immintrin.h>
28 : #endif
29 : #ifdef HAVE_SSE2
30 : #include "gdalsse_priv.h"
31 : #endif
32 :
33 : //! @cond Doxygen_Suppress
34 :
35 : #ifndef _
36 : #define _(x) (x)
37 : #endif
38 :
39 : constexpr const char *SRC_OVER = "src-over";
40 : constexpr const char *HSV_VALUE = "hsv-value";
41 :
42 : /************************************************************************/
43 : /* GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm() */
44 : /************************************************************************/
45 :
46 88 : GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm(bool standaloneStep)
47 : : GDALRasterPipelineStepAlgorithm(
48 : NAME, DESCRIPTION, HELP_URL,
49 0 : ConstructorOptions()
50 88 : .SetStandaloneStep(standaloneStep)
51 88 : .SetAddDefaultArguments(false)
52 176 : .SetInputDatasetHelpMsg(_("Input raster dataset"))
53 176 : .SetInputDatasetAlias("color-input")
54 176 : .SetInputDatasetMetaVar("COLOR-INPUT")
55 264 : .SetOutputDatasetHelpMsg(_("Output raster dataset")))
56 : {
57 88 : const auto AddOverlayDatasetArg = [this]()
58 : {
59 : auto &arg = AddArg("overlay", 0, _("Overlay dataset"),
60 176 : &m_overlayDataset, GDAL_OF_RASTER)
61 88 : .SetPositional()
62 88 : .SetRequired();
63 :
64 88 : SetAutoCompleteFunctionForFilename(arg, GDAL_OF_RASTER);
65 88 : };
66 :
67 88 : if (standaloneStep)
68 : {
69 46 : AddRasterInputArgs(false, false);
70 46 : AddOverlayDatasetArg();
71 46 : AddProgressArg();
72 46 : AddRasterOutputArgs(false);
73 : }
74 : else
75 : {
76 42 : AddRasterHiddenInputDatasetArg();
77 42 : AddOverlayDatasetArg();
78 : }
79 :
80 176 : AddArg("operator", 0, _("Composition operator"), &m_operator)
81 88 : .SetChoices(SRC_OVER, HSV_VALUE)
82 88 : .SetDefault(SRC_OVER);
83 : AddArg("opacity", 0,
84 : _("Opacity percentage to apply to the overlay dataset (0=fully "
85 : "transparent, 100=full use of overlay opacity)"),
86 176 : &m_opacity)
87 88 : .SetDefault(m_opacity)
88 88 : .SetMinValueIncluded(0)
89 88 : .SetMaxValueIncluded(OPACITY_INPUT_RANGE);
90 :
91 132 : AddValidationAction([this]() { return ValidateGlobal(); });
92 88 : }
93 :
94 : namespace
95 : {
96 :
97 : /************************************************************************/
98 : /* BlendDataset */
99 : /************************************************************************/
100 :
101 : class BlendDataset final : public GDALDataset
102 : {
103 : public:
104 : BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,
105 : const std::string &sOperator, int nOpacity255Scale);
106 : ~BlendDataset() override;
107 :
108 2 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
109 : {
110 2 : return m_oColorDS.GetGeoTransform(gt);
111 : }
112 :
113 2 : const OGRSpatialReference *GetSpatialRef() const override
114 : {
115 2 : return m_oColorDS.GetSpatialRef();
116 : }
117 :
118 : bool AcquireSourcePixels(int nXOff, int nYOff, int nXSize, int nYSize,
119 : int nBufXSize, int nBufYSize,
120 : GDALRasterIOExtraArg *psExtraArg);
121 :
122 : protected:
123 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
124 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
125 : GDALDataType eBufType, int nBandCount,
126 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
127 : GSpacing nLineSpace, GSpacing nBandSpace,
128 : GDALRasterIOExtraArg *psExtraArg) override;
129 :
130 : private:
131 : friend class BlendBand;
132 : GDALDataset &m_oColorDS;
133 : GDALDataset &m_oOverlayDS;
134 : const std::string m_operator;
135 : const int m_opacity255Scale;
136 : std::vector<std::unique_ptr<BlendDataset>> m_apoOverviews{};
137 : int m_nCachedXOff = 0;
138 : int m_nCachedYOff = 0;
139 : int m_nCachedXSize = 0;
140 : int m_nCachedYSize = 0;
141 : int m_nCachedBufXSize = 0;
142 : int m_nCachedBufYSize = 0;
143 : GDALRasterIOExtraArg m_sCachedExtraArg{};
144 : std::vector<GByte> m_abyBuffer{};
145 : bool m_ioError = false;
146 : };
147 :
148 : /************************************************************************/
149 : /* rgb_to_hs() */
150 : /************************************************************************/
151 :
152 : // rgb comes in as [r,g,b] with values in the range [0,255]. The returned
153 : // values will be with hue and saturation in the range [0,1].
154 :
155 : // Derived from hsv_merge.py
156 :
157 4982560 : static void rgb_to_hs(int r, int g, int b, float *h, float *s)
158 : {
159 : int minc, maxc;
160 4982560 : if (r <= g)
161 : {
162 2531110 : if (r <= b)
163 : {
164 1700760 : minc = r;
165 1700760 : maxc = std::max(g, b);
166 : }
167 : else /* b < r */
168 : {
169 830353 : minc = b;
170 830353 : maxc = g;
171 : }
172 : }
173 : else /* g < r */
174 : {
175 2451460 : if (g <= b)
176 : {
177 1659840 : minc = g;
178 1659840 : maxc = std::max(r, b);
179 : }
180 : else /* b < g */
181 : {
182 791616 : minc = b;
183 791616 : maxc = r;
184 : }
185 : }
186 4982560 : const int maxc_minus_minc = maxc - minc;
187 4982560 : if (s)
188 4982560 : *s = maxc_minus_minc / static_cast<float>(std::max(1, maxc));
189 4982560 : if (h)
190 : {
191 4982560 : const float maxc_minus_minc_times_6 =
192 4982560 : maxc_minus_minc == 0 ? 1.0f : 6.0f * maxc_minus_minc;
193 4982560 : if (maxc == b)
194 1700100 : *h = 4.0f / 6.0f + (r - g) / maxc_minus_minc_times_6;
195 3282460 : else if (maxc == g)
196 1660930 : *h = 2.0f / 6.0f + (b - r) / maxc_minus_minc_times_6;
197 : else
198 : {
199 1621540 : const float tmp = (g - b) / maxc_minus_minc_times_6;
200 1621540 : *h = tmp < 0.0f ? tmp + 1.0f : tmp;
201 : }
202 : }
203 4982560 : }
204 :
205 : /************************************************************************/
206 : /* choose_among() */
207 : /************************************************************************/
208 :
209 : template <typename T>
210 5509780 : static inline T choose_among(int idx, T a0, T a1, T a2, T a3, T a4, T a5)
211 : {
212 5509780 : switch (idx)
213 : {
214 917280 : case 0:
215 917280 : return a0;
216 918255 : case 1:
217 918255 : return a1;
218 919239 : case 2:
219 919239 : return a2;
220 918978 : case 3:
221 918978 : return a3;
222 918751 : case 4:
223 918751 : return a4;
224 917280 : default:
225 917280 : break;
226 : }
227 917280 : return a5;
228 : }
229 :
230 : /************************************************************************/
231 : /* hsv_to_rgb() */
232 : /************************************************************************/
233 :
234 : // hsv comes in as [h,s,v] with hue and saturation in the range [0,1],
235 : // but value in the range [0,255].
236 :
237 : // Derived from hsv_merge.py
238 :
239 4982560 : static void hsv_to_rgb(float h, float s, GByte v, GByte *r, GByte *g, GByte *b)
240 : {
241 4982560 : const int i = static_cast<int>(6.0f * h);
242 4982560 : const float f = 6.0f * h - i;
243 4982560 : const GByte p = static_cast<GByte>(v * (1.0f - s) + 0.5f);
244 4982560 : const GByte q = static_cast<GByte>(v * (1.0f - s * f) + 0.5f);
245 4982560 : const GByte t = static_cast<GByte>(v * (1.0f - s * (1.0f - f)) + 0.5f);
246 :
247 4982560 : if (r)
248 1836600 : *r = choose_among(i, v, q, p, p, t, v);
249 4982560 : if (g)
250 1836590 : *g = choose_among(i, t, v, v, q, p, p);
251 4982560 : if (b)
252 1836590 : *b = choose_among(i, p, p, t, v, v, q);
253 4982560 : }
254 :
255 : /************************************************************************/
256 : /* XMM_RGB_to_HS() */
257 : /************************************************************************/
258 :
259 : #ifdef HAVE_SSE2
260 : static inline void
261 538808 : XMM_RGB_to_HS(const GByte *CPL_RESTRICT pInR, const GByte *CPL_RESTRICT pInG,
262 : const GByte *CPL_RESTRICT pInB, const XMMReg4Float &zero,
263 : const XMMReg4Float &one, const XMMReg4Float &six,
264 : const XMMReg4Float &two_over_six,
265 : const XMMReg4Float &four_over_six, XMMReg4Float &h,
266 : XMMReg4Float &s)
267 : {
268 538808 : const auto r = XMMReg4Float::Load4Val(pInR);
269 538808 : const auto g = XMMReg4Float::Load4Val(pInG);
270 538808 : const auto b = XMMReg4Float::Load4Val(pInB);
271 538808 : const auto minc = XMMReg4Float::Min(XMMReg4Float::Min(r, g), b);
272 538808 : const auto maxc = XMMReg4Float::Max(XMMReg4Float::Max(r, g), b);
273 538808 : const auto max_minus_min = maxc - minc;
274 538808 : s = max_minus_min / XMMReg4Float::Max(one, maxc);
275 : const auto inv_max_minus_min_times_6_0 =
276 538808 : XMMReg4Float::Ternary(XMMReg4Float::Equals(max_minus_min, zero), one,
277 538808 : six * max_minus_min)
278 538808 : .inverse();
279 538808 : const auto tmp = (g - b) * inv_max_minus_min_times_6_0;
280 538808 : h = XMMReg4Float::Ternary(
281 538808 : XMMReg4Float::Equals(maxc, b),
282 538808 : four_over_six + (r - g) * inv_max_minus_min_times_6_0,
283 538808 : XMMReg4Float::Ternary(
284 538808 : XMMReg4Float::Equals(maxc, g),
285 538808 : two_over_six + (b - r) * inv_max_minus_min_times_6_0,
286 538808 : XMMReg4Float::Ternary(XMMReg4Float::Lesser(tmp, zero), tmp + one,
287 538808 : tmp)));
288 538808 : }
289 : #endif
290 :
291 : /************************************************************************/
292 : /* patch_value_line() */
293 : /************************************************************************/
294 :
295 : static
296 : #ifdef __GNUC__
297 : __attribute__((__noinline__))
298 : #endif
299 : void
300 695 : patch_value_line(int nCount, const GByte *CPL_RESTRICT pInR,
301 : const GByte *CPL_RESTRICT pInG,
302 : const GByte *CPL_RESTRICT pInB,
303 : const GByte *CPL_RESTRICT pInGray,
304 : GByte *CPL_RESTRICT pOutR, GByte *CPL_RESTRICT pOutG,
305 : GByte *CPL_RESTRICT pOutB)
306 : {
307 695 : int i = 0;
308 : #ifdef HAVE_SSE2
309 695 : const auto zero = XMMReg4Float::Zero();
310 695 : const auto one = XMMReg4Float::Set1(1.0f);
311 695 : const auto six = XMMReg4Float::Set1(6.0f);
312 695 : const auto two_over_six = XMMReg4Float::Set1(2.0f / 6.0f);
313 695 : const auto four_over_six = two_over_six + two_over_six;
314 :
315 695 : constexpr int ELTS = 8;
316 270099 : for (; i + (ELTS - 1) < nCount; i += ELTS)
317 : {
318 269404 : XMMReg4Float h0, s0;
319 269404 : XMM_RGB_to_HS(pInR + i, pInG + i, pInB + i, zero, one, six,
320 : two_over_six, four_over_six, h0, s0);
321 269404 : XMMReg4Float h1, s1;
322 269404 : XMM_RGB_to_HS(pInR + i + ELTS / 2, pInG + i + ELTS / 2,
323 269404 : pInB + i + ELTS / 2, zero, one, six, two_over_six,
324 : four_over_six, h1, s1);
325 :
326 269404 : XMMReg4Float v0, v1;
327 269404 : XMMReg4Float::Load8Val(pInGray + i, v0, v1);
328 :
329 269404 : const auto half = XMMReg4Float::Set1(0.5f);
330 269404 : const auto six_h0 = six * h0;
331 269404 : const auto idx0 = six_h0.truncate_to_int();
332 269404 : const auto f0 = six_h0 - idx0.cast_to_float();
333 269404 : const auto p0 = (v0 * (one - s0) + half).truncate_to_int();
334 269404 : const auto q0 = (v0 * (one - s0 * f0) + half).truncate_to_int();
335 269404 : const auto t0 = (v0 * (one - s0 * (one - f0)) + half).truncate_to_int();
336 :
337 269404 : const auto six_h1 = six * h1;
338 269404 : const auto idx1 = six_h1.truncate_to_int();
339 269404 : const auto f1 = six_h1 - idx1.cast_to_float();
340 269404 : const auto p1 = (v1 * (one - s1) + half).truncate_to_int();
341 269404 : const auto q1 = (v1 * (one - s1 * f1) + half).truncate_to_int();
342 269404 : const auto t1 = (v1 * (one - s1 * (one - f1)) + half).truncate_to_int();
343 :
344 269404 : const auto idx = XMMReg8Byte::Pack(idx0, idx1);
345 : const auto v =
346 269404 : XMMReg8Byte::Pack(v0.truncate_to_int(), v1.truncate_to_int());
347 269404 : const auto p = XMMReg8Byte::Pack(p0, p1);
348 269404 : const auto q = XMMReg8Byte::Pack(q0, q1);
349 269404 : const auto t = XMMReg8Byte::Pack(t0, t1);
350 :
351 269404 : const auto equalsTo0 = XMMReg8Byte::Equals(idx, XMMReg8Byte::Zero());
352 269404 : const auto one8Byte = XMMReg8Byte::Set1(1);
353 269404 : const auto equalsTo1 = XMMReg8Byte::Equals(idx, one8Byte);
354 269404 : const auto two8Byte = one8Byte + one8Byte;
355 269404 : const auto equalsTo2 = XMMReg8Byte::Equals(idx, two8Byte);
356 269404 : const auto four8Byte = two8Byte + two8Byte;
357 269404 : const auto equalsTo4 = XMMReg8Byte::Equals(idx, four8Byte);
358 269404 : const auto equalsTo3 = XMMReg8Byte::Equals(idx, four8Byte - one8Byte);
359 : // clang-format off
360 269404 : if (pOutR)
361 : {
362 : const auto out_r =
363 : XMMReg8Byte::Ternary(equalsTo0, v,
364 134702 : XMMReg8Byte::Ternary(equalsTo1, q,
365 134702 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo2, equalsTo3), p,
366 269404 : XMMReg8Byte::Ternary(equalsTo4, t, v))));
367 134702 : out_r.Store8Val(pOutR + i);
368 : }
369 269404 : if (pOutG)
370 : {
371 : const auto out_g =
372 : XMMReg8Byte::Ternary(equalsTo0, t,
373 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo1, equalsTo2), v,
374 236636 : XMMReg8Byte::Ternary(equalsTo3, q, p)));
375 118318 : out_g.Store8Val(pOutG + i);
376 : }
377 269404 : if (pOutB)
378 : {
379 : const auto out_b =
380 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo0, equalsTo1), p,
381 118318 : XMMReg8Byte::Ternary(equalsTo2, t,
382 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo3, equalsTo4),
383 118318 : v, q)));
384 118318 : out_b.Store8Val(pOutB + i);
385 : }
386 : // clang-format on
387 : }
388 : #endif
389 :
390 2524 : for (; i < nCount; ++i)
391 : {
392 : float h, s;
393 1829 : rgb_to_hs(pInR[i], pInG[i], pInB[i], &h, &s);
394 5001 : hsv_to_rgb(h, s, pInGray[i], pOutR ? pOutR + i : nullptr,
395 3172 : pOutG ? pOutG + i : nullptr, pOutB ? pOutB + i : nullptr);
396 : }
397 695 : }
398 :
399 : /************************************************************************/
400 : /* BlendBand */
401 : /************************************************************************/
402 :
403 : class BlendBand final : public GDALRasterBand
404 : {
405 : public:
406 120 : BlendBand(BlendDataset &oBlendDataset, int nBandIn)
407 120 : : m_oBlendDataset(oBlendDataset)
408 : {
409 120 : nBand = nBandIn;
410 120 : nRasterXSize = oBlendDataset.GetRasterXSize();
411 120 : nRasterYSize = oBlendDataset.GetRasterYSize();
412 120 : oBlendDataset.m_oColorDS.GetRasterBand(1)->GetBlockSize(&nBlockXSize,
413 : &nBlockYSize);
414 120 : eDataType = GDT_UInt8;
415 120 : }
416 :
417 10 : GDALColorInterp GetColorInterpretation() override
418 : {
419 10 : if (m_oBlendDataset.GetRasterCount() <= 2 && nBand == 1)
420 0 : return GCI_GrayIndex;
421 10 : else if (m_oBlendDataset.GetRasterCount() == 2 || nBand == 4)
422 0 : return GCI_AlphaBand;
423 : else
424 10 : return static_cast<GDALColorInterp>(GCI_RedBand + nBand - 1);
425 : }
426 :
427 18 : int GetOverviewCount() override
428 : {
429 18 : return static_cast<int>(m_oBlendDataset.m_apoOverviews.size());
430 : }
431 :
432 14 : GDALRasterBand *GetOverview(int idx) override
433 : {
434 13 : return idx >= 0 && idx < GetOverviewCount()
435 27 : ? m_oBlendDataset.m_apoOverviews[idx]->GetRasterBand(nBand)
436 14 : : nullptr;
437 : }
438 :
439 : protected:
440 62 : CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override
441 : {
442 62 : int nReqXSize = 0;
443 62 : int nReqYSize = 0;
444 62 : GetActualBlockSize(nBlockXOff, nBlockYOff, &nReqXSize, &nReqYSize);
445 124 : return RasterIO(GF_Read, nBlockXOff * nBlockXSize,
446 62 : nBlockYOff * nBlockYSize, nReqXSize, nReqYSize, pData,
447 62 : nReqXSize, nReqYSize, GDT_UInt8, 1, nBlockXSize,
448 124 : nullptr);
449 : }
450 :
451 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
452 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
453 : GDALDataType eBufType, GSpacing nPixelSpace,
454 : GSpacing nLineSpace,
455 : GDALRasterIOExtraArg *psExtraArg) override;
456 :
457 : private:
458 : BlendDataset &m_oBlendDataset;
459 : };
460 :
461 : /************************************************************************/
462 : /* BlendDataset::BlendDataset() */
463 : /************************************************************************/
464 :
465 41 : BlendDataset::BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,
466 41 : const std::string &sOperator, int nOpacity255Scale)
467 : : m_oColorDS(oColorDS), m_oOverlayDS(oOverlayDS), m_operator(sOperator),
468 41 : m_opacity255Scale(nOpacity255Scale)
469 : {
470 41 : m_oColorDS.Reference();
471 41 : m_oOverlayDS.Reference();
472 :
473 41 : CPLAssert(oColorDS.GetRasterXSize() == oOverlayDS.GetRasterXSize());
474 41 : CPLAssert(oColorDS.GetRasterYSize() == oOverlayDS.GetRasterYSize());
475 41 : nRasterXSize = oColorDS.GetRasterXSize();
476 41 : nRasterYSize = oColorDS.GetRasterYSize();
477 41 : const int nOvrCount = oOverlayDS.GetRasterBand(1)->GetOverviewCount();
478 41 : bool bCanCreateOvr = true;
479 161 : for (int iBand = 1; iBand <= oColorDS.GetRasterCount(); ++iBand)
480 : {
481 120 : SetBand(iBand, std::make_unique<BlendBand>(*this, iBand));
482 120 : bCanCreateOvr =
483 120 : bCanCreateOvr &&
484 120 : (iBand > oColorDS.GetRasterCount() ||
485 360 : oColorDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount) &&
486 120 : (iBand > oOverlayDS.GetRasterCount() ||
487 67 : oOverlayDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount);
488 : const int nColorBandxIdx =
489 120 : iBand <= oColorDS.GetRasterCount() ? iBand : 1;
490 : const int nOverlayBandIdx =
491 120 : iBand <= oOverlayDS.GetRasterCount() ? iBand : 1;
492 140 : for (int iOvr = 0; iOvr < nOvrCount && bCanCreateOvr; ++iOvr)
493 : {
494 : const auto poColorOvrBand =
495 20 : oColorDS.GetRasterBand(nColorBandxIdx)->GetOverview(iOvr);
496 : const auto poGSOvrBand =
497 20 : oOverlayDS.GetRasterBand(nOverlayBandIdx)->GetOverview(iOvr);
498 20 : bCanCreateOvr =
499 20 : poColorOvrBand->GetDataset() != &oColorDS &&
500 20 : poColorOvrBand->GetDataset() == oColorDS.GetRasterBand(1)
501 20 : ->GetOverview(iOvr)
502 20 : ->GetDataset() &&
503 20 : poGSOvrBand->GetDataset() != &oOverlayDS &&
504 20 : poGSOvrBand->GetDataset() == oOverlayDS.GetRasterBand(1)
505 20 : ->GetOverview(iOvr)
506 20 : ->GetDataset() &&
507 60 : poColorOvrBand->GetXSize() == poGSOvrBand->GetXSize() &&
508 20 : poColorOvrBand->GetYSize() == poGSOvrBand->GetYSize();
509 : }
510 : }
511 :
512 41 : SetDescription(CPLSPrintf("Blend %s width %s", m_oColorDS.GetDescription(),
513 41 : m_oOverlayDS.GetDescription()));
514 41 : if (nBands > 1)
515 : {
516 34 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
517 : }
518 :
519 41 : if (bCanCreateOvr)
520 : {
521 46 : for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
522 : {
523 5 : m_apoOverviews.push_back(std::make_unique<BlendDataset>(
524 5 : *(oColorDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),
525 5 : *(oOverlayDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),
526 5 : m_operator, nOpacity255Scale));
527 : }
528 : }
529 41 : }
530 :
531 : /************************************************************************/
532 : /* ~BlendDataset::BlendDataset() */
533 : /************************************************************************/
534 :
535 82 : BlendDataset::~BlendDataset()
536 : {
537 41 : m_oColorDS.ReleaseRef();
538 41 : m_oOverlayDS.ReleaseRef();
539 82 : }
540 :
541 : /************************************************************************/
542 : /* BlendDataset::AcquireSourcePixels() */
543 : /************************************************************************/
544 :
545 347 : bool BlendDataset::AcquireSourcePixels(int nXOff, int nYOff, int nXSize,
546 : int nYSize, int nBufXSize, int nBufYSize,
547 : GDALRasterIOExtraArg *psExtraArg)
548 : {
549 347 : if (nXOff == m_nCachedXOff && nYOff == m_nCachedYOff &&
550 92 : nXSize == m_nCachedXSize && nYSize == m_nCachedYSize &&
551 55 : nBufXSize == m_nCachedBufXSize && nBufYSize == m_nCachedBufYSize &&
552 55 : psExtraArg->eResampleAlg == m_sCachedExtraArg.eResampleAlg &&
553 55 : psExtraArg->bFloatingPointWindowValidity ==
554 55 : m_sCachedExtraArg.bFloatingPointWindowValidity &&
555 55 : (!psExtraArg->bFloatingPointWindowValidity ||
556 0 : (psExtraArg->dfXOff == m_sCachedExtraArg.dfXOff &&
557 0 : psExtraArg->dfYOff == m_sCachedExtraArg.dfYOff &&
558 0 : psExtraArg->dfXSize == m_sCachedExtraArg.dfXSize &&
559 0 : psExtraArg->dfYSize == m_sCachedExtraArg.dfYSize)))
560 : {
561 55 : return !m_abyBuffer.empty();
562 : }
563 :
564 292 : const int nColorCount = m_oColorDS.GetRasterCount();
565 292 : const int nOverlayCount = m_oOverlayDS.GetRasterCount();
566 292 : const int nCompsInBuffer = nColorCount + nOverlayCount;
567 292 : assert(nCompsInBuffer > 0);
568 :
569 584 : if (static_cast<size_t>(nBufXSize) >
570 292 : std::numeric_limits<size_t>::max() / nBufYSize / nCompsInBuffer)
571 : {
572 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
573 : "Out of memory allocating temporary buffer");
574 0 : m_abyBuffer.clear();
575 0 : m_ioError = true;
576 0 : return false;
577 : }
578 :
579 292 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
580 : try
581 : {
582 292 : if (m_abyBuffer.size() < nPixelCount * nCompsInBuffer)
583 36 : m_abyBuffer.resize(nPixelCount * nCompsInBuffer);
584 : }
585 4 : catch (const std::exception &)
586 : {
587 4 : CPLError(CE_Failure, CPLE_OutOfMemory,
588 : "Out of memory allocating temporary buffer");
589 4 : m_abyBuffer.clear();
590 4 : m_ioError = true;
591 4 : return false;
592 : }
593 :
594 : const bool bOK =
595 288 : (m_oColorDS.RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
596 288 : m_abyBuffer.data(), nBufXSize, nBufYSize,
597 : GDT_UInt8, nColorCount, nullptr, 1, nBufXSize,
598 : static_cast<GSpacing>(nPixelCount),
599 574 : psExtraArg) == CE_None &&
600 286 : m_oOverlayDS.RasterIO(
601 : GF_Read, nXOff, nYOff, nXSize, nYSize,
602 286 : m_abyBuffer.data() + nPixelCount * nColorCount, nBufXSize,
603 : nBufYSize, GDT_UInt8, nOverlayCount, nullptr, 1, nBufXSize,
604 288 : static_cast<GSpacing>(nPixelCount), psExtraArg) == CE_None);
605 288 : if (bOK)
606 : {
607 286 : m_nCachedXOff = nXOff;
608 286 : m_nCachedYOff = nYOff;
609 286 : m_nCachedXSize = nXSize;
610 286 : m_nCachedYSize = nYSize;
611 286 : m_nCachedBufXSize = nBufXSize;
612 286 : m_nCachedBufYSize = nBufYSize;
613 286 : m_sCachedExtraArg = *psExtraArg;
614 : }
615 : else
616 : {
617 2 : m_abyBuffer.clear();
618 2 : m_ioError = true;
619 : }
620 288 : return bOK;
621 : }
622 :
623 : /************************************************************************/
624 : /* gTabInvDstA */
625 : /************************************************************************/
626 :
627 : constexpr int SHIFT_DIV_DSTA = 8;
628 :
629 : // Table of (255 * 256 + k/2) / k values for k in [0,255]
630 : constexpr auto gTabInvDstA = []()
631 : {
632 : std::array<uint16_t, 256> arr{};
633 :
634 : arr[0] = 0;
635 : for (int k = 1; k <= 255; ++k)
636 : {
637 : arr[k] = static_cast<uint16_t>(((255 << SHIFT_DIV_DSTA) + (k / 2)) / k);
638 : }
639 :
640 : return arr;
641 : }();
642 :
643 : /************************************************************************/
644 : /* BlendSrcOverRGBA_SSE2() */
645 : /************************************************************************/
646 :
647 : #ifdef HAVE_SSE2
648 : #if defined(__GNUC__)
649 : __attribute__((noinline))
650 : #endif
651 : static int
652 1 : BlendSrcOverRGBA_SSE2(const GByte *CPL_RESTRICT pabyR,
653 : const GByte *CPL_RESTRICT pabyG,
654 : const GByte *CPL_RESTRICT pabyB,
655 : const GByte *CPL_RESTRICT pabyA,
656 : const GByte *CPL_RESTRICT pabyOverlayR,
657 : const GByte *CPL_RESTRICT pabyOverlayG,
658 : const GByte *CPL_RESTRICT pabyOverlayB,
659 : const GByte *CPL_RESTRICT pabyOverlayA,
660 : GByte *CPL_RESTRICT pabyDst, GSpacing nBandSpace, int N,
661 : GByte nOpacity)
662 : {
663 : // See scalar code after call to BlendSrcOverRGBA_SSE2() below for the
664 : // non-obfuscated formulas...
665 :
666 0 : const auto load_and_unpack = [](const void *p)
667 : {
668 0 : const auto zero = _mm_setzero_si128();
669 0 : auto overlayA = _mm_loadu_si128(reinterpret_cast<const __m128i *>(p));
670 0 : return std::make_pair(_mm_unpacklo_epi8(overlayA, zero),
671 0 : _mm_unpackhi_epi8(overlayA, zero));
672 : };
673 :
674 0 : const auto pack_and_store = [](void *p, __m128i lo, __m128i hi) {
675 0 : _mm_storeu_si128(reinterpret_cast<__m128i *>(p),
676 : _mm_packus_epi16(lo, hi));
677 0 : };
678 :
679 0 : const auto mul16bit_8bit_result = [](__m128i a, __m128i b)
680 : {
681 0 : const auto r255 = _mm_set1_epi16(255);
682 0 : return _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(a, b), r255), 8);
683 : };
684 :
685 2 : const auto opacity = _mm_set1_epi16(nOpacity);
686 1 : const auto r255 = _mm_set1_epi16(255);
687 : const int16_t *tabInvDstASigned =
688 1 : reinterpret_cast<const int16_t *>(gTabInvDstA.data());
689 1 : constexpr int REG_WIDTH = static_cast<int>(sizeof(opacity));
690 :
691 1 : int i = 0;
692 1 : for (; i <= N - REG_WIDTH; i += REG_WIDTH)
693 : {
694 0 : auto [overlayA_lo, overlayA_hi] = load_and_unpack(pabyOverlayA + i);
695 0 : auto [srcA_lo, srcA_hi] = load_and_unpack(pabyA + i);
696 0 : overlayA_lo = mul16bit_8bit_result(overlayA_lo, opacity);
697 0 : overlayA_hi = mul16bit_8bit_result(overlayA_hi, opacity);
698 : auto srcAMul255MinusOverlayA_lo =
699 0 : mul16bit_8bit_result(srcA_lo, _mm_sub_epi16(r255, overlayA_lo));
700 : auto srcAMul255MinusOverlayA_hi =
701 0 : mul16bit_8bit_result(srcA_hi, _mm_sub_epi16(r255, overlayA_hi));
702 0 : auto dstA_lo = _mm_add_epi16(overlayA_lo, srcAMul255MinusOverlayA_lo);
703 0 : auto dstA_hi = _mm_add_epi16(overlayA_hi, srcAMul255MinusOverlayA_hi);
704 :
705 : // The & 0xff should not be necessary. This is mostly a safety
706 : // belt if the above math yields a result outside [0, 255]...
707 0 : dstA_lo = _mm_and_si128(dstA_lo, r255);
708 0 : dstA_hi = _mm_and_si128(dstA_hi, r255);
709 :
710 : // This would be the equivalent of a "_mm_i16gather_epi16" operation
711 : // which does not exist...
712 : // invDstA_{i} = [tabInvDstASigned[dstA_{i}] for i in range(8)]
713 0 : auto invDstA_lo = _mm_undefined_si128();
714 0 : auto invDstA_hi = _mm_undefined_si128();
715 : #define SET_INVDSTA(k) \
716 : do \
717 : { \
718 : const int idxLo = _mm_extract_epi16(dstA_lo, k); \
719 : const int idxHi = _mm_extract_epi16(dstA_hi, k); \
720 : invDstA_lo = _mm_insert_epi16(invDstA_lo, tabInvDstASigned[idxLo], k); \
721 : invDstA_hi = _mm_insert_epi16(invDstA_hi, tabInvDstASigned[idxHi], k); \
722 : } while (0)
723 :
724 0 : SET_INVDSTA(0);
725 0 : SET_INVDSTA(1);
726 0 : SET_INVDSTA(2);
727 0 : SET_INVDSTA(3);
728 0 : SET_INVDSTA(4);
729 0 : SET_INVDSTA(5);
730 0 : SET_INVDSTA(6);
731 0 : SET_INVDSTA(7);
732 :
733 0 : pack_and_store(pabyDst + i + 3 * nBandSpace, dstA_lo, dstA_hi);
734 :
735 : #define PROCESS_COMPONENT(pabySrc, pabyOverlay, iBand) \
736 : do \
737 : { \
738 : auto [src_lo, src_hi] = load_and_unpack((pabySrc) + i); \
739 : auto [overlay_lo, overlay_hi] = load_and_unpack((pabyOverlay) + i); \
740 : auto dst_lo = _mm_srli_epi16( \
741 : _mm_add_epi16( \
742 : _mm_add_epi16( \
743 : _mm_mullo_epi16(overlay_lo, overlayA_lo), \
744 : _mm_mullo_epi16(src_lo, srcAMul255MinusOverlayA_lo)), \
745 : r255), \
746 : 8); \
747 : auto dst_hi = _mm_srli_epi16( \
748 : _mm_add_epi16( \
749 : _mm_add_epi16( \
750 : _mm_mullo_epi16(overlay_hi, overlayA_hi), \
751 : _mm_mullo_epi16(src_hi, srcAMul255MinusOverlayA_hi)), \
752 : r255), \
753 : 8); \
754 : dst_lo = mul16bit_8bit_result(dst_lo, invDstA_lo); \
755 : dst_hi = mul16bit_8bit_result(dst_hi, invDstA_hi); \
756 : pack_and_store(pabyDst + i + (iBand)*nBandSpace, dst_lo, dst_hi); \
757 : } while (0)
758 :
759 0 : PROCESS_COMPONENT(pabyR, pabyOverlayR, 0);
760 0 : PROCESS_COMPONENT(pabyG, pabyOverlayG, 1);
761 0 : PROCESS_COMPONENT(pabyB, pabyOverlayB, 2);
762 : }
763 1 : return i;
764 : }
765 : #endif
766 :
767 : template <bool bPixelSpaceIsOne>
768 : #if defined(__GNUC__)
769 : __attribute__((noinline))
770 : #endif
771 : static void
772 1 : BlendSrcOverRGBA_Generic(const GByte *CPL_RESTRICT pabyR,
773 : const GByte *CPL_RESTRICT pabyG,
774 : const GByte *CPL_RESTRICT pabyB,
775 : const GByte *CPL_RESTRICT pabyA,
776 : const GByte *CPL_RESTRICT pabyOverlayR,
777 : const GByte *CPL_RESTRICT pabyOverlayG,
778 : const GByte *CPL_RESTRICT pabyOverlayB,
779 : const GByte *CPL_RESTRICT pabyOverlayA,
780 : GByte *CPL_RESTRICT pabyDst,
781 : [[maybe_unused]] GSpacing nPixelSpace,
782 : GSpacing nBandSpace, int i, int N, GByte nOpacity)
783 : {
784 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
785 : if constexpr (!bPixelSpaceIsOne)
786 : {
787 : assert(nPixelSpace != 1);
788 : }
789 : #endif
790 1 : [[maybe_unused]] GSpacing nDstOffset = 0;
791 : #if defined(HAVE_SSE2) && defined(__clang__) && !defined(__INTEL_CLANG_COMPILER)
792 : // No need to vectorize for SSE2 and clang
793 : #pragma clang loop vectorize(disable)
794 : #endif
795 2 : for (; i < N; ++i)
796 : {
797 1 : const GByte nOverlayR = pabyOverlayR[i];
798 1 : const GByte nOverlayG = pabyOverlayG[i];
799 1 : const GByte nOverlayB = pabyOverlayB[i];
800 1 : const GByte nOverlayA =
801 1 : static_cast<GByte>((pabyOverlayA[i] * nOpacity + 255) / 256);
802 1 : const GByte nR = pabyR[i];
803 1 : const GByte nG = pabyG[i];
804 1 : const GByte nB = pabyB[i];
805 1 : const GByte nA = pabyA[i];
806 1 : const GByte nSrcAMul255MinusOverlayA =
807 1 : static_cast<GByte>((nA * (255 - nOverlayA) + 255) / 256);
808 1 : const GByte nDstA =
809 : static_cast<GByte>(nOverlayA + nSrcAMul255MinusOverlayA);
810 1 : GByte nDstR = static_cast<GByte>(
811 1 : (nOverlayR * nOverlayA + nR * nSrcAMul255MinusOverlayA + 255) /
812 : 256);
813 1 : GByte nDstG = static_cast<GByte>(
814 1 : (nOverlayG * nOverlayA + nG * nSrcAMul255MinusOverlayA + 255) /
815 : 256);
816 1 : GByte nDstB = static_cast<GByte>(
817 1 : (nOverlayB * nOverlayA + nB * nSrcAMul255MinusOverlayA + 255) /
818 : 256);
819 : // nInvDstA = (255 << SHIFT_DIV_DSTA) / nDstA;
820 1 : const uint16_t nInvDstA = gTabInvDstA[nDstA];
821 1 : constexpr unsigned ROUND_OFFSET_DIV_DSTA = ((1 << SHIFT_DIV_DSTA) - 1);
822 1 : nDstR = static_cast<GByte>((nDstR * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
823 : SHIFT_DIV_DSTA);
824 1 : nDstG = static_cast<GByte>((nDstG * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
825 : SHIFT_DIV_DSTA);
826 1 : nDstB = static_cast<GByte>((nDstB * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
827 : SHIFT_DIV_DSTA);
828 : if constexpr (bPixelSpaceIsOne)
829 : {
830 : pabyDst[i] = nDstR;
831 : pabyDst[i + nBandSpace] = nDstG;
832 : pabyDst[i + 2 * nBandSpace] = nDstB;
833 : pabyDst[i + 3 * nBandSpace] = nDstA;
834 : }
835 : else
836 : {
837 1 : pabyDst[nDstOffset] = nDstR;
838 1 : pabyDst[nDstOffset + nBandSpace] = nDstG;
839 1 : pabyDst[nDstOffset + 2 * nBandSpace] = nDstB;
840 1 : pabyDst[nDstOffset + 3 * nBandSpace] = nDstA;
841 1 : nDstOffset += nPixelSpace;
842 : }
843 : }
844 1 : }
845 :
846 : /************************************************************************/
847 : /* BlendDataset::IRasterIO() */
848 : /************************************************************************/
849 :
850 229 : CPLErr BlendDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
851 : int nXSize, int nYSize, void *pData,
852 : int nBufXSize, int nBufYSize,
853 : GDALDataType eBufType, int nBandCount,
854 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
855 : GSpacing nLineSpace, GSpacing nBandSpace,
856 : GDALRasterIOExtraArg *psExtraArg)
857 : {
858 : // Try to pass the request to the most appropriate overview dataset.
859 229 : if (nBufXSize < nXSize && nBufYSize < nYSize)
860 : {
861 2 : int bTried = FALSE;
862 2 : const CPLErr eErr = TryOverviewRasterIO(
863 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
864 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
865 : nBandSpace, psExtraArg, &bTried);
866 2 : if (bTried)
867 2 : return eErr;
868 : }
869 :
870 227 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
871 227 : const int nColorCount = m_oColorDS.GetRasterCount();
872 227 : const int nOverlayCount = m_oOverlayDS.GetRasterCount();
873 419 : if (nOverlayCount == 1 && m_opacity255Scale == 255 &&
874 410 : m_operator == HSV_VALUE && eRWFlag == GF_Read &&
875 203 : eBufType == GDT_UInt8 && nBandCount == nBands &&
876 643 : IsAllBands(nBands, panBandMap) &&
877 203 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
878 : psExtraArg))
879 : {
880 197 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
881 197 : const GByte *pabyR = m_abyBuffer.data();
882 197 : const GByte *pabyG = m_abyBuffer.data() + nPixelCount;
883 197 : const GByte *pabyB = m_abyBuffer.data() + nPixelCount * 2;
884 197 : const GByte *pabyValue = m_abyBuffer.data() + nPixelCount * nColorCount;
885 197 : size_t nSrcIdx = 0;
886 516 : for (int j = 0; j < nBufYSize; ++j)
887 : {
888 319 : auto nDstOffset = j * nLineSpace;
889 319 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize &&
890 317 : nBandSpace >= nLineSpace * nBufYSize)
891 : {
892 317 : patch_value_line(nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
893 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
894 317 : pabyDst + nDstOffset,
895 317 : pabyDst + nDstOffset + nBandSpace,
896 317 : pabyDst + nDstOffset + 2 * nBandSpace);
897 317 : nSrcIdx += nBufXSize;
898 : }
899 : else
900 : {
901 262146 : for (int i = 0; i < nBufXSize;
902 262144 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
903 : {
904 : float h, s;
905 262144 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
906 : &h, &s);
907 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
908 262144 : &pabyDst[nDstOffset + 0 * nBandSpace],
909 262144 : &pabyDst[nDstOffset + 1 * nBandSpace],
910 262144 : &pabyDst[nDstOffset + 2 * nBandSpace]);
911 : }
912 : }
913 : }
914 197 : if (nColorCount == 4)
915 : {
916 394 : for (int j = 0; j < nBufYSize; ++j)
917 : {
918 198 : auto nDstOffset = 3 * nBandSpace + j * nLineSpace;
919 198 : const GByte *pabyA = m_abyBuffer.data() + nPixelCount * 3;
920 198 : GDALCopyWords64(pabyA, GDT_UInt8, 1, pabyDst + nDstOffset,
921 : GDT_UInt8, static_cast<int>(nPixelSpace),
922 : nBufXSize);
923 : }
924 : }
925 :
926 197 : return CE_None;
927 : }
928 3 : else if (nOverlayCount == 4 && nColorCount == 4 && m_operator == SRC_OVER &&
929 1 : eRWFlag == GF_Read && eBufType == GDT_UInt8 &&
930 34 : nBandCount == nBands && IsAllBands(nBands, panBandMap) &&
931 1 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize,
932 : nBufYSize, psExtraArg))
933 : {
934 1 : const GByte nOpacity = static_cast<GByte>(m_opacity255Scale);
935 1 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
936 1 : const GByte *CPL_RESTRICT pabyR = m_abyBuffer.data();
937 1 : const GByte *CPL_RESTRICT pabyG = m_abyBuffer.data() + nPixelCount;
938 1 : const GByte *CPL_RESTRICT pabyB = m_abyBuffer.data() + nPixelCount * 2;
939 1 : const GByte *CPL_RESTRICT pabyA = m_abyBuffer.data() + nPixelCount * 3;
940 : const GByte *CPL_RESTRICT pabyOverlayR =
941 1 : m_abyBuffer.data() + nPixelCount * nColorCount;
942 : const GByte *CPL_RESTRICT pabyOverlayG =
943 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
944 : const GByte *CPL_RESTRICT pabyOverlayB =
945 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 2);
946 : const GByte *CPL_RESTRICT pabyOverlayA =
947 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 3);
948 1 : size_t nSrcIdx = 0;
949 2 : for (int j = 0; j < nBufYSize; ++j)
950 : {
951 1 : auto nDstOffset = j * nLineSpace;
952 1 : int i = 0;
953 : #ifdef HAVE_SSE2
954 1 : if (nPixelSpace == 1)
955 : {
956 2 : i = BlendSrcOverRGBA_SSE2(
957 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
958 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
959 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
960 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nBandSpace,
961 : nBufXSize, nOpacity);
962 1 : nSrcIdx += i;
963 1 : nDstOffset += i;
964 : }
965 : #endif
966 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
967 : if (nPixelSpace == 1)
968 : {
969 : // Note: clang 20 does a rather good job at autovectorizing
970 : // for SSE2, but BlendSrcOverRGBA_SSE2() performs better.
971 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ true>(
972 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
973 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
974 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
975 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset,
976 : /*nPixelSpace = */ 1, nBandSpace, i, nBufXSize, nOpacity);
977 : }
978 : else
979 : #endif
980 : {
981 1 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ false>(
982 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
983 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
984 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
985 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nPixelSpace,
986 : nBandSpace, i, nBufXSize, nOpacity);
987 : }
988 1 : nSrcIdx += nBufXSize - i;
989 : }
990 1 : return CE_None;
991 : }
992 29 : else if (m_ioError)
993 : {
994 6 : return CE_Failure;
995 : }
996 : else
997 : {
998 23 : const CPLErr eErr = GDALDataset::IRasterIO(
999 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1000 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1001 : nBandSpace, psExtraArg);
1002 23 : m_ioError = eErr != CE_None;
1003 23 : return eErr;
1004 : }
1005 : }
1006 :
1007 : /************************************************************************/
1008 : /* SrcOverRGBOneComponent() */
1009 : /************************************************************************/
1010 :
1011 : // GCC and clang do a god job a auto vectorizing the below function
1012 : #if defined(__GNUC__) && !defined(__clang__)
1013 : __attribute__((optimize("tree-vectorize")))
1014 : #endif
1015 : static void
1016 12 : SrcOverRGB(const uint8_t *const __restrict pabyOverlay,
1017 : const uint8_t *const __restrict pabySrc,
1018 : uint8_t *const __restrict pabyDst, const size_t N,
1019 : const uint8_t nOpacity)
1020 : {
1021 915 : for (size_t i = 0; i < N; ++i)
1022 : {
1023 903 : const uint8_t nOverlay = pabyOverlay[i];
1024 903 : const uint8_t nSrc = pabySrc[i];
1025 903 : pabyDst[i] = static_cast<uint8_t>(
1026 903 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) / 256);
1027 : }
1028 12 : }
1029 :
1030 : /************************************************************************/
1031 : /* BlendBand::IRasterIO() */
1032 : /************************************************************************/
1033 :
1034 212 : CPLErr BlendBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
1035 : int nXSize, int nYSize, void *pData, int nBufXSize,
1036 : int nBufYSize, GDALDataType eBufType,
1037 : GSpacing nPixelSpace, GSpacing nLineSpace,
1038 : GDALRasterIOExtraArg *psExtraArg)
1039 : {
1040 : // Try to pass the request to the most appropriate overview dataset.
1041 212 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1042 : {
1043 1 : int bTried = FALSE;
1044 1 : const CPLErr eErr = TryOverviewRasterIO(
1045 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1046 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
1047 1 : if (bTried)
1048 1 : return eErr;
1049 : }
1050 :
1051 211 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1052 211 : const int nColorCount = m_oBlendDataset.m_oColorDS.GetRasterCount();
1053 211 : const int nOverlayCount = m_oBlendDataset.m_oOverlayDS.GetRasterCount();
1054 211 : if (nBand == 4 && m_oBlendDataset.m_operator == HSV_VALUE)
1055 : {
1056 9 : if (nColorCount == 3)
1057 : {
1058 0 : GByte ch = 255;
1059 0 : for (int iY = 0; iY < nBufYSize; ++iY)
1060 : {
1061 0 : GDALCopyWords64(&ch, GDT_UInt8, 0,
1062 0 : static_cast<GByte *>(pData) + iY * nLineSpace,
1063 : eBufType, static_cast<int>(nPixelSpace),
1064 : nBufXSize);
1065 : }
1066 0 : return CE_None;
1067 : }
1068 : else
1069 : {
1070 9 : CPLAssert(nColorCount == 4);
1071 9 : return m_oBlendDataset.m_oColorDS.GetRasterBand(4)->RasterIO(
1072 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
1073 9 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg);
1074 : }
1075 : }
1076 43 : else if (nOverlayCount == 3 && nColorCount == 3 &&
1077 30 : m_oBlendDataset.m_operator == SRC_OVER && eRWFlag == GF_Read &&
1078 227 : eBufType == GDT_UInt8 &&
1079 12 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
1080 : nBufXSize, nBufYSize,
1081 : psExtraArg))
1082 : {
1083 12 : const int nOpacity = m_oBlendDataset.m_opacity255Scale;
1084 : const GByte *const CPL_RESTRICT pabySrc =
1085 12 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * (nBand - 1);
1086 : const GByte *const CPL_RESTRICT pabyOverlay =
1087 12 : m_oBlendDataset.m_abyBuffer.data() +
1088 12 : nPixelCount * (nColorCount + nBand - 1);
1089 12 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
1090 12 : size_t nSrcIdx = 0;
1091 24 : for (int j = 0; j < nBufYSize; ++j)
1092 : {
1093 12 : auto nDstOffset = j * nLineSpace;
1094 12 : if (nPixelSpace == 1)
1095 : {
1096 12 : SrcOverRGB(pabyOverlay + nSrcIdx, pabySrc + nSrcIdx,
1097 12 : pabyDst + nDstOffset, nBufXSize,
1098 : static_cast<uint8_t>(nOpacity));
1099 12 : nSrcIdx += nBufXSize;
1100 : }
1101 : else
1102 : {
1103 0 : for (int i = 0; i < nBufXSize;
1104 0 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1105 : {
1106 0 : const int nOverlay = pabyOverlay[nSrcIdx];
1107 0 : const int nSrc = pabySrc[nSrcIdx];
1108 0 : pabyDst[nDstOffset] = static_cast<GByte>(
1109 0 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) /
1110 : 256);
1111 : }
1112 : }
1113 : }
1114 12 : return CE_None;
1115 : }
1116 321 : else if (eRWFlag == GF_Read && eBufType == GDT_UInt8 &&
1117 131 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
1118 : nBufXSize, nBufYSize,
1119 : psExtraArg))
1120 : {
1121 131 : GByte *pabyDst = static_cast<GByte *>(pData);
1122 131 : if (m_oBlendDataset.m_operator == SRC_OVER)
1123 : {
1124 3 : const auto RGBToGrayScale = [](int R, int G, int B)
1125 : {
1126 : // Equivalent to R * 0.299 + G * 0.587 + B * 0.114
1127 : // but using faster computation
1128 3 : return (R * 306 + G * 601 + B * 117) / 1024;
1129 : };
1130 :
1131 : const GByte *paby =
1132 83 : (nBand <= nColorCount) ? m_oBlendDataset.m_abyBuffer.data() +
1133 83 : nPixelCount * (nBand - 1)
1134 0 : : (nBand == 4 && nColorCount == 2)
1135 0 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
1136 83 : : nullptr;
1137 : const GByte *pabyA =
1138 : (nColorCount == 4)
1139 87 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 3
1140 : : nColorCount == 2
1141 4 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
1142 83 : : nullptr;
1143 : const GByte *pabyOverlay =
1144 83 : (nBand <= nOverlayCount)
1145 91 : ? m_oBlendDataset.m_abyBuffer.data() +
1146 75 : nPixelCount * (nColorCount + nBand - 1)
1147 8 : : (nBand <= 3) ? m_oBlendDataset.m_abyBuffer.data() +
1148 7 : nPixelCount * nColorCount
1149 83 : : nullptr;
1150 : const GByte *pabyOverlayA =
1151 80 : (nOverlayCount == 2 || nOverlayCount == 4)
1152 163 : ? m_oBlendDataset.m_abyBuffer.data() +
1153 62 : nPixelCount * (nColorCount + nOverlayCount - 1)
1154 83 : : nullptr;
1155 : const GByte *pabyOverlayR =
1156 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1157 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1158 3 : nPixelCount * nColorCount
1159 83 : : nullptr;
1160 : const GByte *pabyOverlayG =
1161 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1162 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1163 3 : nPixelCount * (nColorCount + 1)
1164 83 : : nullptr;
1165 : const GByte *pabyOverlayB =
1166 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1167 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1168 3 : nPixelCount * (nColorCount + 2)
1169 83 : : nullptr;
1170 :
1171 83 : size_t nSrcIdx = 0;
1172 714 : for (int j = 0; j < nBufYSize; ++j)
1173 : {
1174 631 : auto nDstOffset = j * nLineSpace;
1175 98753 : for (int i = 0; i < nBufXSize;
1176 98122 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1177 : {
1178 : // Corrected to take into account m_opacity255Scale
1179 98122 : const int nOverlayA =
1180 98122 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
1181 97606 : m_oBlendDataset.m_opacity255Scale +
1182 : 255) /
1183 : 256)
1184 516 : : m_oBlendDataset.m_opacity255Scale;
1185 :
1186 98122 : const int nSrcA = pabyA ? pabyA[nSrcIdx] : 255;
1187 :
1188 98122 : const int nSrcAMul255MinusOverlayA =
1189 98122 : (nSrcA * (255 - nOverlayA) + 255) / 256;
1190 98122 : const int nDstA = nOverlayA + nSrcAMul255MinusOverlayA;
1191 98122 : if (nBand != 4)
1192 : {
1193 : const int nOverlay =
1194 3 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
1195 147439 : ? RGBToGrayScale(pabyOverlayR[nSrcIdx],
1196 3 : pabyOverlayG[nSrcIdx],
1197 3 : pabyOverlayB[nSrcIdx])
1198 73718 : : pabyOverlay ? pabyOverlay[nSrcIdx]
1199 73721 : : 255;
1200 :
1201 73721 : const int nSrc = paby ? paby[nSrcIdx] : 255;
1202 73721 : int nDst = (nOverlay * nOverlayA +
1203 73721 : nSrc * nSrcAMul255MinusOverlayA + 255) /
1204 : 256;
1205 73721 : if (nDstA != 0 && nDstA != 255)
1206 19672 : nDst = (nDst * 255 + (nDstA / 2)) / nDstA;
1207 73721 : pabyDst[nDstOffset] =
1208 73721 : static_cast<GByte>(std::min(nDst, 255));
1209 : }
1210 : else
1211 : {
1212 24401 : pabyDst[nDstOffset] = static_cast<GByte>(nDstA);
1213 : }
1214 : }
1215 : }
1216 : }
1217 48 : else if (nOverlayCount == 1 && m_oBlendDataset.m_opacity255Scale == 255)
1218 : {
1219 18 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
1220 : const GByte *pabyG =
1221 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
1222 : const GByte *pabyB =
1223 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
1224 18 : CPLAssert(m_oBlendDataset.m_operator == HSV_VALUE);
1225 18 : size_t nSrcIdx = 0;
1226 : const GByte *pabyValue =
1227 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
1228 402 : for (int j = 0; j < nBufYSize; ++j)
1229 : {
1230 384 : auto nDstOffset = j * nLineSpace;
1231 384 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize)
1232 : {
1233 884 : patch_value_line(
1234 : nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
1235 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
1236 378 : nBand == 1 ? pabyDst + nDstOffset : nullptr,
1237 378 : nBand == 2 ? pabyDst + nDstOffset : nullptr,
1238 378 : nBand == 3 ? pabyDst + nDstOffset : nullptr);
1239 378 : nSrcIdx += nBufXSize;
1240 : }
1241 : else
1242 : {
1243 786438 : for (int i = 0; i < nBufXSize;
1244 786432 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1245 : {
1246 : float h, s;
1247 786432 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx],
1248 786432 : pabyB[nSrcIdx], &h, &s);
1249 786432 : if (nBand == 1)
1250 : {
1251 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
1252 262144 : &pabyDst[nDstOffset], nullptr, nullptr);
1253 : }
1254 524288 : else if (nBand == 2)
1255 : {
1256 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
1257 262144 : &pabyDst[nDstOffset], nullptr);
1258 : }
1259 : else
1260 : {
1261 262144 : CPLAssert(nBand == 3);
1262 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
1263 262144 : nullptr, &pabyDst[nDstOffset]);
1264 : }
1265 : }
1266 : }
1267 18 : }
1268 : }
1269 : else
1270 : {
1271 30 : CPLAssert(m_oBlendDataset.m_operator == HSV_VALUE);
1272 30 : CPLAssert(nBand <= 3);
1273 30 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
1274 : const GByte *pabyG =
1275 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
1276 : const GByte *pabyB =
1277 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
1278 : const GByte *pabyValue =
1279 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
1280 : const GByte *pabyOverlayR =
1281 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1282 12 : nPixelCount * nColorCount
1283 30 : : nullptr;
1284 : const GByte *pabyOverlayG =
1285 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1286 12 : nPixelCount * (nColorCount + 1)
1287 30 : : nullptr;
1288 : const GByte *pabyOverlayB =
1289 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1290 12 : nPixelCount * (nColorCount + 2)
1291 30 : : nullptr;
1292 : const GByte *pabyOverlayA =
1293 24 : (nOverlayCount == 2 || nOverlayCount == 4)
1294 54 : ? m_oBlendDataset.m_abyBuffer.data() +
1295 12 : nPixelCount * (nColorCount + nOverlayCount - 1)
1296 30 : : nullptr;
1297 :
1298 30 : size_t nSrcIdx = 0;
1299 60 : for (int j = 0; j < nBufYSize; ++j)
1300 : {
1301 30 : auto nDstOffset = j * nLineSpace;
1302 3932190 : for (int i = 0; i < nBufXSize;
1303 3932160 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1304 : {
1305 3932160 : const int nColorR = pabyR[nSrcIdx];
1306 3932160 : const int nColorG = pabyG[nSrcIdx];
1307 3932160 : const int nColorB = pabyB[nSrcIdx];
1308 : const int nOverlayV =
1309 1572860 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
1310 5505020 : ? std::max({pabyOverlayR[nSrcIdx],
1311 1572860 : pabyOverlayG[nSrcIdx],
1312 1572860 : pabyOverlayB[nSrcIdx]})
1313 2359300 : : pabyValue[nSrcIdx];
1314 3932160 : const int nOverlayA =
1315 3932160 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
1316 1572860 : m_oBlendDataset.m_opacity255Scale +
1317 : 255) /
1318 : 256)
1319 2359300 : : m_oBlendDataset.m_opacity255Scale;
1320 : const int nColorValue =
1321 3932160 : std::max({nColorR, nColorG, nColorB});
1322 :
1323 : float h, s;
1324 3932160 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
1325 : &h, &s);
1326 :
1327 3932160 : const GByte nTargetValue = static_cast<GByte>(
1328 3932160 : (nOverlayV * nOverlayA +
1329 3932160 : nColorValue * (255 - nOverlayA) + 255) /
1330 : 256);
1331 :
1332 3932160 : if (nBand == 1)
1333 : {
1334 1310720 : hsv_to_rgb(h, s, nTargetValue, &pabyDst[nDstOffset],
1335 : nullptr, nullptr);
1336 : }
1337 2621440 : else if (nBand == 2)
1338 : {
1339 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr,
1340 1310720 : &pabyDst[nDstOffset], nullptr);
1341 : }
1342 : else
1343 : {
1344 1310720 : CPLAssert(nBand == 3);
1345 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr, nullptr,
1346 1310720 : &pabyDst[nDstOffset]);
1347 : }
1348 : }
1349 : }
1350 : }
1351 :
1352 131 : return CE_None;
1353 : }
1354 59 : else if (m_oBlendDataset.m_ioError)
1355 : {
1356 0 : return CE_Failure;
1357 : }
1358 : else
1359 : {
1360 59 : const CPLErr eErr = GDALRasterBand::IRasterIO(
1361 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1362 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
1363 59 : m_oBlendDataset.m_ioError = eErr != CE_None;
1364 59 : return eErr;
1365 : }
1366 : }
1367 :
1368 : } // namespace
1369 :
1370 : /************************************************************************/
1371 : /* GDALRasterBlendAlgorithm::ValidateGlobal() */
1372 : /************************************************************************/
1373 :
1374 80 : bool GDALRasterBlendAlgorithm::ValidateGlobal()
1375 : {
1376 : auto poSrcDS =
1377 80 : m_inputDataset.empty() ? nullptr : m_inputDataset[0].GetDatasetRef();
1378 80 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
1379 80 : if (poSrcDS)
1380 : {
1381 154 : if (poSrcDS->GetRasterCount() == 0 || poSrcDS->GetRasterCount() > 4 ||
1382 77 : poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_UInt8)
1383 : {
1384 1 : ReportError(CE_Failure, CPLE_NotSupported,
1385 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
1386 : "supported as input");
1387 1 : return false;
1388 : }
1389 : }
1390 79 : if (poOverlayDS)
1391 : {
1392 79 : if (poOverlayDS->GetRasterCount() == 0 ||
1393 158 : poOverlayDS->GetRasterCount() > 4 ||
1394 79 : poOverlayDS->GetRasterBand(1)->GetRasterDataType() != GDT_UInt8)
1395 : {
1396 1 : ReportError(CE_Failure, CPLE_NotSupported,
1397 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
1398 : "supported as overlay");
1399 1 : return false;
1400 : }
1401 : }
1402 :
1403 78 : if (poSrcDS && poOverlayDS)
1404 : {
1405 149 : if (poSrcDS->GetRasterXSize() != poOverlayDS->GetRasterXSize() ||
1406 74 : poSrcDS->GetRasterYSize() != poOverlayDS->GetRasterYSize())
1407 : {
1408 2 : ReportError(CE_Failure, CPLE_IllegalArg,
1409 : "Input dataset and overlay dataset must have "
1410 : "the same dimensions");
1411 2 : return false;
1412 : }
1413 :
1414 84 : if (m_operator == HSV_VALUE && poSrcDS->GetRasterCount() != 3 &&
1415 11 : poSrcDS->GetRasterCount() != 4)
1416 : {
1417 1 : ReportError(CE_Failure, CPLE_AppDefined,
1418 : "Operator %s requires a 3-band or 4-band input dataset",
1419 : HSV_VALUE);
1420 1 : return false;
1421 : }
1422 : }
1423 :
1424 75 : return true;
1425 : }
1426 :
1427 : /************************************************************************/
1428 : /* GDALRasterBlendAlgorithm::RunStep() */
1429 : /************************************************************************/
1430 :
1431 36 : bool GDALRasterBlendAlgorithm::RunStep(GDALPipelineStepRunContext &)
1432 : {
1433 36 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
1434 36 : CPLAssert(poSrcDS);
1435 :
1436 36 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
1437 36 : CPLAssert(poOverlayDS);
1438 :
1439 36 : if (!ValidateGlobal())
1440 0 : return false;
1441 :
1442 36 : const int nOpacity255Scale =
1443 36 : (m_opacity * 255 + OPACITY_INPUT_RANGE / 2) / OPACITY_INPUT_RANGE;
1444 :
1445 36 : m_outputDataset.Set(std::make_unique<BlendDataset>(
1446 36 : *poSrcDS, *poOverlayDS, m_operator, nOpacity255Scale));
1447 :
1448 36 : return true;
1449 : }
1450 :
1451 : GDALRasterBlendAlgorithmStandalone::~GDALRasterBlendAlgorithmStandalone() =
1452 : default;
1453 :
1454 : //! @endcond
|
