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 84 : GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm(bool standaloneStep)
47 : : GDALRasterPipelineStepAlgorithm(
48 : NAME, DESCRIPTION, HELP_URL,
49 0 : ConstructorOptions()
50 84 : .SetStandaloneStep(standaloneStep)
51 84 : .SetAddDefaultArguments(false)
52 168 : .SetInputDatasetHelpMsg(_("Input raster dataset"))
53 168 : .SetInputDatasetAlias("color-input")
54 168 : .SetInputDatasetMetaVar("COLOR-INPUT")
55 252 : .SetOutputDatasetHelpMsg(_("Output raster dataset")))
56 : {
57 84 : const auto AddOverlayDatasetArg = [this]()
58 : {
59 : auto &arg = AddArg("overlay", 0, _("Overlay dataset"),
60 168 : &m_overlayDataset, GDAL_OF_RASTER)
61 84 : .SetPositional()
62 84 : .SetRequired();
63 :
64 84 : SetAutoCompleteFunctionForFilename(arg, GDAL_OF_RASTER);
65 84 : };
66 :
67 84 : if (standaloneStep)
68 : {
69 42 : AddRasterInputArgs(false, false);
70 42 : AddOverlayDatasetArg();
71 42 : AddProgressArg();
72 42 : AddRasterOutputArgs(false);
73 : }
74 : else
75 : {
76 42 : AddRasterHiddenInputDatasetArg();
77 42 : AddOverlayDatasetArg();
78 : }
79 :
80 168 : AddArg("operator", 0, _("Composition operator"), &m_operator)
81 84 : .SetChoices(SRC_OVER, HSV_VALUE)
82 84 : .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 168 : &m_opacity)
87 84 : .SetDefault(m_opacity)
88 84 : .SetMinValueIncluded(0)
89 84 : .SetMaxValueIncluded(OPACITY_INPUT_RANGE);
90 :
91 128 : AddValidationAction([this]() { return ValidateGlobal(); });
92 84 : }
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.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.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_Byte;
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_Byte, 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(
596 288 : GF_Read, nXOff, nYOff, nXSize, nYSize, m_abyBuffer.data(),
597 : nBufXSize, nBufYSize, GDT_Byte, nColorCount, nullptr, 1, nBufXSize,
598 574 : static_cast<GSpacing>(nPixelCount), psExtraArg) == CE_None &&
599 286 : m_oOverlayDS.RasterIO(
600 : GF_Read, nXOff, nYOff, nXSize, nYSize,
601 286 : m_abyBuffer.data() + nPixelCount * nColorCount, nBufXSize,
602 : nBufYSize, GDT_Byte, nOverlayCount, nullptr, 1, nBufXSize,
603 288 : static_cast<GSpacing>(nPixelCount), psExtraArg) == CE_None);
604 288 : if (bOK)
605 : {
606 286 : m_nCachedXOff = nXOff;
607 286 : m_nCachedYOff = nYOff;
608 286 : m_nCachedXSize = nXSize;
609 286 : m_nCachedYSize = nYSize;
610 286 : m_nCachedBufXSize = nBufXSize;
611 286 : m_nCachedBufYSize = nBufYSize;
612 286 : m_sCachedExtraArg = *psExtraArg;
613 : }
614 : else
615 : {
616 2 : m_abyBuffer.clear();
617 2 : m_ioError = true;
618 : }
619 288 : return bOK;
620 : }
621 :
622 : /************************************************************************/
623 : /* gTabInvDstA */
624 : /************************************************************************/
625 :
626 : constexpr int SHIFT_DIV_DSTA = 8;
627 :
628 : // Table of (255 * 256 + k/2) / k values for k in [0,255]
629 : constexpr auto gTabInvDstA = []()
630 : {
631 : std::array<uint16_t, 256> arr{};
632 :
633 : arr[0] = 0;
634 : for (int k = 1; k <= 255; ++k)
635 : {
636 : arr[k] = static_cast<uint16_t>(((255 << SHIFT_DIV_DSTA) + (k / 2)) / k);
637 : }
638 :
639 : return arr;
640 : }();
641 :
642 : /************************************************************************/
643 : /* BlendSrcOverRGBA_SSE2() */
644 : /************************************************************************/
645 :
646 : #ifdef HAVE_SSE2
647 : #if defined(__GNUC__)
648 : __attribute__((noinline))
649 : #endif
650 : static int
651 1 : BlendSrcOverRGBA_SSE2(const GByte *CPL_RESTRICT pabyR,
652 : const GByte *CPL_RESTRICT pabyG,
653 : const GByte *CPL_RESTRICT pabyB,
654 : const GByte *CPL_RESTRICT pabyA,
655 : const GByte *CPL_RESTRICT pabyOverlayR,
656 : const GByte *CPL_RESTRICT pabyOverlayG,
657 : const GByte *CPL_RESTRICT pabyOverlayB,
658 : const GByte *CPL_RESTRICT pabyOverlayA,
659 : GByte *CPL_RESTRICT pabyDst, GSpacing nBandSpace, int N,
660 : GByte nOpacity)
661 : {
662 : // See scalar code after call to BlendSrcOverRGBA_SSE2() below for the
663 : // non-obfuscated formulas...
664 :
665 0 : const auto load_and_unpack = [](const void *p)
666 : {
667 0 : const auto zero = _mm_setzero_si128();
668 0 : auto overlayA = _mm_loadu_si128(reinterpret_cast<const __m128i *>(p));
669 0 : return std::make_pair(_mm_unpacklo_epi8(overlayA, zero),
670 0 : _mm_unpackhi_epi8(overlayA, zero));
671 : };
672 :
673 0 : const auto pack_and_store = [](void *p, __m128i lo, __m128i hi) {
674 0 : _mm_storeu_si128(reinterpret_cast<__m128i *>(p),
675 : _mm_packus_epi16(lo, hi));
676 0 : };
677 :
678 0 : const auto mul16bit_8bit_result = [](__m128i a, __m128i b)
679 : {
680 0 : const auto r255 = _mm_set1_epi16(255);
681 0 : return _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(a, b), r255), 8);
682 : };
683 :
684 2 : const auto opacity = _mm_set1_epi16(nOpacity);
685 1 : const auto r255 = _mm_set1_epi16(255);
686 : const int16_t *tabInvDstASigned =
687 1 : reinterpret_cast<const int16_t *>(gTabInvDstA.data());
688 1 : constexpr int REG_WIDTH = static_cast<int>(sizeof(opacity));
689 :
690 1 : int i = 0;
691 1 : for (; i <= N - REG_WIDTH; i += REG_WIDTH)
692 : {
693 0 : auto [overlayA_lo, overlayA_hi] = load_and_unpack(pabyOverlayA + i);
694 0 : auto [srcA_lo, srcA_hi] = load_and_unpack(pabyA + i);
695 0 : overlayA_lo = mul16bit_8bit_result(overlayA_lo, opacity);
696 0 : overlayA_hi = mul16bit_8bit_result(overlayA_hi, opacity);
697 : auto srcAMul255MinusOverlayA_lo =
698 0 : mul16bit_8bit_result(srcA_lo, _mm_sub_epi16(r255, overlayA_lo));
699 : auto srcAMul255MinusOverlayA_hi =
700 0 : mul16bit_8bit_result(srcA_hi, _mm_sub_epi16(r255, overlayA_hi));
701 0 : auto dstA_lo = _mm_add_epi16(overlayA_lo, srcAMul255MinusOverlayA_lo);
702 0 : auto dstA_hi = _mm_add_epi16(overlayA_hi, srcAMul255MinusOverlayA_hi);
703 :
704 : // The & 0xff should not be necessary. This is mostly a safety
705 : // belt if the above math yields a result outside [0, 255]...
706 0 : dstA_lo = _mm_and_si128(dstA_lo, r255);
707 0 : dstA_hi = _mm_and_si128(dstA_hi, r255);
708 :
709 : // This would be the equivalent of a "_mm_i16gather_epi16" operation
710 : // which does not exist...
711 : // invDstA_{i} = [tabInvDstASigned[dstA_{i}] for i in range(8)]
712 0 : auto invDstA_lo = _mm_undefined_si128();
713 0 : auto invDstA_hi = _mm_undefined_si128();
714 : #define SET_INVDSTA(k) \
715 : do \
716 : { \
717 : const int idxLo = _mm_extract_epi16(dstA_lo, k); \
718 : const int idxHi = _mm_extract_epi16(dstA_hi, k); \
719 : invDstA_lo = _mm_insert_epi16(invDstA_lo, tabInvDstASigned[idxLo], k); \
720 : invDstA_hi = _mm_insert_epi16(invDstA_hi, tabInvDstASigned[idxHi], k); \
721 : } while (0)
722 :
723 0 : SET_INVDSTA(0);
724 0 : SET_INVDSTA(1);
725 0 : SET_INVDSTA(2);
726 0 : SET_INVDSTA(3);
727 0 : SET_INVDSTA(4);
728 0 : SET_INVDSTA(5);
729 0 : SET_INVDSTA(6);
730 0 : SET_INVDSTA(7);
731 :
732 0 : pack_and_store(pabyDst + i + 3 * nBandSpace, dstA_lo, dstA_hi);
733 :
734 : #define PROCESS_COMPONENT(pabySrc, pabyOverlay, iBand) \
735 : do \
736 : { \
737 : auto [src_lo, src_hi] = load_and_unpack((pabySrc) + i); \
738 : auto [overlay_lo, overlay_hi] = load_and_unpack((pabyOverlay) + i); \
739 : auto dst_lo = _mm_srli_epi16( \
740 : _mm_add_epi16( \
741 : _mm_add_epi16( \
742 : _mm_mullo_epi16(overlay_lo, overlayA_lo), \
743 : _mm_mullo_epi16(src_lo, srcAMul255MinusOverlayA_lo)), \
744 : r255), \
745 : 8); \
746 : auto dst_hi = _mm_srli_epi16( \
747 : _mm_add_epi16( \
748 : _mm_add_epi16( \
749 : _mm_mullo_epi16(overlay_hi, overlayA_hi), \
750 : _mm_mullo_epi16(src_hi, srcAMul255MinusOverlayA_hi)), \
751 : r255), \
752 : 8); \
753 : dst_lo = mul16bit_8bit_result(dst_lo, invDstA_lo); \
754 : dst_hi = mul16bit_8bit_result(dst_hi, invDstA_hi); \
755 : pack_and_store(pabyDst + i + (iBand)*nBandSpace, dst_lo, dst_hi); \
756 : } while (0)
757 :
758 0 : PROCESS_COMPONENT(pabyR, pabyOverlayR, 0);
759 0 : PROCESS_COMPONENT(pabyG, pabyOverlayG, 1);
760 0 : PROCESS_COMPONENT(pabyB, pabyOverlayB, 2);
761 : }
762 1 : return i;
763 : }
764 : #endif
765 :
766 : template <bool bPixelSpaceIsOne>
767 : #if defined(__GNUC__)
768 : __attribute__((noinline))
769 : #endif
770 : static void
771 1 : BlendSrcOverRGBA_Generic(const GByte *CPL_RESTRICT pabyR,
772 : const GByte *CPL_RESTRICT pabyG,
773 : const GByte *CPL_RESTRICT pabyB,
774 : const GByte *CPL_RESTRICT pabyA,
775 : const GByte *CPL_RESTRICT pabyOverlayR,
776 : const GByte *CPL_RESTRICT pabyOverlayG,
777 : const GByte *CPL_RESTRICT pabyOverlayB,
778 : const GByte *CPL_RESTRICT pabyOverlayA,
779 : GByte *CPL_RESTRICT pabyDst,
780 : [[maybe_unused]] GSpacing nPixelSpace,
781 : GSpacing nBandSpace, int i, int N, GByte nOpacity)
782 : {
783 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
784 : if constexpr (!bPixelSpaceIsOne)
785 : {
786 : assert(nPixelSpace != 1);
787 : }
788 : #endif
789 1 : [[maybe_unused]] GSpacing nDstOffset = 0;
790 : #if defined(HAVE_SSE2) && defined(__clang__) && !defined(__INTEL_CLANG_COMPILER)
791 : // No need to vectorize for SSE2 and clang
792 : #pragma clang loop vectorize(disable)
793 : #endif
794 2 : for (; i < N; ++i)
795 : {
796 1 : const GByte nOverlayR = pabyOverlayR[i];
797 1 : const GByte nOverlayG = pabyOverlayG[i];
798 1 : const GByte nOverlayB = pabyOverlayB[i];
799 1 : const GByte nOverlayA =
800 1 : static_cast<GByte>((pabyOverlayA[i] * nOpacity + 255) / 256);
801 1 : const GByte nR = pabyR[i];
802 1 : const GByte nG = pabyG[i];
803 1 : const GByte nB = pabyB[i];
804 1 : const GByte nA = pabyA[i];
805 1 : const GByte nSrcAMul255MinusOverlayA =
806 1 : static_cast<GByte>((nA * (255 - nOverlayA) + 255) / 256);
807 1 : const GByte nDstA =
808 : static_cast<GByte>(nOverlayA + nSrcAMul255MinusOverlayA);
809 1 : GByte nDstR = static_cast<GByte>(
810 1 : (nOverlayR * nOverlayA + nR * nSrcAMul255MinusOverlayA + 255) /
811 : 256);
812 1 : GByte nDstG = static_cast<GByte>(
813 1 : (nOverlayG * nOverlayA + nG * nSrcAMul255MinusOverlayA + 255) /
814 : 256);
815 1 : GByte nDstB = static_cast<GByte>(
816 1 : (nOverlayB * nOverlayA + nB * nSrcAMul255MinusOverlayA + 255) /
817 : 256);
818 : // nInvDstA = (255 << SHIFT_DIV_DSTA) / nDstA;
819 1 : const uint16_t nInvDstA = gTabInvDstA[nDstA];
820 1 : constexpr unsigned ROUND_OFFSET_DIV_DSTA = ((1 << SHIFT_DIV_DSTA) - 1);
821 1 : nDstR = static_cast<GByte>((nDstR * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
822 : SHIFT_DIV_DSTA);
823 1 : nDstG = static_cast<GByte>((nDstG * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
824 : SHIFT_DIV_DSTA);
825 1 : nDstB = static_cast<GByte>((nDstB * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
826 : SHIFT_DIV_DSTA);
827 : if constexpr (bPixelSpaceIsOne)
828 : {
829 : pabyDst[i] = nDstR;
830 : pabyDst[i + nBandSpace] = nDstG;
831 : pabyDst[i + 2 * nBandSpace] = nDstB;
832 : pabyDst[i + 3 * nBandSpace] = nDstA;
833 : }
834 : else
835 : {
836 1 : pabyDst[nDstOffset] = nDstR;
837 1 : pabyDst[nDstOffset + nBandSpace] = nDstG;
838 1 : pabyDst[nDstOffset + 2 * nBandSpace] = nDstB;
839 1 : pabyDst[nDstOffset + 3 * nBandSpace] = nDstA;
840 1 : nDstOffset += nPixelSpace;
841 : }
842 : }
843 1 : }
844 :
845 : /************************************************************************/
846 : /* BlendDataset::IRasterIO() */
847 : /************************************************************************/
848 :
849 229 : CPLErr BlendDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
850 : int nXSize, int nYSize, void *pData,
851 : int nBufXSize, int nBufYSize,
852 : GDALDataType eBufType, int nBandCount,
853 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
854 : GSpacing nLineSpace, GSpacing nBandSpace,
855 : GDALRasterIOExtraArg *psExtraArg)
856 : {
857 : // Try to pass the request to the most appropriate overview dataset.
858 229 : if (nBufXSize < nXSize && nBufYSize < nYSize)
859 : {
860 2 : int bTried = FALSE;
861 2 : const CPLErr eErr = TryOverviewRasterIO(
862 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
863 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
864 : nBandSpace, psExtraArg, &bTried);
865 2 : if (bTried)
866 2 : return eErr;
867 : }
868 :
869 227 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
870 227 : const int nColorCount = m_oColorDS.GetRasterCount();
871 227 : const int nOverlayCount = m_oOverlayDS.GetRasterCount();
872 419 : if (nOverlayCount == 1 && m_opacity255Scale == 255 &&
873 410 : m_operator == HSV_VALUE && eRWFlag == GF_Read && eBufType == GDT_Byte &&
874 643 : nBandCount == nBands && IsAllBands(nBands, panBandMap) &&
875 203 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
876 : psExtraArg))
877 : {
878 197 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
879 197 : const GByte *pabyR = m_abyBuffer.data();
880 197 : const GByte *pabyG = m_abyBuffer.data() + nPixelCount;
881 197 : const GByte *pabyB = m_abyBuffer.data() + nPixelCount * 2;
882 197 : const GByte *pabyValue = m_abyBuffer.data() + nPixelCount * nColorCount;
883 197 : size_t nSrcIdx = 0;
884 516 : for (int j = 0; j < nBufYSize; ++j)
885 : {
886 319 : auto nDstOffset = j * nLineSpace;
887 319 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize &&
888 317 : nBandSpace >= nLineSpace * nBufYSize)
889 : {
890 317 : patch_value_line(nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
891 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
892 317 : pabyDst + nDstOffset,
893 317 : pabyDst + nDstOffset + nBandSpace,
894 317 : pabyDst + nDstOffset + 2 * nBandSpace);
895 317 : nSrcIdx += nBufXSize;
896 : }
897 : else
898 : {
899 262146 : for (int i = 0; i < nBufXSize;
900 262144 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
901 : {
902 : float h, s;
903 262144 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
904 : &h, &s);
905 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
906 262144 : &pabyDst[nDstOffset + 0 * nBandSpace],
907 262144 : &pabyDst[nDstOffset + 1 * nBandSpace],
908 262144 : &pabyDst[nDstOffset + 2 * nBandSpace]);
909 : }
910 : }
911 : }
912 197 : if (nColorCount == 4)
913 : {
914 394 : for (int j = 0; j < nBufYSize; ++j)
915 : {
916 198 : auto nDstOffset = 3 * nBandSpace + j * nLineSpace;
917 198 : const GByte *pabyA = m_abyBuffer.data() + nPixelCount * 3;
918 198 : GDALCopyWords64(pabyA, GDT_Byte, 1, pabyDst + nDstOffset,
919 : GDT_Byte, static_cast<int>(nPixelSpace),
920 : nBufXSize);
921 : }
922 : }
923 :
924 197 : return CE_None;
925 : }
926 3 : else if (nOverlayCount == 4 && nColorCount == 4 && m_operator == SRC_OVER &&
927 1 : eRWFlag == GF_Read && eBufType == GDT_Byte &&
928 34 : nBandCount == nBands && IsAllBands(nBands, panBandMap) &&
929 1 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize,
930 : nBufYSize, psExtraArg))
931 : {
932 1 : const GByte nOpacity = static_cast<GByte>(m_opacity255Scale);
933 1 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
934 1 : const GByte *CPL_RESTRICT pabyR = m_abyBuffer.data();
935 1 : const GByte *CPL_RESTRICT pabyG = m_abyBuffer.data() + nPixelCount;
936 1 : const GByte *CPL_RESTRICT pabyB = m_abyBuffer.data() + nPixelCount * 2;
937 1 : const GByte *CPL_RESTRICT pabyA = m_abyBuffer.data() + nPixelCount * 3;
938 : const GByte *CPL_RESTRICT pabyOverlayR =
939 1 : m_abyBuffer.data() + nPixelCount * nColorCount;
940 : const GByte *CPL_RESTRICT pabyOverlayG =
941 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
942 : const GByte *CPL_RESTRICT pabyOverlayB =
943 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 2);
944 : const GByte *CPL_RESTRICT pabyOverlayA =
945 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 3);
946 1 : size_t nSrcIdx = 0;
947 2 : for (int j = 0; j < nBufYSize; ++j)
948 : {
949 1 : auto nDstOffset = j * nLineSpace;
950 1 : int i = 0;
951 : #ifdef HAVE_SSE2
952 1 : if (nPixelSpace == 1)
953 : {
954 2 : i = BlendSrcOverRGBA_SSE2(
955 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
956 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
957 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
958 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nBandSpace,
959 : nBufXSize, nOpacity);
960 1 : nSrcIdx += i;
961 1 : nDstOffset += i;
962 : }
963 : #endif
964 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
965 : if (nPixelSpace == 1)
966 : {
967 : // Note: clang 20 does a rather good job at autovectorizing
968 : // for SSE2, but BlendSrcOverRGBA_SSE2() performs better.
969 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ true>(
970 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
971 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
972 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
973 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset,
974 : /*nPixelSpace = */ 1, nBandSpace, i, nBufXSize, nOpacity);
975 : }
976 : else
977 : #endif
978 : {
979 1 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ false>(
980 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
981 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
982 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
983 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nPixelSpace,
984 : nBandSpace, i, nBufXSize, nOpacity);
985 : }
986 1 : nSrcIdx += nBufXSize - i;
987 : }
988 1 : return CE_None;
989 : }
990 29 : else if (m_ioError)
991 : {
992 6 : return CE_Failure;
993 : }
994 : else
995 : {
996 23 : const CPLErr eErr = GDALDataset::IRasterIO(
997 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
998 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
999 : nBandSpace, psExtraArg);
1000 23 : m_ioError = eErr != CE_None;
1001 23 : return eErr;
1002 : }
1003 : }
1004 :
1005 : /************************************************************************/
1006 : /* SrcOverRGBOneComponent() */
1007 : /************************************************************************/
1008 :
1009 : // GCC and clang do a god job a auto vectorizing the below function
1010 : #if defined(__GNUC__) && !defined(__clang__)
1011 : __attribute__((optimize("tree-vectorize")))
1012 : #endif
1013 : static void
1014 12 : SrcOverRGB(const uint8_t *const __restrict pabyOverlay,
1015 : const uint8_t *const __restrict pabySrc,
1016 : uint8_t *const __restrict pabyDst, const size_t N,
1017 : const uint8_t nOpacity)
1018 : {
1019 915 : for (size_t i = 0; i < N; ++i)
1020 : {
1021 903 : const uint8_t nOverlay = pabyOverlay[i];
1022 903 : const uint8_t nSrc = pabySrc[i];
1023 903 : pabyDst[i] = static_cast<uint8_t>(
1024 903 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) / 256);
1025 : }
1026 12 : }
1027 :
1028 : /************************************************************************/
1029 : /* BlendBand::IRasterIO() */
1030 : /************************************************************************/
1031 :
1032 212 : CPLErr BlendBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
1033 : int nXSize, int nYSize, void *pData, int nBufXSize,
1034 : int nBufYSize, GDALDataType eBufType,
1035 : GSpacing nPixelSpace, GSpacing nLineSpace,
1036 : GDALRasterIOExtraArg *psExtraArg)
1037 : {
1038 : // Try to pass the request to the most appropriate overview dataset.
1039 212 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1040 : {
1041 1 : int bTried = FALSE;
1042 1 : const CPLErr eErr = TryOverviewRasterIO(
1043 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1044 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
1045 1 : if (bTried)
1046 1 : return eErr;
1047 : }
1048 :
1049 211 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1050 211 : const int nColorCount = m_oBlendDataset.m_oColorDS.GetRasterCount();
1051 211 : const int nOverlayCount = m_oBlendDataset.m_oOverlayDS.GetRasterCount();
1052 211 : if (nBand == 4 && m_oBlendDataset.m_operator == HSV_VALUE)
1053 : {
1054 9 : if (nColorCount == 3)
1055 : {
1056 0 : GByte ch = 255;
1057 0 : for (int iY = 0; iY < nBufYSize; ++iY)
1058 : {
1059 0 : GDALCopyWords64(&ch, GDT_Byte, 0,
1060 0 : static_cast<GByte *>(pData) + iY * nLineSpace,
1061 : eBufType, static_cast<int>(nPixelSpace),
1062 : nBufXSize);
1063 : }
1064 0 : return CE_None;
1065 : }
1066 : else
1067 : {
1068 9 : CPLAssert(nColorCount == 4);
1069 9 : return m_oBlendDataset.m_oColorDS.GetRasterBand(4)->RasterIO(
1070 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
1071 9 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg);
1072 : }
1073 : }
1074 43 : else if (nOverlayCount == 3 && nColorCount == 3 &&
1075 30 : m_oBlendDataset.m_operator == SRC_OVER && eRWFlag == GF_Read &&
1076 227 : eBufType == GDT_Byte &&
1077 12 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
1078 : nBufXSize, nBufYSize,
1079 : psExtraArg))
1080 : {
1081 12 : const int nOpacity = m_oBlendDataset.m_opacity255Scale;
1082 : const GByte *const CPL_RESTRICT pabySrc =
1083 12 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * (nBand - 1);
1084 : const GByte *const CPL_RESTRICT pabyOverlay =
1085 12 : m_oBlendDataset.m_abyBuffer.data() +
1086 12 : nPixelCount * (nColorCount + nBand - 1);
1087 12 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
1088 12 : size_t nSrcIdx = 0;
1089 24 : for (int j = 0; j < nBufYSize; ++j)
1090 : {
1091 12 : auto nDstOffset = j * nLineSpace;
1092 12 : if (nPixelSpace == 1)
1093 : {
1094 12 : SrcOverRGB(pabyOverlay + nSrcIdx, pabySrc + nSrcIdx,
1095 12 : pabyDst + nDstOffset, nBufXSize,
1096 : static_cast<uint8_t>(nOpacity));
1097 12 : nSrcIdx += nBufXSize;
1098 : }
1099 : else
1100 : {
1101 0 : for (int i = 0; i < nBufXSize;
1102 0 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1103 : {
1104 0 : const int nOverlay = pabyOverlay[nSrcIdx];
1105 0 : const int nSrc = pabySrc[nSrcIdx];
1106 0 : pabyDst[nDstOffset] = static_cast<GByte>(
1107 0 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) /
1108 : 256);
1109 : }
1110 : }
1111 : }
1112 12 : return CE_None;
1113 : }
1114 321 : else if (eRWFlag == GF_Read && eBufType == GDT_Byte &&
1115 131 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
1116 : nBufXSize, nBufYSize,
1117 : psExtraArg))
1118 : {
1119 131 : GByte *pabyDst = static_cast<GByte *>(pData);
1120 131 : if (m_oBlendDataset.m_operator == SRC_OVER)
1121 : {
1122 3 : const auto RGBToGrayScale = [](int R, int G, int B)
1123 : {
1124 : // Equivalent to R * 0.299 + G * 0.587 + B * 0.114
1125 : // but using faster computation
1126 3 : return (R * 306 + G * 601 + B * 117) / 1024;
1127 : };
1128 :
1129 : const GByte *paby =
1130 83 : (nBand <= nColorCount) ? m_oBlendDataset.m_abyBuffer.data() +
1131 83 : nPixelCount * (nBand - 1)
1132 0 : : (nBand == 4 && nColorCount == 2)
1133 0 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
1134 83 : : nullptr;
1135 : const GByte *pabyA =
1136 : (nColorCount == 4)
1137 87 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 3
1138 : : nColorCount == 2
1139 4 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
1140 83 : : nullptr;
1141 : const GByte *pabyOverlay =
1142 83 : (nBand <= nOverlayCount)
1143 91 : ? m_oBlendDataset.m_abyBuffer.data() +
1144 75 : nPixelCount * (nColorCount + nBand - 1)
1145 8 : : (nBand <= 3) ? m_oBlendDataset.m_abyBuffer.data() +
1146 7 : nPixelCount * nColorCount
1147 83 : : nullptr;
1148 : const GByte *pabyOverlayA =
1149 80 : (nOverlayCount == 2 || nOverlayCount == 4)
1150 163 : ? m_oBlendDataset.m_abyBuffer.data() +
1151 62 : nPixelCount * (nColorCount + nOverlayCount - 1)
1152 83 : : nullptr;
1153 : const GByte *pabyOverlayR =
1154 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1155 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1156 3 : nPixelCount * nColorCount
1157 83 : : nullptr;
1158 : const GByte *pabyOverlayG =
1159 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1160 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1161 3 : nPixelCount * (nColorCount + 1)
1162 83 : : nullptr;
1163 : const GByte *pabyOverlayB =
1164 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
1165 149 : ? m_oBlendDataset.m_abyBuffer.data() +
1166 3 : nPixelCount * (nColorCount + 2)
1167 83 : : nullptr;
1168 :
1169 83 : size_t nSrcIdx = 0;
1170 714 : for (int j = 0; j < nBufYSize; ++j)
1171 : {
1172 631 : auto nDstOffset = j * nLineSpace;
1173 98753 : for (int i = 0; i < nBufXSize;
1174 98122 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1175 : {
1176 : // Corrected to take into account m_opacity255Scale
1177 98122 : const int nOverlayA =
1178 98122 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
1179 97606 : m_oBlendDataset.m_opacity255Scale +
1180 : 255) /
1181 : 256)
1182 516 : : m_oBlendDataset.m_opacity255Scale;
1183 :
1184 98122 : const int nSrcA = pabyA ? pabyA[nSrcIdx] : 255;
1185 :
1186 98122 : const int nSrcAMul255MinusOverlayA =
1187 98122 : (nSrcA * (255 - nOverlayA) + 255) / 256;
1188 98122 : const int nDstA = nOverlayA + nSrcAMul255MinusOverlayA;
1189 98122 : if (nBand != 4)
1190 : {
1191 : const int nOverlay =
1192 3 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
1193 147439 : ? RGBToGrayScale(pabyOverlayR[nSrcIdx],
1194 3 : pabyOverlayG[nSrcIdx],
1195 3 : pabyOverlayB[nSrcIdx])
1196 73718 : : pabyOverlay ? pabyOverlay[nSrcIdx]
1197 73721 : : 255;
1198 :
1199 73721 : const int nSrc = paby ? paby[nSrcIdx] : 255;
1200 73721 : int nDst = (nOverlay * nOverlayA +
1201 73721 : nSrc * nSrcAMul255MinusOverlayA + 255) /
1202 : 256;
1203 73721 : if (nDstA != 0 && nDstA != 255)
1204 19672 : nDst = (nDst * 255 + (nDstA / 2)) / nDstA;
1205 73721 : pabyDst[nDstOffset] =
1206 73721 : static_cast<GByte>(std::min(nDst, 255));
1207 : }
1208 : else
1209 : {
1210 24401 : pabyDst[nDstOffset] = static_cast<GByte>(nDstA);
1211 : }
1212 : }
1213 : }
1214 : }
1215 48 : else if (nOverlayCount == 1 && m_oBlendDataset.m_opacity255Scale == 255)
1216 : {
1217 18 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
1218 : const GByte *pabyG =
1219 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
1220 : const GByte *pabyB =
1221 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
1222 18 : CPLAssert(m_oBlendDataset.m_operator == HSV_VALUE);
1223 18 : size_t nSrcIdx = 0;
1224 : const GByte *pabyValue =
1225 18 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
1226 402 : for (int j = 0; j < nBufYSize; ++j)
1227 : {
1228 384 : auto nDstOffset = j * nLineSpace;
1229 384 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize)
1230 : {
1231 884 : patch_value_line(
1232 : nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
1233 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
1234 378 : nBand == 1 ? pabyDst + nDstOffset : nullptr,
1235 378 : nBand == 2 ? pabyDst + nDstOffset : nullptr,
1236 378 : nBand == 3 ? pabyDst + nDstOffset : nullptr);
1237 378 : nSrcIdx += nBufXSize;
1238 : }
1239 : else
1240 : {
1241 786438 : for (int i = 0; i < nBufXSize;
1242 786432 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1243 : {
1244 : float h, s;
1245 786432 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx],
1246 786432 : pabyB[nSrcIdx], &h, &s);
1247 786432 : if (nBand == 1)
1248 : {
1249 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
1250 262144 : &pabyDst[nDstOffset], nullptr, nullptr);
1251 : }
1252 524288 : else if (nBand == 2)
1253 : {
1254 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
1255 262144 : &pabyDst[nDstOffset], nullptr);
1256 : }
1257 : else
1258 : {
1259 262144 : CPLAssert(nBand == 3);
1260 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
1261 262144 : nullptr, &pabyDst[nDstOffset]);
1262 : }
1263 : }
1264 : }
1265 18 : }
1266 : }
1267 : else
1268 : {
1269 30 : CPLAssert(m_oBlendDataset.m_operator == HSV_VALUE);
1270 30 : CPLAssert(nBand <= 3);
1271 30 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
1272 : const GByte *pabyG =
1273 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
1274 : const GByte *pabyB =
1275 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
1276 : const GByte *pabyValue =
1277 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
1278 : const GByte *pabyOverlayR =
1279 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1280 12 : nPixelCount * nColorCount
1281 30 : : nullptr;
1282 : const GByte *pabyOverlayG =
1283 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1284 12 : nPixelCount * (nColorCount + 1)
1285 30 : : nullptr;
1286 : const GByte *pabyOverlayB =
1287 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
1288 12 : nPixelCount * (nColorCount + 2)
1289 30 : : nullptr;
1290 : const GByte *pabyOverlayA =
1291 24 : (nOverlayCount == 2 || nOverlayCount == 4)
1292 54 : ? m_oBlendDataset.m_abyBuffer.data() +
1293 12 : nPixelCount * (nColorCount + nOverlayCount - 1)
1294 30 : : nullptr;
1295 :
1296 30 : size_t nSrcIdx = 0;
1297 60 : for (int j = 0; j < nBufYSize; ++j)
1298 : {
1299 30 : auto nDstOffset = j * nLineSpace;
1300 3932190 : for (int i = 0; i < nBufXSize;
1301 3932160 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1302 : {
1303 3932160 : const int nColorR = pabyR[nSrcIdx];
1304 3932160 : const int nColorG = pabyG[nSrcIdx];
1305 3932160 : const int nColorB = pabyB[nSrcIdx];
1306 : const int nOverlayV =
1307 1572860 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
1308 5505020 : ? std::max({pabyOverlayR[nSrcIdx],
1309 1572860 : pabyOverlayG[nSrcIdx],
1310 1572860 : pabyOverlayB[nSrcIdx]})
1311 2359300 : : pabyValue[nSrcIdx];
1312 3932160 : const int nOverlayA =
1313 3932160 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
1314 1572860 : m_oBlendDataset.m_opacity255Scale +
1315 : 255) /
1316 : 256)
1317 2359300 : : m_oBlendDataset.m_opacity255Scale;
1318 : const int nColorValue =
1319 3932160 : std::max({nColorR, nColorG, nColorB});
1320 :
1321 : float h, s;
1322 3932160 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
1323 : &h, &s);
1324 :
1325 3932160 : const GByte nTargetValue = static_cast<GByte>(
1326 3932160 : (nOverlayV * nOverlayA +
1327 3932160 : nColorValue * (255 - nOverlayA) + 255) /
1328 : 256);
1329 :
1330 3932160 : if (nBand == 1)
1331 : {
1332 1310720 : hsv_to_rgb(h, s, nTargetValue, &pabyDst[nDstOffset],
1333 : nullptr, nullptr);
1334 : }
1335 2621440 : else if (nBand == 2)
1336 : {
1337 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr,
1338 1310720 : &pabyDst[nDstOffset], nullptr);
1339 : }
1340 : else
1341 : {
1342 1310720 : CPLAssert(nBand == 3);
1343 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr, nullptr,
1344 1310720 : &pabyDst[nDstOffset]);
1345 : }
1346 : }
1347 : }
1348 : }
1349 :
1350 131 : return CE_None;
1351 : }
1352 59 : else if (m_oBlendDataset.m_ioError)
1353 : {
1354 0 : return CE_Failure;
1355 : }
1356 : else
1357 : {
1358 59 : const CPLErr eErr = GDALRasterBand::IRasterIO(
1359 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1360 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
1361 59 : m_oBlendDataset.m_ioError = eErr != CE_None;
1362 59 : return eErr;
1363 : }
1364 : }
1365 :
1366 : } // namespace
1367 :
1368 : /************************************************************************/
1369 : /* GDALRasterBlendAlgorithm::ValidateGlobal() */
1370 : /************************************************************************/
1371 :
1372 80 : bool GDALRasterBlendAlgorithm::ValidateGlobal()
1373 : {
1374 : auto poSrcDS =
1375 80 : m_inputDataset.empty() ? nullptr : m_inputDataset[0].GetDatasetRef();
1376 80 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
1377 80 : if (poSrcDS)
1378 : {
1379 154 : if (poSrcDS->GetRasterCount() == 0 || poSrcDS->GetRasterCount() > 4 ||
1380 77 : poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte)
1381 : {
1382 1 : ReportError(CE_Failure, CPLE_NotSupported,
1383 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
1384 : "supported as input");
1385 1 : return false;
1386 : }
1387 : }
1388 79 : if (poOverlayDS)
1389 : {
1390 79 : if (poOverlayDS->GetRasterCount() == 0 ||
1391 158 : poOverlayDS->GetRasterCount() > 4 ||
1392 79 : poOverlayDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte)
1393 : {
1394 1 : ReportError(CE_Failure, CPLE_NotSupported,
1395 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
1396 : "supported as overlay");
1397 1 : return false;
1398 : }
1399 : }
1400 :
1401 78 : if (poSrcDS && poOverlayDS)
1402 : {
1403 149 : if (poSrcDS->GetRasterXSize() != poOverlayDS->GetRasterXSize() ||
1404 74 : poSrcDS->GetRasterYSize() != poOverlayDS->GetRasterYSize())
1405 : {
1406 2 : ReportError(CE_Failure, CPLE_IllegalArg,
1407 : "Input dataset and overlay dataset must have "
1408 : "the same dimensions");
1409 2 : return false;
1410 : }
1411 :
1412 84 : if (m_operator == HSV_VALUE && poSrcDS->GetRasterCount() != 3 &&
1413 11 : poSrcDS->GetRasterCount() != 4)
1414 : {
1415 1 : ReportError(CE_Failure, CPLE_AppDefined,
1416 : "Operator %s requires a 3-band or 4-band input dataset",
1417 : HSV_VALUE);
1418 1 : return false;
1419 : }
1420 : }
1421 :
1422 75 : return true;
1423 : }
1424 :
1425 : /************************************************************************/
1426 : /* GDALRasterBlendAlgorithm::RunStep() */
1427 : /************************************************************************/
1428 :
1429 36 : bool GDALRasterBlendAlgorithm::RunStep(GDALPipelineStepRunContext &)
1430 : {
1431 36 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
1432 36 : CPLAssert(poSrcDS);
1433 :
1434 36 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
1435 36 : CPLAssert(poOverlayDS);
1436 :
1437 36 : if (!ValidateGlobal())
1438 0 : return false;
1439 :
1440 36 : const int nOpacity255Scale =
1441 36 : (m_opacity * 255 + OPACITY_INPUT_RANGE / 2) / OPACITY_INPUT_RANGE;
1442 :
1443 36 : m_outputDataset.Set(std::make_unique<BlendDataset>(
1444 36 : *poSrcDS, *poOverlayDS, m_operator, nOpacity255Scale));
1445 :
1446 36 : return true;
1447 : }
1448 :
1449 : GDALRasterBlendAlgorithmStandalone::~GDALRasterBlendAlgorithmStandalone() =
1450 : default;
1451 :
1452 : //! @endcond
|
