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 : #include "gdal_utils.h"
19 :
20 : #include <algorithm>
21 : #include <array>
22 : #include <cassert>
23 : #include <limits>
24 : #include <utility>
25 :
26 : #if defined(__x86_64) || defined(_M_X64)
27 : #define HAVE_SSE2
28 : #include <immintrin.h>
29 : #endif
30 : #ifdef HAVE_SSE2
31 : #include "gdalsse_priv.h"
32 : #endif
33 :
34 : //! @cond Doxygen_Suppress
35 :
36 : #ifndef _
37 : #define _(x) (x)
38 : #endif
39 :
40 : /************************************************************************/
41 : /* CompositionModes */
42 : /************************************************************************/
43 614 : std::map<CompositionMode, std::string> CompositionModes()
44 : {
45 : return {
46 1228 : {CompositionMode::SRC_OVER, "src-over"},
47 0 : {CompositionMode::HSV_VALUE, "hsv-value"},
48 0 : {CompositionMode::MULTIPLY, "multiply"},
49 0 : {CompositionMode::SCREEN, "screen"},
50 0 : {CompositionMode::OVERLAY, "overlay"},
51 0 : {CompositionMode::HARD_LIGHT, "hard-light"},
52 0 : {CompositionMode::DARKEN, "darken"},
53 0 : {CompositionMode::LIGHTEN, "lighten"},
54 0 : {CompositionMode::COLOR_BURN, "color-burn"},
55 0 : {CompositionMode::COLOR_DODGE, "color-dodge"},
56 6754 : };
57 : }
58 :
59 : /************************************************************************/
60 : /* CompositionModeToString() */
61 : /************************************************************************/
62 :
63 229 : std::string CompositionModeToString(CompositionMode mode)
64 : {
65 458 : const auto &modes = CompositionModes();
66 229 : const auto &iter = modes.find(mode);
67 229 : if (iter != modes.end())
68 : {
69 229 : return iter->second;
70 : }
71 0 : CPLError(CE_Failure, CPLE_IllegalArg,
72 : "Invalid composition mode value: %d, returning 'src-over'",
73 : static_cast<int>(mode));
74 0 : return "src-over";
75 : }
76 :
77 : /************************************************************************/
78 : /* CompositionModesIdentifiers() */
79 : /************************************************************************/
80 :
81 228 : std::vector<std::string> CompositionModesIdentifiers()
82 : {
83 456 : const auto &modes = CompositionModes();
84 228 : std::vector<std::string> identifiers;
85 2508 : for (const auto &pair : modes)
86 : {
87 2280 : identifiers.push_back(pair.second);
88 : }
89 456 : return identifiers;
90 : }
91 :
92 : /************************************************************************/
93 : /* CompositionModeFromString() */
94 : /************************************************************************/
95 :
96 157 : CompositionMode CompositionModeFromString(const std::string &str)
97 : {
98 314 : const auto &modes = CompositionModes();
99 : auto iter =
100 : std::find_if(modes.begin(), modes.end(),
101 965 : [&str](const auto &pair) { return pair.second == str; });
102 157 : if (iter != modes.end())
103 : {
104 157 : return iter->first;
105 : }
106 0 : CPLError(CE_Failure, CPLE_IllegalArg,
107 : "Invalid composition identifier: %s, returning SRC_OVER",
108 : str.c_str());
109 0 : return CompositionMode::SRC_OVER;
110 : }
111 :
112 : /************************************************************************/
113 : /* MinBandCountForCompositionMode() */
114 : /************************************************************************/
115 :
116 : //! Returns the minimum number of bands required for the given composition mode
117 334 : int MinBandCountForCompositionMode(CompositionMode mode)
118 : {
119 334 : switch (mode)
120 : {
121 36 : case CompositionMode::HSV_VALUE:
122 36 : return 3;
123 298 : case CompositionMode::SRC_OVER:
124 : case CompositionMode::MULTIPLY:
125 : case CompositionMode::SCREEN:
126 : case CompositionMode::OVERLAY:
127 : case CompositionMode::HARD_LIGHT:
128 : case CompositionMode::DARKEN:
129 : case CompositionMode::LIGHTEN:
130 : case CompositionMode::COLOR_BURN:
131 : case CompositionMode::COLOR_DODGE:
132 298 : return 1;
133 : }
134 : // unreachable...
135 0 : return 1;
136 : }
137 :
138 : /************************************************************************/
139 : /* MaxBandCountForCompositionMode() */
140 : /************************************************************************/
141 :
142 : /**
143 : * Returns the maximum number of bands allowed for the given composition mode
144 : */
145 334 : int MaxBandCountForCompositionMode(CompositionMode mode)
146 : {
147 334 : switch (mode)
148 : {
149 334 : case CompositionMode::SRC_OVER:
150 : case CompositionMode::HSV_VALUE:
151 : case CompositionMode::MULTIPLY:
152 : case CompositionMode::SCREEN:
153 : case CompositionMode::OVERLAY:
154 : case CompositionMode::HARD_LIGHT:
155 : case CompositionMode::DARKEN:
156 : case CompositionMode::LIGHTEN:
157 : case CompositionMode::COLOR_BURN:
158 : case CompositionMode::COLOR_DODGE:
159 334 : return 4;
160 : }
161 : // unreachable...
162 0 : return 4;
163 : }
164 :
165 : /************************************************************************/
166 : /* BandCountIsCompatibleWithCompositionMode() */
167 : /************************************************************************/
168 :
169 : //! Checks whether the number of bands is compatible with the given composition mode
170 333 : bool BandCountIsCompatibleWithCompositionMode(int bandCount,
171 : CompositionMode mode)
172 : {
173 333 : const int minBands = MinBandCountForCompositionMode(mode);
174 333 : const int maxBands = MaxBandCountForCompositionMode(mode);
175 333 : return minBands <= bandCount && bandCount <= maxBands;
176 : }
177 :
178 : /************************************************************************/
179 : /* MulScale255() */
180 : /************************************************************************/
181 :
182 : /** Multiply 2 bytes considering them as ratios with 255 = 100%, and return their product unscaled to [0, 255], by ceiling */
183 10964 : inline GByte MulScale255(GByte a, GByte b)
184 : {
185 10964 : return static_cast<GByte>((a * b + 255) / 256);
186 : }
187 :
188 : /************************************************************************/
189 : /* ProcessAlphaChannels */
190 : /************************************************************************/
191 :
192 824 : static inline void ProcessAlphaChannels(size_t i,
193 : const GByte *CPL_RESTRICT pabyA,
194 : const GByte *CPL_RESTRICT pabyOverlayA,
195 : int nOpacity, bool bSwappedOpacity,
196 : GByte &outA, GByte &outOverlaA,
197 : GByte &outFinalAlpha)
198 : {
199 : // Apply opacity depending on whether overlay and base were swapped
200 824 : const GByte byOpacity = static_cast<GByte>(nOpacity);
201 824 : if (!bSwappedOpacity)
202 : {
203 756 : outOverlaA =
204 756 : pabyOverlayA ? MulScale255(pabyOverlayA[i], byOpacity) : byOpacity;
205 :
206 756 : outA = pabyA ? pabyA[i] : 255;
207 : }
208 : else
209 : {
210 68 : outOverlaA = pabyOverlayA ? pabyOverlayA[i] : 255;
211 68 : if (outOverlaA != 255)
212 : {
213 0 : outOverlaA = pabyOverlayA ? pabyOverlayA[i] : 255;
214 : }
215 68 : outA = pabyA ? MulScale255(pabyA[i], byOpacity) : byOpacity;
216 : }
217 :
218 : // Da' = Sa + Da - Sa.Da
219 824 : outFinalAlpha =
220 824 : static_cast<GByte>(outOverlaA + outA - MulScale255(outOverlaA, outA));
221 824 : }
222 :
223 : /************************************************************************/
224 : /* DivScale255() */
225 : /************************************************************************/
226 :
227 : /** Divide 2 bytes considering them as ratios with 255 = 100%, and return their quotient unscaled to [0, 255], by flooring
228 : * \warning Caution: this function does not check that the result actually fits on a byte, and just casts the computed value to byte.
229 : */
230 2456 : inline GByte DivScale255(GByte a, GByte b)
231 : {
232 2456 : if (a == 0)
233 : {
234 272 : return 0;
235 : }
236 2184 : else if (b == 0)
237 : {
238 0 : return 255;
239 : }
240 : else
241 : {
242 2184 : const int nRes = (a * 255) / b;
243 : #ifdef DEBUG
244 2184 : CPLAssert(nRes <= 255);
245 : #endif
246 2184 : return static_cast<GByte>(nRes);
247 : }
248 : }
249 :
250 : /************************************************************************/
251 : /* PremultiplyChannels */
252 : /************************************************************************/
253 :
254 : //! Premultiply RGB channels by alpha (A)
255 1648 : static inline void PremultiplyChannels(size_t i, const GByte *pabyR,
256 : const GByte *pabyG, const GByte *pabyB,
257 : GByte &outR, GByte &outG, GByte &outB,
258 : const GByte &A)
259 :
260 : {
261 1648 : if (A == 255)
262 : {
263 960 : outR = pabyR ? pabyR[i] : 255;
264 960 : outG = pabyG ? pabyG[i] : outR; // in case only R is present
265 960 : outB = pabyB ? pabyB[i] : outR; // in case only R is present
266 : }
267 : else
268 : {
269 688 : outR = pabyR ? MulScale255(pabyR[i], A) : A;
270 688 : outG = pabyG ? MulScale255(pabyG[i], A)
271 : : outR; // in case only R is present
272 688 : outB = pabyB ? MulScale255(pabyB[i], A)
273 : : outR; // in case only R is present
274 : }
275 1648 : }
276 :
277 : /************************************************************************/
278 : /* GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm() */
279 : /************************************************************************/
280 :
281 228 : GDALRasterBlendAlgorithm::GDALRasterBlendAlgorithm(bool standaloneStep)
282 : : GDALRasterPipelineStepAlgorithm(
283 : NAME, DESCRIPTION, HELP_URL,
284 0 : ConstructorOptions()
285 228 : .SetStandaloneStep(standaloneStep)
286 228 : .SetAddDefaultArguments(false)
287 456 : .SetInputDatasetHelpMsg(_("Input raster dataset"))
288 456 : .SetInputDatasetAlias("color-input")
289 456 : .SetInputDatasetMetaVar("COLOR-INPUT")
290 684 : .SetOutputDatasetHelpMsg(_("Output raster dataset")))
291 : {
292 228 : const auto AddOverlayDatasetArg = [this]()
293 : {
294 : auto &arg = AddArg("overlay", 0, _("Overlay dataset"),
295 456 : &m_overlayDataset, GDAL_OF_RASTER)
296 228 : .SetPositional()
297 228 : .SetRequired();
298 :
299 228 : SetAutoCompleteFunctionForFilename(arg, GDAL_OF_RASTER);
300 228 : };
301 :
302 228 : if (standaloneStep)
303 : {
304 186 : AddRasterInputArgs(false, false);
305 186 : AddOverlayDatasetArg();
306 186 : AddProgressArg();
307 186 : AddRasterOutputArgs(false);
308 : }
309 : else
310 : {
311 42 : AddRasterHiddenInputDatasetArg();
312 42 : AddOverlayDatasetArg();
313 : }
314 :
315 : const std::vector<std::string> compositionModeChoices{
316 228 : CompositionModesIdentifiers()};
317 456 : AddArg("operator", 0, _("Composition operator"), &m_operatorIdentifier)
318 228 : .SetChoices(compositionModeChoices)
319 456 : .SetDefault(CompositionModeToString(CompositionMode::SRC_OVER))
320 : .AddAction(
321 157 : [this]()
322 385 : { m_operator = CompositionModeFromString(m_operatorIdentifier); });
323 :
324 : AddArg("opacity", 0,
325 : _("Opacity percentage to apply to the overlay dataset (0=fully "
326 : "transparent, 100=full use of overlay opacity)"),
327 456 : &m_opacity)
328 228 : .SetDefault(m_opacity)
329 228 : .SetMinValueIncluded(0)
330 228 : .SetMaxValueIncluded(OPACITY_INPUT_RANGE);
331 :
332 410 : AddValidationAction([this]() { return ValidateGlobal(); });
333 228 : }
334 :
335 : namespace
336 : {
337 :
338 : /************************************************************************/
339 : /* BlendDataset */
340 : /************************************************************************/
341 :
342 : class BlendDataset final : public GDALDataset
343 : {
344 : public:
345 : BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,
346 : const CompositionMode eOperator, int nOpacity255Scale,
347 : bool bSwappedOpacity);
348 : ~BlendDataset() override;
349 :
350 2 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
351 : {
352 2 : return m_oColorDS.GetGeoTransform(gt);
353 : }
354 :
355 2 : const OGRSpatialReference *GetSpatialRef() const override
356 : {
357 2 : return m_oColorDS.GetSpatialRef();
358 : }
359 :
360 : bool AcquireSourcePixels(int nXOff, int nYOff, int nXSize, int nYSize,
361 : int nBufXSize, int nBufYSize,
362 : GDALRasterIOExtraArg *psExtraArg);
363 :
364 : protected:
365 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
366 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
367 : GDALDataType eBufType, int nBandCount,
368 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
369 : GSpacing nLineSpace, GSpacing nBandSpace,
370 : GDALRasterIOExtraArg *psExtraArg) override;
371 :
372 : private:
373 : friend class BlendBand;
374 : GDALDataset &m_oColorDS;
375 : GDALDataset &m_oOverlayDS;
376 : const CompositionMode m_operator;
377 : const int m_opacity255Scale;
378 : std::vector<std::unique_ptr<BlendDataset>> m_apoOverviews{};
379 : int m_nCachedXOff = 0;
380 : int m_nCachedYOff = 0;
381 : int m_nCachedXSize = 0;
382 : int m_nCachedYSize = 0;
383 : int m_nCachedBufXSize = 0;
384 : int m_nCachedBufYSize = 0;
385 : GDALRasterIOExtraArg m_sCachedExtraArg{};
386 : std::vector<GByte> m_abyBuffer{};
387 : bool m_ioError = false;
388 : bool m_bSwappedOpacity = false;
389 : };
390 :
391 : /************************************************************************/
392 : /* rgb_to_hs() */
393 : /************************************************************************/
394 :
395 : // rgb comes in as [r,g,b] with values in the range [0,255]. The returned
396 : // values will be with hue and saturation in the range [0,1].
397 :
398 : // Derived from hsv_merge.py
399 :
400 4982560 : static void rgb_to_hs(int r, int g, int b, float *h, float *s)
401 : {
402 : int minc, maxc;
403 4982560 : if (r <= g)
404 : {
405 2531110 : if (r <= b)
406 : {
407 1700760 : minc = r;
408 1700760 : maxc = std::max(g, b);
409 : }
410 : else /* b < r */
411 : {
412 830353 : minc = b;
413 830353 : maxc = g;
414 : }
415 : }
416 : else /* g < r */
417 : {
418 2451460 : if (g <= b)
419 : {
420 1659840 : minc = g;
421 1659840 : maxc = std::max(r, b);
422 : }
423 : else /* b < g */
424 : {
425 791616 : minc = b;
426 791616 : maxc = r;
427 : }
428 : }
429 4982560 : const int maxc_minus_minc = maxc - minc;
430 4982560 : if (s)
431 4982560 : *s = maxc_minus_minc / static_cast<float>(std::max(1, maxc));
432 4982560 : if (h)
433 : {
434 4982560 : const float maxc_minus_minc_times_6 =
435 4982560 : maxc_minus_minc == 0 ? 1.0f : 6.0f * maxc_minus_minc;
436 4982560 : if (maxc == b)
437 1700100 : *h = 4.0f / 6.0f + (r - g) / maxc_minus_minc_times_6;
438 3282460 : else if (maxc == g)
439 1660930 : *h = 2.0f / 6.0f + (b - r) / maxc_minus_minc_times_6;
440 : else
441 : {
442 1621540 : const float tmp = (g - b) / maxc_minus_minc_times_6;
443 1621540 : *h = tmp < 0.0f ? tmp + 1.0f : tmp;
444 : }
445 : }
446 4982560 : }
447 :
448 : /************************************************************************/
449 : /* choose_among() */
450 : /************************************************************************/
451 :
452 : template <typename T>
453 5509780 : static inline T choose_among(int idx, T a0, T a1, T a2, T a3, T a4, T a5)
454 : {
455 5509780 : switch (idx)
456 : {
457 917280 : case 0:
458 917280 : return a0;
459 918255 : case 1:
460 918255 : return a1;
461 919239 : case 2:
462 919239 : return a2;
463 918978 : case 3:
464 918978 : return a3;
465 918751 : case 4:
466 918751 : return a4;
467 917280 : default:
468 917280 : break;
469 : }
470 917280 : return a5;
471 : }
472 :
473 : /************************************************************************/
474 : /* hsv_to_rgb() */
475 : /************************************************************************/
476 :
477 : // hsv comes in as [h,s,v] with hue and saturation in the range [0,1],
478 : // but value in the range [0,255].
479 :
480 : // Derived from hsv_merge.py
481 :
482 4982560 : static void hsv_to_rgb(float h, float s, GByte v, GByte *r, GByte *g, GByte *b)
483 : {
484 4982560 : const int i = static_cast<int>(6.0f * h);
485 4982560 : const float f = 6.0f * h - i;
486 4982560 : const GByte p = static_cast<GByte>(v * (1.0f - s) + 0.5f);
487 4982560 : const GByte q = static_cast<GByte>(v * (1.0f - s * f) + 0.5f);
488 4982560 : const GByte t = static_cast<GByte>(v * (1.0f - s * (1.0f - f)) + 0.5f);
489 :
490 4982560 : if (r)
491 1836600 : *r = choose_among(i, v, q, p, p, t, v);
492 4982560 : if (g)
493 1836590 : *g = choose_among(i, t, v, v, q, p, p);
494 4982560 : if (b)
495 1836590 : *b = choose_among(i, p, p, t, v, v, q);
496 4982560 : }
497 :
498 : /************************************************************************/
499 : /* XMM_RGB_to_HS() */
500 : /************************************************************************/
501 :
502 : #ifdef HAVE_SSE2
503 : static inline void
504 538808 : XMM_RGB_to_HS(const GByte *CPL_RESTRICT pInR, const GByte *CPL_RESTRICT pInG,
505 : const GByte *CPL_RESTRICT pInB, const XMMReg4Float &zero,
506 : const XMMReg4Float &one, const XMMReg4Float &six,
507 : const XMMReg4Float &two_over_six,
508 : const XMMReg4Float &four_over_six, XMMReg4Float &h,
509 : XMMReg4Float &s)
510 : {
511 538808 : const auto r = XMMReg4Float::Load4Val(pInR);
512 538808 : const auto g = XMMReg4Float::Load4Val(pInG);
513 538808 : const auto b = XMMReg4Float::Load4Val(pInB);
514 538808 : const auto minc = XMMReg4Float::Min(XMMReg4Float::Min(r, g), b);
515 538808 : const auto maxc = XMMReg4Float::Max(XMMReg4Float::Max(r, g), b);
516 538808 : const auto max_minus_min = maxc - minc;
517 538808 : s = max_minus_min / XMMReg4Float::Max(one, maxc);
518 : const auto inv_max_minus_min_times_6_0 =
519 538808 : XMMReg4Float::Ternary(XMMReg4Float::Equals(max_minus_min, zero), one,
520 538808 : six * max_minus_min)
521 538808 : .inverse();
522 538808 : const auto tmp = (g - b) * inv_max_minus_min_times_6_0;
523 538808 : h = XMMReg4Float::Ternary(
524 538808 : XMMReg4Float::Equals(maxc, b),
525 538808 : four_over_six + (r - g) * inv_max_minus_min_times_6_0,
526 538808 : XMMReg4Float::Ternary(
527 538808 : XMMReg4Float::Equals(maxc, g),
528 538808 : two_over_six + (b - r) * inv_max_minus_min_times_6_0,
529 538808 : XMMReg4Float::Ternary(XMMReg4Float::Lesser(tmp, zero), tmp + one,
530 538808 : tmp)));
531 538808 : }
532 : #endif
533 :
534 : /************************************************************************/
535 : /* patch_value_line() */
536 : /************************************************************************/
537 :
538 : static
539 : #ifdef __GNUC__
540 : __attribute__((__noinline__))
541 : #endif
542 : void
543 695 : patch_value_line(int nCount, const GByte *CPL_RESTRICT pInR,
544 : const GByte *CPL_RESTRICT pInG,
545 : const GByte *CPL_RESTRICT pInB,
546 : const GByte *CPL_RESTRICT pInGray,
547 : GByte *CPL_RESTRICT pOutR, GByte *CPL_RESTRICT pOutG,
548 : GByte *CPL_RESTRICT pOutB)
549 : {
550 695 : int i = 0;
551 : #ifdef HAVE_SSE2
552 695 : const auto zero = XMMReg4Float::Zero();
553 695 : const auto one = XMMReg4Float::Set1(1.0f);
554 695 : const auto six = XMMReg4Float::Set1(6.0f);
555 695 : const auto two_over_six = XMMReg4Float::Set1(2.0f / 6.0f);
556 695 : const auto four_over_six = two_over_six + two_over_six;
557 :
558 695 : constexpr int ELTS = 8;
559 270099 : for (; i + (ELTS - 1) < nCount; i += ELTS)
560 : {
561 269404 : XMMReg4Float h0, s0;
562 269404 : XMM_RGB_to_HS(pInR + i, pInG + i, pInB + i, zero, one, six,
563 : two_over_six, four_over_six, h0, s0);
564 269404 : XMMReg4Float h1, s1;
565 269404 : XMM_RGB_to_HS(pInR + i + ELTS / 2, pInG + i + ELTS / 2,
566 269404 : pInB + i + ELTS / 2, zero, one, six, two_over_six,
567 : four_over_six, h1, s1);
568 :
569 269404 : XMMReg4Float v0, v1;
570 269404 : XMMReg4Float::Load8Val(pInGray + i, v0, v1);
571 :
572 269404 : const auto half = XMMReg4Float::Set1(0.5f);
573 269404 : const auto six_h0 = six * h0;
574 269404 : const auto idx0 = six_h0.truncate_to_int();
575 269404 : const auto f0 = six_h0 - idx0.cast_to_float();
576 269404 : const auto p0 = (v0 * (one - s0) + half).truncate_to_int();
577 269404 : const auto q0 = (v0 * (one - s0 * f0) + half).truncate_to_int();
578 269404 : const auto t0 = (v0 * (one - s0 * (one - f0)) + half).truncate_to_int();
579 :
580 269404 : const auto six_h1 = six * h1;
581 269404 : const auto idx1 = six_h1.truncate_to_int();
582 269404 : const auto f1 = six_h1 - idx1.cast_to_float();
583 269404 : const auto p1 = (v1 * (one - s1) + half).truncate_to_int();
584 269404 : const auto q1 = (v1 * (one - s1 * f1) + half).truncate_to_int();
585 269404 : const auto t1 = (v1 * (one - s1 * (one - f1)) + half).truncate_to_int();
586 :
587 269404 : const auto idx = XMMReg8Byte::Pack(idx0, idx1);
588 : const auto v =
589 269404 : XMMReg8Byte::Pack(v0.truncate_to_int(), v1.truncate_to_int());
590 269404 : const auto p = XMMReg8Byte::Pack(p0, p1);
591 269404 : const auto q = XMMReg8Byte::Pack(q0, q1);
592 269404 : const auto t = XMMReg8Byte::Pack(t0, t1);
593 :
594 269404 : const auto equalsTo0 = XMMReg8Byte::Equals(idx, XMMReg8Byte::Zero());
595 269404 : const auto one8Byte = XMMReg8Byte::Set1(1);
596 269404 : const auto equalsTo1 = XMMReg8Byte::Equals(idx, one8Byte);
597 269404 : const auto two8Byte = one8Byte + one8Byte;
598 269404 : const auto equalsTo2 = XMMReg8Byte::Equals(idx, two8Byte);
599 269404 : const auto four8Byte = two8Byte + two8Byte;
600 269404 : const auto equalsTo4 = XMMReg8Byte::Equals(idx, four8Byte);
601 269404 : const auto equalsTo3 = XMMReg8Byte::Equals(idx, four8Byte - one8Byte);
602 : // clang-format off
603 269404 : if (pOutR)
604 : {
605 : const auto out_r =
606 : XMMReg8Byte::Ternary(equalsTo0, v,
607 134702 : XMMReg8Byte::Ternary(equalsTo1, q,
608 134702 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo2, equalsTo3), p,
609 269404 : XMMReg8Byte::Ternary(equalsTo4, t, v))));
610 134702 : out_r.Store8Val(pOutR + i);
611 : }
612 269404 : if (pOutG)
613 : {
614 : const auto out_g =
615 : XMMReg8Byte::Ternary(equalsTo0, t,
616 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo1, equalsTo2), v,
617 236636 : XMMReg8Byte::Ternary(equalsTo3, q, p)));
618 118318 : out_g.Store8Val(pOutG + i);
619 : }
620 269404 : if (pOutB)
621 : {
622 : const auto out_b =
623 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo0, equalsTo1), p,
624 118318 : XMMReg8Byte::Ternary(equalsTo2, t,
625 118318 : XMMReg8Byte::Ternary(XMMReg8Byte::Or(equalsTo3, equalsTo4),
626 118318 : v, q)));
627 118318 : out_b.Store8Val(pOutB + i);
628 : }
629 : // clang-format on
630 : }
631 : #endif
632 :
633 2524 : for (; i < nCount; ++i)
634 : {
635 : float h, s;
636 1829 : rgb_to_hs(pInR[i], pInG[i], pInB[i], &h, &s);
637 5001 : hsv_to_rgb(h, s, pInGray[i], pOutR ? pOutR + i : nullptr,
638 3172 : pOutG ? pOutG + i : nullptr, pOutB ? pOutB + i : nullptr);
639 : }
640 695 : }
641 :
642 : /************************************************************************/
643 : /* BlendBand */
644 : /************************************************************************/
645 :
646 : class BlendBand final : public GDALRasterBand
647 : {
648 : public:
649 542 : BlendBand(BlendDataset &oBlendDataset, int nBandIn)
650 542 : : m_oBlendDataset(oBlendDataset)
651 : {
652 542 : nBand = nBandIn;
653 542 : nRasterXSize = oBlendDataset.GetRasterXSize();
654 542 : nRasterYSize = oBlendDataset.GetRasterYSize();
655 542 : oBlendDataset.m_oColorDS.GetRasterBand(1)->GetBlockSize(&nBlockXSize,
656 : &nBlockYSize);
657 542 : eDataType = GDT_UInt8;
658 542 : }
659 :
660 10 : GDALColorInterp GetColorInterpretation() override
661 : {
662 10 : if (m_oBlendDataset.GetRasterCount() <= 2 && nBand == 1)
663 0 : return GCI_GrayIndex;
664 10 : else if (m_oBlendDataset.GetRasterCount() == 2 || nBand == 4)
665 0 : return GCI_AlphaBand;
666 : else
667 10 : return static_cast<GDALColorInterp>(GCI_RedBand + nBand - 1);
668 : }
669 :
670 18 : int GetOverviewCount() override
671 : {
672 18 : return static_cast<int>(m_oBlendDataset.m_apoOverviews.size());
673 : }
674 :
675 14 : GDALRasterBand *GetOverview(int idx) override
676 : {
677 13 : return idx >= 0 && idx < GetOverviewCount()
678 27 : ? m_oBlendDataset.m_apoOverviews[idx]->GetRasterBand(nBand)
679 14 : : nullptr;
680 : }
681 :
682 : protected:
683 71 : CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override
684 : {
685 71 : int nReqXSize = 0;
686 71 : int nReqYSize = 0;
687 71 : GetActualBlockSize(nBlockXOff, nBlockYOff, &nReqXSize, &nReqYSize);
688 142 : return RasterIO(GF_Read, nBlockXOff * nBlockXSize,
689 71 : nBlockYOff * nBlockYSize, nReqXSize, nReqYSize, pData,
690 71 : nReqXSize, nReqYSize, GDT_UInt8, 1, nBlockXSize,
691 142 : nullptr);
692 : }
693 :
694 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
695 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
696 : GDALDataType eBufType, GSpacing nPixelSpace,
697 : GSpacing nLineSpace,
698 : GDALRasterIOExtraArg *psExtraArg) override;
699 :
700 : private:
701 : BlendDataset &m_oBlendDataset;
702 : };
703 :
704 : /************************************************************************/
705 : /* BlendDataset::BlendDataset() */
706 : /************************************************************************/
707 :
708 163 : BlendDataset::BlendDataset(GDALDataset &oColorDS, GDALDataset &oOverlayDS,
709 : const CompositionMode eOperator,
710 163 : int nOpacity255Scale, bool bSwappedOpacity)
711 : : m_oColorDS(oColorDS), m_oOverlayDS(oOverlayDS), m_operator(eOperator),
712 163 : m_opacity255Scale(nOpacity255Scale), m_bSwappedOpacity(bSwappedOpacity)
713 : {
714 163 : m_oColorDS.Reference();
715 163 : m_oOverlayDS.Reference();
716 :
717 163 : CPLAssert(oColorDS.GetRasterXSize() == oOverlayDS.GetRasterXSize());
718 163 : CPLAssert(oColorDS.GetRasterYSize() == oOverlayDS.GetRasterYSize());
719 163 : nRasterXSize = oColorDS.GetRasterXSize();
720 163 : nRasterYSize = oColorDS.GetRasterYSize();
721 163 : const int nOvrCount = oOverlayDS.GetRasterBand(1)->GetOverviewCount();
722 163 : bool bCanCreateOvr = true;
723 705 : for (int iBand = 1; iBand <= oColorDS.GetRasterCount(); ++iBand)
724 : {
725 542 : SetBand(iBand, std::make_unique<BlendBand>(*this, iBand));
726 542 : bCanCreateOvr =
727 542 : bCanCreateOvr &&
728 542 : (iBand > oColorDS.GetRasterCount() ||
729 1626 : oColorDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount) &&
730 542 : (iBand > oOverlayDS.GetRasterCount() ||
731 449 : oOverlayDS.GetRasterBand(iBand)->GetOverviewCount() == nOvrCount);
732 : const int nColorBandxIdx =
733 542 : iBand <= oColorDS.GetRasterCount() ? iBand : 1;
734 : const int nOverlayBandIdx =
735 542 : iBand <= oOverlayDS.GetRasterCount() ? iBand : 1;
736 562 : for (int iOvr = 0; iOvr < nOvrCount && bCanCreateOvr; ++iOvr)
737 : {
738 : const auto poColorOvrBand =
739 20 : oColorDS.GetRasterBand(nColorBandxIdx)->GetOverview(iOvr);
740 : const auto poGSOvrBand =
741 20 : oOverlayDS.GetRasterBand(nOverlayBandIdx)->GetOverview(iOvr);
742 20 : bCanCreateOvr =
743 20 : poColorOvrBand->GetDataset() != &oColorDS &&
744 20 : poColorOvrBand->GetDataset() == oColorDS.GetRasterBand(1)
745 20 : ->GetOverview(iOvr)
746 20 : ->GetDataset() &&
747 20 : poGSOvrBand->GetDataset() != &oOverlayDS &&
748 20 : poGSOvrBand->GetDataset() == oOverlayDS.GetRasterBand(1)
749 20 : ->GetOverview(iOvr)
750 20 : ->GetDataset() &&
751 60 : poColorOvrBand->GetXSize() == poGSOvrBand->GetXSize() &&
752 20 : poColorOvrBand->GetYSize() == poGSOvrBand->GetYSize();
753 : }
754 : }
755 :
756 163 : SetDescription(CPLSPrintf("Blend %s width %s", m_oColorDS.GetDescription(),
757 163 : m_oOverlayDS.GetDescription()));
758 163 : if (nBands > 1)
759 : {
760 150 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
761 : }
762 :
763 163 : if (bCanCreateOvr)
764 : {
765 168 : for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
766 : {
767 5 : m_apoOverviews.push_back(std::make_unique<BlendDataset>(
768 5 : *(oColorDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),
769 5 : *(oOverlayDS.GetRasterBand(1)->GetOverview(iOvr)->GetDataset()),
770 5 : m_operator, nOpacity255Scale, bSwappedOpacity));
771 : }
772 : }
773 163 : }
774 :
775 : /************************************************************************/
776 : /* ~BlendDataset::BlendDataset() */
777 : /************************************************************************/
778 :
779 326 : BlendDataset::~BlendDataset()
780 : {
781 163 : m_oColorDS.ReleaseRef();
782 163 : m_oOverlayDS.ReleaseRef();
783 326 : }
784 :
785 : /************************************************************************/
786 : /* BlendDataset::AcquireSourcePixels() */
787 : /************************************************************************/
788 :
789 1181 : bool BlendDataset::AcquireSourcePixels(int nXOff, int nYOff, int nXSize,
790 : int nYSize, int nBufXSize, int nBufYSize,
791 : GDALRasterIOExtraArg *psExtraArg)
792 : {
793 1181 : if (nXOff == m_nCachedXOff && nYOff == m_nCachedYOff &&
794 464 : nXSize == m_nCachedXSize && nYSize == m_nCachedYSize &&
795 305 : nBufXSize == m_nCachedBufXSize && nBufYSize == m_nCachedBufYSize &&
796 305 : psExtraArg->eResampleAlg == m_sCachedExtraArg.eResampleAlg &&
797 305 : psExtraArg->bFloatingPointWindowValidity ==
798 305 : m_sCachedExtraArg.bFloatingPointWindowValidity &&
799 305 : (!psExtraArg->bFloatingPointWindowValidity ||
800 0 : (psExtraArg->dfXOff == m_sCachedExtraArg.dfXOff &&
801 0 : psExtraArg->dfYOff == m_sCachedExtraArg.dfYOff &&
802 0 : psExtraArg->dfXSize == m_sCachedExtraArg.dfXSize &&
803 0 : psExtraArg->dfYSize == m_sCachedExtraArg.dfYSize)))
804 : {
805 305 : return !m_abyBuffer.empty();
806 : }
807 :
808 876 : const int nColorCount = m_oColorDS.GetRasterCount();
809 876 : const int nOverlayCount = m_oOverlayDS.GetRasterCount();
810 876 : const int nCompsInBuffer = nColorCount + nOverlayCount;
811 876 : assert(nCompsInBuffer > 0);
812 :
813 1752 : if (static_cast<size_t>(nBufXSize) >
814 876 : std::numeric_limits<size_t>::max() / nBufYSize / nCompsInBuffer)
815 : {
816 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
817 : "Out of memory allocating temporary buffer");
818 0 : m_abyBuffer.clear();
819 0 : m_ioError = true;
820 0 : return false;
821 : }
822 :
823 876 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
824 : try
825 : {
826 876 : if (m_abyBuffer.size() < nPixelCount * nCompsInBuffer)
827 158 : m_abyBuffer.resize(nPixelCount * nCompsInBuffer);
828 : }
829 4 : catch (const std::exception &)
830 : {
831 4 : CPLError(CE_Failure, CPLE_OutOfMemory,
832 : "Out of memory allocating temporary buffer");
833 4 : m_abyBuffer.clear();
834 4 : m_ioError = true;
835 4 : return false;
836 : }
837 :
838 : const bool bOK =
839 872 : (m_oColorDS.RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
840 872 : m_abyBuffer.data(), nBufXSize, nBufYSize,
841 : GDT_UInt8, nColorCount, nullptr, 1, nBufXSize,
842 : static_cast<GSpacing>(nPixelCount),
843 1742 : psExtraArg) == CE_None &&
844 870 : m_oOverlayDS.RasterIO(
845 : GF_Read, nXOff, nYOff, nXSize, nYSize,
846 870 : m_abyBuffer.data() + nPixelCount * nColorCount, nBufXSize,
847 : nBufYSize, GDT_UInt8, nOverlayCount, nullptr, 1, nBufXSize,
848 872 : static_cast<GSpacing>(nPixelCount), psExtraArg) == CE_None);
849 872 : if (bOK)
850 : {
851 870 : m_nCachedXOff = nXOff;
852 870 : m_nCachedYOff = nYOff;
853 870 : m_nCachedXSize = nXSize;
854 870 : m_nCachedYSize = nYSize;
855 870 : m_nCachedBufXSize = nBufXSize;
856 870 : m_nCachedBufYSize = nBufYSize;
857 870 : m_sCachedExtraArg = *psExtraArg;
858 : }
859 : else
860 : {
861 2 : m_abyBuffer.clear();
862 2 : m_ioError = true;
863 : }
864 872 : return bOK;
865 : }
866 :
867 : /************************************************************************/
868 : /* gTabInvDstA */
869 : /************************************************************************/
870 :
871 : constexpr int SHIFT_DIV_DSTA = 8;
872 :
873 : // Table of (255 * 256 + k/2) / k values for k in [0,255]
874 : constexpr auto gTabInvDstA = []()
875 : {
876 : std::array<uint16_t, 256> arr{};
877 :
878 : arr[0] = 0;
879 : for (int k = 1; k <= 255; ++k)
880 : {
881 : arr[k] = static_cast<uint16_t>(((255 << SHIFT_DIV_DSTA) + (k / 2)) / k);
882 : }
883 :
884 : return arr;
885 : }();
886 :
887 : /************************************************************************/
888 : /* BlendMultiply_Generic */
889 : /************************************************************************/
890 :
891 328 : static void BlendMultiply_Generic(
892 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
893 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
894 : const GByte *CPL_RESTRICT pabyOverlayR,
895 : const GByte *CPL_RESTRICT pabyOverlayG,
896 : const GByte *CPL_RESTRICT pabyOverlayB,
897 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
898 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
899 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
900 : {
901 :
902 : // Generic formulas from Mapserver
903 : // Dca' = Sca.Dca + Sca.(1 - Da) + Dca.(1 - Sa)
904 : // Da' = Sa + Da - Sa.Da
905 :
906 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
907 : // TODO: optimize mathematically to avoid redundant computations
908 :
909 328 : GSpacing nDstOffset = 0;
910 :
911 656 : for (; i < N; ++i)
912 : {
913 :
914 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
915 : GByte nR, nG, nB, nA;
916 : GByte nFinalAlpha;
917 328 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
918 : nA, nOverlayA, nFinalAlpha);
919 328 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
920 328 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
921 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
922 :
923 : // Result
924 840 : auto processComponent = [&nFinalAlpha](GByte C, GByte A, GByte OverlayC,
925 840 : GByte OverlayA) -> GByte
926 : {
927 3360 : return DivScale255((MulScale255(C, OverlayC) +
928 1680 : MulScale255(C, 255 - OverlayA) +
929 840 : MulScale255(OverlayC, 255 - A)),
930 1680 : nFinalAlpha);
931 328 : };
932 :
933 328 : pabyDst[nDstOffset] = processComponent(nR, nA, nOverlayR, nOverlayA);
934 :
935 : // Grayscale with alpha
936 328 : if (nOutputBands == 2)
937 : {
938 48 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
939 : }
940 : else
941 : {
942 : // RBG and RGBA
943 280 : if (nOutputBands >= 3)
944 : {
945 512 : pabyDst[nDstOffset + nBandSpace] =
946 256 : processComponent(nG, nA, nOverlayG, nOverlayA);
947 512 : pabyDst[nDstOffset + 2 * nBandSpace] =
948 256 : processComponent(nB, nA, nOverlayB, nOverlayA);
949 : }
950 :
951 : // RGBA
952 280 : if (nOutputBands == 4)
953 : {
954 256 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
955 : }
956 : }
957 328 : nDstOffset += nPixelSpace;
958 : }
959 328 : }
960 :
961 : /************************************************************************/
962 : /* BlendScreen_Generic */
963 : /************************************************************************/
964 :
965 80 : static void BlendScreen_Generic(
966 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
967 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
968 : const GByte *CPL_RESTRICT pabyOverlayR,
969 : const GByte *CPL_RESTRICT pabyOverlayG,
970 : const GByte *CPL_RESTRICT pabyOverlayB,
971 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
972 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
973 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
974 : {
975 :
976 : // Generic formulas from Mapserver
977 : // Dca' = Sca + Dca - Sca.Dca
978 : // Da' = Sa + Da - Sa.Da
979 :
980 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
981 : // TODO: optimize mathematically to avoid redundant computations
982 :
983 80 : GSpacing nDstOffset = 0;
984 :
985 160 : for (; i < N; ++i)
986 : {
987 :
988 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
989 : GByte nR, nG, nB, nA;
990 : GByte nFinalAlpha;
991 80 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
992 : nA, nOverlayA, nFinalAlpha);
993 80 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
994 80 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
995 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
996 :
997 : // Result
998 :
999 240 : auto processComponent = [&nFinalAlpha](GByte C, GByte OverlayC) -> GByte
1000 : {
1001 240 : return DivScale255(C + OverlayC - MulScale255(C, OverlayC),
1002 480 : nFinalAlpha);
1003 80 : };
1004 :
1005 80 : pabyDst[nDstOffset] = processComponent(nR, nOverlayR);
1006 :
1007 : // Grayscale with alpha
1008 80 : if (nOutputBands == 2)
1009 : {
1010 0 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1011 : }
1012 : else
1013 : {
1014 : // RBG and RGBA
1015 80 : if (nOutputBands >= 3)
1016 : {
1017 160 : pabyDst[nDstOffset + nBandSpace] =
1018 80 : processComponent(nG, nOverlayG);
1019 160 : pabyDst[nDstOffset + 2 * nBandSpace] =
1020 80 : processComponent(nB, nOverlayB);
1021 : }
1022 :
1023 : // RGBA
1024 80 : if (nOutputBands == 4)
1025 : {
1026 80 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1027 : }
1028 : }
1029 80 : nDstOffset += nPixelSpace;
1030 : }
1031 80 : }
1032 :
1033 : /************************************************************************/
1034 : /* BlendOverlay_Generic */
1035 : /************************************************************************/
1036 :
1037 72 : static void BlendOverlay_Generic(
1038 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1039 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1040 : const GByte *CPL_RESTRICT pabyOverlayR,
1041 : const GByte *CPL_RESTRICT pabyOverlayG,
1042 : const GByte *CPL_RESTRICT pabyOverlayB,
1043 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1044 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1045 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1046 : {
1047 :
1048 : // Generic formulas from Mapserver
1049 : // Where "D" is destination, "S" is source (overlay)
1050 : // if 2.Dca < Da
1051 : // Dca' = 2.Sca.Dca + Sca.(1 - Da) + Dca.(1 - Sa)
1052 : // otherwise
1053 : // Dca' = Sa.Da - 2.(Da - Dca).(Sa - Sca) + Sca.(1 - Da) + Dca.(1 - Sa)
1054 : //
1055 : // Da' = Sa + Da - Sa.Da
1056 :
1057 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
1058 : // TODO: optimize mathematically to avoid redundant computations
1059 :
1060 72 : GSpacing nDstOffset = 0;
1061 :
1062 144 : for (; i < N; ++i)
1063 : {
1064 :
1065 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
1066 : GByte nR, nG, nB, nA;
1067 : GByte nFinalAlpha;
1068 72 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
1069 : nA, nOverlayA, nFinalAlpha);
1070 72 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
1071 72 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
1072 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
1073 :
1074 : // Sa.Da
1075 72 : const GByte nAlphaMul{MulScale255(nOverlayA, nA)};
1076 :
1077 112 : auto processComponent_lessThan = [&nFinalAlpha](GByte C, GByte A,
1078 : GByte OverlayC,
1079 112 : GByte OverlayA) -> GByte
1080 : {
1081 560 : return DivScale255(2 * MulScale255(C, OverlayC) +
1082 224 : MulScale255(C, 255 - OverlayA) +
1083 112 : MulScale255(OverlayC, 255 - A),
1084 224 : nFinalAlpha);
1085 72 : };
1086 :
1087 : auto processComponent_greaterEqual =
1088 104 : [&nFinalAlpha, &nAlphaMul](GByte C, GByte A, GByte OverlayC,
1089 104 : GByte OverlayA) -> GByte
1090 : {
1091 416 : return DivScale255(nAlphaMul -
1092 208 : 2 * MulScale255(A - C, OverlayA - OverlayC) +
1093 208 : MulScale255(C, 255 - OverlayA) +
1094 104 : MulScale255(OverlayC, 255 - A),
1095 208 : nFinalAlpha);
1096 72 : };
1097 :
1098 72 : if (2 * nR < nA)
1099 : {
1100 64 : pabyDst[nDstOffset] =
1101 32 : processComponent_lessThan(nR, nA, nOverlayR, nOverlayA);
1102 : }
1103 : else
1104 : {
1105 80 : pabyDst[nDstOffset] =
1106 40 : processComponent_greaterEqual(nR, nA, nOverlayR, nOverlayA);
1107 : }
1108 :
1109 : // Grayscale with alpha
1110 72 : if (nOutputBands == 2)
1111 : {
1112 0 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1113 : }
1114 : else
1115 : {
1116 : // RBG and RGBA
1117 72 : if (nOutputBands >= 3)
1118 : {
1119 :
1120 72 : if (2 * nG < nA)
1121 : {
1122 80 : pabyDst[nDstOffset + nBandSpace] =
1123 40 : processComponent_lessThan(nG, nA, nOverlayG, nOverlayA);
1124 : }
1125 : else
1126 : {
1127 64 : pabyDst[nDstOffset + nBandSpace] =
1128 32 : processComponent_greaterEqual(nG, nA, nOverlayG,
1129 : nOverlayA);
1130 : }
1131 :
1132 72 : if (2 * nB < nA)
1133 : {
1134 80 : pabyDst[nDstOffset + 2 * nBandSpace] =
1135 40 : processComponent_lessThan(nB, nA, nOverlayB, nOverlayA);
1136 : }
1137 : else
1138 : {
1139 64 : pabyDst[nDstOffset + 2 * nBandSpace] =
1140 32 : processComponent_greaterEqual(nB, nA, nOverlayB,
1141 : nOverlayA);
1142 : }
1143 : }
1144 :
1145 : // RGBA
1146 72 : if (nOutputBands == 4)
1147 : {
1148 72 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1149 : }
1150 : }
1151 72 : nDstOffset += nPixelSpace;
1152 : }
1153 72 : }
1154 :
1155 : /************************************************************************/
1156 : /* BlendHardLight_Generic */
1157 : /************************************************************************/
1158 :
1159 16 : static void BlendHardLight_Generic(
1160 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1161 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1162 : const GByte *CPL_RESTRICT pabyOverlayR,
1163 : const GByte *CPL_RESTRICT pabyOverlayG,
1164 : const GByte *CPL_RESTRICT pabyOverlayB,
1165 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1166 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1167 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1168 : {
1169 : // Hard Light is Overlay with roles of source and overlay swapped
1170 16 : BlendOverlay_Generic(pabyOverlayR, pabyOverlayG, pabyOverlayB, pabyOverlayA,
1171 : pabyR, pabyG, pabyB, pabyA, pabyDst, nPixelSpace,
1172 : nBandSpace, i, N, nOpacity, nOutputBands,
1173 16 : !bSwappedOpacity);
1174 16 : }
1175 :
1176 : /************************************************************************/
1177 : /* BlendDarken_Generic */
1178 : /************************************************************************/
1179 :
1180 88 : static void BlendDarken_Generic(
1181 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1182 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1183 : const GByte *CPL_RESTRICT pabyOverlayR,
1184 : const GByte *CPL_RESTRICT pabyOverlayG,
1185 : const GByte *CPL_RESTRICT pabyOverlayB,
1186 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1187 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1188 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1189 : {
1190 : // Generic formulas from Mapserver
1191 : // Dca' = min(Sca.Da, Dca.Sa) + Sca.(1 - Da) + Dca.(1 - Sa)
1192 : // Da' = Sa + Da - Sa.Da
1193 :
1194 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
1195 : // TODO: optimize mathematically to avoid redundant computations
1196 :
1197 88 : GSpacing nDstOffset = 0;
1198 :
1199 176 : for (; i < N; ++i)
1200 : {
1201 :
1202 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
1203 : GByte nR, nG, nB, nA;
1204 : GByte nFinalAlpha;
1205 88 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
1206 : nA, nOverlayA, nFinalAlpha);
1207 88 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
1208 88 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
1209 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
1210 :
1211 : // Result
1212 :
1213 216 : auto processComponent = [&nFinalAlpha](GByte C, GByte A, GByte OverlayC,
1214 216 : GByte OverlayA) -> GByte
1215 : {
1216 648 : return DivScale255(
1217 432 : std::min(MulScale255(OverlayC, A), MulScale255(C, OverlayA)) +
1218 432 : MulScale255(C, 255 - OverlayA) +
1219 216 : MulScale255(OverlayC, 255 - A),
1220 648 : nFinalAlpha);
1221 88 : };
1222 :
1223 88 : pabyDst[nDstOffset] = processComponent(nR, nA, nOverlayR, nOverlayA);
1224 :
1225 : // Grayscale with alpha
1226 88 : if (nOutputBands == 2)
1227 : {
1228 24 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1229 : }
1230 : else
1231 : {
1232 : // RBG and RGBA
1233 64 : if (nOutputBands >= 3)
1234 : {
1235 128 : pabyDst[nDstOffset + nBandSpace] =
1236 64 : processComponent(nG, nA, nOverlayG, nOverlayA);
1237 128 : pabyDst[nDstOffset + 2 * nBandSpace] =
1238 64 : processComponent(nB, nA, nOverlayB, nOverlayA);
1239 : }
1240 :
1241 : // RGBA
1242 64 : if (nOutputBands == 4)
1243 : {
1244 64 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1245 : }
1246 : }
1247 88 : nDstOffset += nPixelSpace;
1248 : }
1249 88 : }
1250 :
1251 : /************************************************************************/
1252 : /* BlendLighten_Generic */
1253 : /************************************************************************/
1254 :
1255 136 : static void BlendLighten_Generic(
1256 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1257 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1258 : const GByte *CPL_RESTRICT pabyOverlayR,
1259 : const GByte *CPL_RESTRICT pabyOverlayG,
1260 : const GByte *CPL_RESTRICT pabyOverlayB,
1261 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1262 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1263 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1264 : {
1265 :
1266 : // Generic formulas from Mapserver
1267 : // Dca' = max(Sca.Da, Dca.Sa) + Sca.(1 - Da) + Dca.(1 - Sa)
1268 : // Da' = Sa + Da - Sa.Da
1269 :
1270 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
1271 : // TODO: optimize mathematically to avoid redundant computations
1272 :
1273 136 : GSpacing nDstOffset = 0;
1274 :
1275 272 : for (; i < N; ++i)
1276 : {
1277 :
1278 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
1279 : GByte nR, nG, nB, nA;
1280 : GByte nFinalAlpha;
1281 136 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
1282 : nA, nOverlayA, nFinalAlpha);
1283 136 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
1284 136 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
1285 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
1286 :
1287 : // Result
1288 :
1289 360 : auto processComponent = [&nFinalAlpha](GByte C, GByte A, GByte OverlayC,
1290 360 : GByte OverlayA) -> GByte
1291 : {
1292 1080 : return DivScale255(
1293 720 : std::max(MulScale255(OverlayC, A), MulScale255(C, OverlayA)) +
1294 720 : MulScale255(C, 255 - OverlayA) +
1295 360 : MulScale255(OverlayC, 255 - A),
1296 1080 : nFinalAlpha);
1297 136 : };
1298 :
1299 136 : pabyDst[nDstOffset] = processComponent(nR, nA, nOverlayR, nOverlayA);
1300 :
1301 : // Grayscale with alpha
1302 136 : if (nOutputBands == 2)
1303 : {
1304 24 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1305 : }
1306 : else
1307 : {
1308 : // RBG and RGBA
1309 112 : if (nOutputBands >= 3)
1310 : {
1311 224 : pabyDst[nDstOffset + nBandSpace] =
1312 112 : processComponent(nG, nA, nOverlayG, nOverlayA);
1313 224 : pabyDst[nDstOffset + 2 * nBandSpace] =
1314 112 : processComponent(nB, nA, nOverlayB, nOverlayA);
1315 : }
1316 :
1317 : // RGBA
1318 112 : if (nOutputBands == 4)
1319 : {
1320 112 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1321 : }
1322 : }
1323 136 : nDstOffset += nPixelSpace;
1324 : }
1325 136 : }
1326 :
1327 : /************************************************************************/
1328 : /* BlendColorDodge_Generic */
1329 : /************************************************************************/
1330 :
1331 56 : static void BlendColorDodge_Generic(
1332 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1333 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1334 : const GByte *CPL_RESTRICT pabyOverlayR,
1335 : const GByte *CPL_RESTRICT pabyOverlayG,
1336 : const GByte *CPL_RESTRICT pabyOverlayB,
1337 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1338 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1339 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1340 : {
1341 :
1342 : // Generic formulas from Mapserver
1343 : // if Sca.Da + Dca.Sa >= Sa.Da
1344 : // Dca' = Sa.Da + Sca.(1 - Da) + Dca.(1 - Sa)
1345 : // otherwise
1346 : // Dca' = Dca.Sa/(1-Sca/Sa) + Sca.(1 - Da) + Dca.(1 - Sa)
1347 : //
1348 : // Da' = Sa + Da - Sa.Da
1349 :
1350 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
1351 : // TODO: optimize mathematically to avoid redundant computations
1352 :
1353 56 : GSpacing nDstOffset = 0;
1354 :
1355 112 : for (; i < N; ++i)
1356 : {
1357 :
1358 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
1359 : GByte nR, nG, nB, nA;
1360 : GByte nFinalAlpha;
1361 56 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
1362 : nA, nOverlayA, nFinalAlpha);
1363 56 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
1364 56 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
1365 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
1366 :
1367 : // Result
1368 :
1369 56 : const GByte alphaMul255{MulScale255(nOverlayA, nA)};
1370 :
1371 : // Dca' = Sa.Da + Sca.(1 - Da) + Dca.(1 - Sa)
1372 : auto processComponent_greaterEqual =
1373 104 : [&nFinalAlpha, &alphaMul255](GByte C, GByte A, GByte OverlayC,
1374 104 : GByte OverlayA) -> GByte
1375 : {
1376 104 : return DivScale255(alphaMul255 + MulScale255(C, 255 - OverlayA) +
1377 104 : MulScale255(OverlayC, 255 - A),
1378 208 : nFinalAlpha);
1379 56 : };
1380 :
1381 : // Dca.Sa/(1-Sca/Sa) + Sca.(1 - Da) + Dca.(1 - Sa)
1382 64 : auto processComponent_lessThan = [&nFinalAlpha](GByte C, GByte A,
1383 : GByte OverlayC,
1384 64 : GByte OverlayA) -> GByte
1385 : {
1386 192 : return DivScale255(
1387 64 : DivScale255(MulScale255(C, OverlayA),
1388 64 : 255 - DivScale255(OverlayC, OverlayA)) +
1389 128 : MulScale255(C, 255 - OverlayA) +
1390 64 : MulScale255(OverlayC, 255 - A),
1391 128 : nFinalAlpha);
1392 56 : };
1393 :
1394 168 : auto branchingCondition = [&alphaMul255](GByte C, GByte A,
1395 : GByte OverlayC,
1396 168 : GByte OverlayA) -> bool {
1397 168 : return MulScale255(OverlayC, A) + MulScale255(C, OverlayA) >=
1398 168 : alphaMul255;
1399 56 : };
1400 :
1401 56 : if (branchingCondition(nR, nA, nOverlayR, nOverlayA))
1402 : {
1403 64 : pabyDst[nDstOffset] =
1404 32 : processComponent_greaterEqual(nR, nA, nOverlayR, nOverlayA);
1405 : }
1406 : else
1407 : {
1408 48 : pabyDst[nDstOffset] =
1409 24 : processComponent_lessThan(nR, nA, nOverlayR, nOverlayA);
1410 : }
1411 :
1412 : // Grayscale with alpha
1413 56 : if (nOutputBands == 2)
1414 : {
1415 0 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1416 : }
1417 : else
1418 : {
1419 : // RBG and RGBA
1420 56 : if (nOutputBands >= 3)
1421 : {
1422 56 : if (branchingCondition(nG, nA, nOverlayG, nOverlayA))
1423 : {
1424 64 : pabyDst[nDstOffset + nBandSpace] =
1425 32 : processComponent_greaterEqual(nG, nA, nOverlayG,
1426 : nOverlayA);
1427 : }
1428 : else
1429 : {
1430 48 : pabyDst[nDstOffset + nBandSpace] =
1431 24 : processComponent_lessThan(nG, nA, nOverlayG, nOverlayA);
1432 : }
1433 :
1434 56 : if (branchingCondition(nB, nA, nOverlayB, nOverlayA))
1435 : {
1436 80 : pabyDst[nDstOffset + 2 * nBandSpace] =
1437 40 : processComponent_greaterEqual(nB, nA, nOverlayB,
1438 : nOverlayA);
1439 : }
1440 : else
1441 : {
1442 32 : pabyDst[nDstOffset + 2 * nBandSpace] =
1443 16 : processComponent_lessThan(nB, nA, nOverlayB, nOverlayA);
1444 : }
1445 : }
1446 :
1447 : // RGBA
1448 56 : if (nOutputBands == 4)
1449 : {
1450 56 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1451 : }
1452 : }
1453 56 : nDstOffset += nPixelSpace;
1454 : }
1455 56 : }
1456 :
1457 : /************************************************************************/
1458 : /* BlendColorBurn_Generic */
1459 : /************************************************************************/
1460 :
1461 64 : static void BlendColorBurn_Generic(
1462 : const GByte *CPL_RESTRICT pabyR, const GByte *CPL_RESTRICT pabyG,
1463 : const GByte *CPL_RESTRICT pabyB, const GByte *CPL_RESTRICT pabyA,
1464 : const GByte *CPL_RESTRICT pabyOverlayR,
1465 : const GByte *CPL_RESTRICT pabyOverlayG,
1466 : const GByte *CPL_RESTRICT pabyOverlayB,
1467 : const GByte *CPL_RESTRICT pabyOverlayA, GByte *CPL_RESTRICT pabyDst,
1468 : GSpacing nPixelSpace, GSpacing nBandSpace, size_t i, size_t N,
1469 : GByte nOpacity, int nOutputBands, bool bSwappedOpacity)
1470 : {
1471 :
1472 : // if Sca.Da + Dca.Sa <= Sa.Da
1473 : // Dca' = Sca.(1 - Da) + Dca.(1 - Sa)
1474 : // otherwise
1475 : // Dca' = Sa.(Sca.Da + Dca.Sa - Sa.Da)/Sca + Sca.(1 - Da) + Dca.(1 - Sa)
1476 : // Da' = Sa + Da - Sa.Da
1477 :
1478 : // TODO: optimize for the various cases (with/without alpha, grayscale/RGB)
1479 : // TODO: optimize mathematically to avoid redundant computations
1480 :
1481 64 : GSpacing nDstOffset = 0;
1482 :
1483 128 : for (; i < N; ++i)
1484 : {
1485 :
1486 : GByte nOverlayR, nOverlayG, nOverlayB, nOverlayA;
1487 : GByte nR, nG, nB, nA;
1488 : GByte nFinalAlpha;
1489 64 : ProcessAlphaChannels(i, pabyA, pabyOverlayA, nOpacity, bSwappedOpacity,
1490 : nA, nOverlayA, nFinalAlpha);
1491 64 : PremultiplyChannels(i, pabyR, pabyG, pabyB, nR, nG, nB, nA);
1492 64 : PremultiplyChannels(i, pabyOverlayR, pabyOverlayG, pabyOverlayB,
1493 : nOverlayR, nOverlayG, nOverlayB, nOverlayA);
1494 :
1495 : // Result
1496 :
1497 64 : const GByte alphaMul255{MulScale255(nOverlayA, nA)};
1498 :
1499 : auto processComponent_lessEqual =
1500 144 : [&nFinalAlpha](GByte C, GByte A, GByte OverlayC,
1501 144 : GByte OverlayA) -> GByte
1502 : {
1503 144 : return DivScale255(MulScale255(C, 255 - OverlayA) +
1504 144 : MulScale255(OverlayC, 255 - A),
1505 288 : nFinalAlpha);
1506 64 : };
1507 :
1508 : // The simplified formula below was tested and verified with the floating point equivalent [0..1] normalized version
1509 : auto processComponent_greater =
1510 48 : [&nFinalAlpha, &alphaMul255](GByte C, GByte A, GByte OverlayC,
1511 48 : GByte OverlayA) -> GByte
1512 : {
1513 48 : const GByte C_unpremultiplied = DivScale255(C, A);
1514 : const GByte Overlay_C_unpremultiplied =
1515 48 : DivScale255(OverlayC, OverlayA);
1516 144 : return DivScale255(
1517 48 : MulScale255(alphaMul255, (C_unpremultiplied +
1518 96 : Overlay_C_unpremultiplied - 255)) +
1519 96 : MulScale255(C, 255 - OverlayA) +
1520 48 : MulScale255(OverlayC, 255 - A),
1521 96 : nFinalAlpha);
1522 64 : };
1523 :
1524 192 : auto branchingCondition = [&alphaMul255](GByte C, GByte A,
1525 : GByte OverlayC,
1526 192 : GByte OverlayA) -> bool {
1527 192 : return MulScale255(OverlayC, A) + MulScale255(C, OverlayA) <=
1528 192 : alphaMul255;
1529 64 : };
1530 :
1531 64 : if (branchingCondition(nR, nA, nOverlayR, nOverlayA))
1532 : {
1533 96 : pabyDst[nDstOffset] =
1534 48 : processComponent_lessEqual(nR, nA, nOverlayR, nOverlayA);
1535 : }
1536 : else
1537 : {
1538 32 : pabyDst[nDstOffset] =
1539 16 : processComponent_greater(nR, nA, nOverlayR, nOverlayA);
1540 : }
1541 :
1542 : // Grayscale with alpha
1543 64 : if (nOutputBands == 2)
1544 : {
1545 0 : pabyDst[nDstOffset + nBandSpace] = nFinalAlpha;
1546 : }
1547 : else
1548 : {
1549 : // RBG and RGBA
1550 64 : if (nOutputBands >= 3)
1551 : {
1552 64 : if (branchingCondition(nG, nA, nOverlayG, nOverlayA))
1553 : {
1554 96 : pabyDst[nDstOffset + nBandSpace] =
1555 48 : processComponent_lessEqual(nG, nA, nOverlayG,
1556 : nOverlayA);
1557 : }
1558 : else
1559 : {
1560 32 : pabyDst[nDstOffset + nBandSpace] =
1561 16 : processComponent_greater(nG, nA, nOverlayG, nOverlayA);
1562 : }
1563 :
1564 64 : if (branchingCondition(nB, nA, nOverlayB, nOverlayA))
1565 : {
1566 96 : pabyDst[nDstOffset + 2 * nBandSpace] =
1567 48 : processComponent_lessEqual(nB, nA, nOverlayB,
1568 : nOverlayA);
1569 : }
1570 : else
1571 : {
1572 32 : pabyDst[nDstOffset + 2 * nBandSpace] =
1573 16 : processComponent_greater(nB, nA, nOverlayB, nOverlayA);
1574 : }
1575 : }
1576 :
1577 : // RGBA
1578 64 : if (nOutputBands == 4)
1579 : {
1580 64 : pabyDst[nDstOffset + 3 * nBandSpace] = nFinalAlpha;
1581 : }
1582 : }
1583 64 : nDstOffset += nPixelSpace;
1584 : }
1585 64 : }
1586 :
1587 : /************************************************************************/
1588 : /* BlendSrcOverRGBA_SSE2() */
1589 : /************************************************************************/
1590 :
1591 : #ifdef HAVE_SSE2
1592 : #if defined(__GNUC__)
1593 : __attribute__((noinline))
1594 : #endif
1595 : static int
1596 1 : BlendSrcOverRGBA_SSE2(const GByte *CPL_RESTRICT pabyR,
1597 : const GByte *CPL_RESTRICT pabyG,
1598 : const GByte *CPL_RESTRICT pabyB,
1599 : const GByte *CPL_RESTRICT pabyA,
1600 : const GByte *CPL_RESTRICT pabyOverlayR,
1601 : const GByte *CPL_RESTRICT pabyOverlayG,
1602 : const GByte *CPL_RESTRICT pabyOverlayB,
1603 : const GByte *CPL_RESTRICT pabyOverlayA,
1604 : GByte *CPL_RESTRICT pabyDst, GSpacing nBandSpace, int N,
1605 : GByte nOpacity)
1606 : {
1607 : // See scalar code after call to BlendSrcOverRGBA_SSE2() below for the
1608 : // non-obfuscated formulas...
1609 :
1610 0 : const auto load_and_unpack = [](const void *p)
1611 : {
1612 0 : const auto zero = _mm_setzero_si128();
1613 0 : auto overlayA = _mm_loadu_si128(reinterpret_cast<const __m128i *>(p));
1614 0 : return std::make_pair(_mm_unpacklo_epi8(overlayA, zero),
1615 0 : _mm_unpackhi_epi8(overlayA, zero));
1616 : };
1617 :
1618 0 : const auto pack_and_store = [](void *p, __m128i lo, __m128i hi) {
1619 0 : _mm_storeu_si128(reinterpret_cast<__m128i *>(p),
1620 : _mm_packus_epi16(lo, hi));
1621 0 : };
1622 :
1623 0 : const auto mul16bit_8bit_result = [](__m128i a, __m128i b)
1624 : {
1625 0 : const auto r255 = _mm_set1_epi16(255);
1626 0 : return _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(a, b), r255), 8);
1627 : };
1628 :
1629 2 : const auto opacity = _mm_set1_epi16(nOpacity);
1630 1 : const auto r255 = _mm_set1_epi16(255);
1631 : const int16_t *tabInvDstASigned =
1632 1 : reinterpret_cast<const int16_t *>(gTabInvDstA.data());
1633 1 : constexpr int REG_WIDTH = static_cast<int>(sizeof(opacity));
1634 :
1635 1 : int i = 0;
1636 1 : for (; i <= N - REG_WIDTH; i += REG_WIDTH)
1637 : {
1638 0 : auto [overlayA_lo, overlayA_hi] = load_and_unpack(pabyOverlayA + i);
1639 0 : auto [srcA_lo, srcA_hi] = load_and_unpack(pabyA + i);
1640 0 : overlayA_lo = mul16bit_8bit_result(overlayA_lo, opacity);
1641 0 : overlayA_hi = mul16bit_8bit_result(overlayA_hi, opacity);
1642 : auto srcAMul255MinusOverlayA_lo =
1643 0 : mul16bit_8bit_result(srcA_lo, _mm_sub_epi16(r255, overlayA_lo));
1644 : auto srcAMul255MinusOverlayA_hi =
1645 0 : mul16bit_8bit_result(srcA_hi, _mm_sub_epi16(r255, overlayA_hi));
1646 0 : auto dstA_lo = _mm_add_epi16(overlayA_lo, srcAMul255MinusOverlayA_lo);
1647 0 : auto dstA_hi = _mm_add_epi16(overlayA_hi, srcAMul255MinusOverlayA_hi);
1648 :
1649 : // This would be the equivalent of a "_mm_i16gather_epi16" operation
1650 : // which does not exist...
1651 : // invDstA_{i} = [tabInvDstASigned[dstA_{i}] for i in range(8)]
1652 0 : auto invDstA_lo = _mm_undefined_si128();
1653 0 : auto invDstA_hi = _mm_undefined_si128();
1654 : #define SET_INVDSTA(k) \
1655 : do \
1656 : { \
1657 : const int idxLo = _mm_extract_epi16(dstA_lo, k); \
1658 : const int idxHi = _mm_extract_epi16(dstA_hi, k); \
1659 : invDstA_lo = _mm_insert_epi16(invDstA_lo, tabInvDstASigned[idxLo], k); \
1660 : invDstA_hi = _mm_insert_epi16(invDstA_hi, tabInvDstASigned[idxHi], k); \
1661 : } while (0)
1662 :
1663 0 : SET_INVDSTA(0);
1664 0 : SET_INVDSTA(1);
1665 0 : SET_INVDSTA(2);
1666 0 : SET_INVDSTA(3);
1667 0 : SET_INVDSTA(4);
1668 0 : SET_INVDSTA(5);
1669 0 : SET_INVDSTA(6);
1670 0 : SET_INVDSTA(7);
1671 :
1672 0 : pack_and_store(pabyDst + i + 3 * nBandSpace, dstA_lo, dstA_hi);
1673 :
1674 : #define PROCESS_COMPONENT(pabySrc, pabyOverlay, iBand) \
1675 : do \
1676 : { \
1677 : auto [src_lo, src_hi] = load_and_unpack((pabySrc) + i); \
1678 : auto [overlay_lo, overlay_hi] = load_and_unpack((pabyOverlay) + i); \
1679 : auto dst_lo = _mm_srli_epi16( \
1680 : _mm_add_epi16( \
1681 : _mm_add_epi16( \
1682 : _mm_mullo_epi16(overlay_lo, overlayA_lo), \
1683 : _mm_mullo_epi16(src_lo, srcAMul255MinusOverlayA_lo)), \
1684 : r255), \
1685 : 8); \
1686 : auto dst_hi = _mm_srli_epi16( \
1687 : _mm_add_epi16( \
1688 : _mm_add_epi16( \
1689 : _mm_mullo_epi16(overlay_hi, overlayA_hi), \
1690 : _mm_mullo_epi16(src_hi, srcAMul255MinusOverlayA_hi)), \
1691 : r255), \
1692 : 8); \
1693 : dst_lo = mul16bit_8bit_result(dst_lo, invDstA_lo); \
1694 : dst_hi = mul16bit_8bit_result(dst_hi, invDstA_hi); \
1695 : pack_and_store(pabyDst + i + (iBand)*nBandSpace, dst_lo, dst_hi); \
1696 : } while (0)
1697 :
1698 0 : PROCESS_COMPONENT(pabyR, pabyOverlayR, 0);
1699 0 : PROCESS_COMPONENT(pabyG, pabyOverlayG, 1);
1700 0 : PROCESS_COMPONENT(pabyB, pabyOverlayB, 2);
1701 : }
1702 1 : return i;
1703 : }
1704 : #endif
1705 :
1706 : template <bool bPixelSpaceIsOne>
1707 : #if defined(__GNUC__)
1708 : __attribute__((noinline))
1709 : #endif
1710 : static void
1711 1 : BlendSrcOverRGBA_Generic(const GByte *CPL_RESTRICT pabyR,
1712 : const GByte *CPL_RESTRICT pabyG,
1713 : const GByte *CPL_RESTRICT pabyB,
1714 : const GByte *CPL_RESTRICT pabyA,
1715 : const GByte *CPL_RESTRICT pabyOverlayR,
1716 : const GByte *CPL_RESTRICT pabyOverlayG,
1717 : const GByte *CPL_RESTRICT pabyOverlayB,
1718 : const GByte *CPL_RESTRICT pabyOverlayA,
1719 : GByte *CPL_RESTRICT pabyDst,
1720 : [[maybe_unused]] GSpacing nPixelSpace,
1721 : GSpacing nBandSpace, int i, int N, GByte nOpacity)
1722 : {
1723 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
1724 : if constexpr (!bPixelSpaceIsOne)
1725 : {
1726 : assert(nPixelSpace != 1);
1727 : }
1728 : #endif
1729 1 : [[maybe_unused]] GSpacing nDstOffset = 0;
1730 : #if defined(HAVE_SSE2) && defined(__clang__) && !defined(__INTEL_CLANG_COMPILER)
1731 : // No need to vectorize for SSE2 and clang
1732 : #pragma clang loop vectorize(disable)
1733 : #endif
1734 2 : for (; i < N; ++i)
1735 : {
1736 1 : const GByte nOverlayR = pabyOverlayR[i];
1737 1 : const GByte nOverlayG = pabyOverlayG[i];
1738 1 : const GByte nOverlayB = pabyOverlayB[i];
1739 1 : const GByte nOverlayA =
1740 1 : static_cast<GByte>((pabyOverlayA[i] * nOpacity + 255) / 256);
1741 1 : const GByte nR = pabyR[i];
1742 1 : const GByte nG = pabyG[i];
1743 1 : const GByte nB = pabyB[i];
1744 1 : const GByte nA = pabyA[i];
1745 1 : const GByte nSrcAMul255MinusOverlayA =
1746 1 : static_cast<GByte>((nA * (255 - nOverlayA) + 255) / 256);
1747 1 : const GByte nDstA =
1748 : static_cast<GByte>(nOverlayA + nSrcAMul255MinusOverlayA);
1749 1 : GByte nDstR = static_cast<GByte>(
1750 1 : (nOverlayR * nOverlayA + nR * nSrcAMul255MinusOverlayA + 255) /
1751 : 256);
1752 1 : GByte nDstG = static_cast<GByte>(
1753 1 : (nOverlayG * nOverlayA + nG * nSrcAMul255MinusOverlayA + 255) /
1754 : 256);
1755 1 : GByte nDstB = static_cast<GByte>(
1756 1 : (nOverlayB * nOverlayA + nB * nSrcAMul255MinusOverlayA + 255) /
1757 : 256);
1758 : // nInvDstA = (255 << SHIFT_DIV_DSTA) / nDstA;
1759 1 : const uint16_t nInvDstA = gTabInvDstA[nDstA];
1760 1 : constexpr unsigned ROUND_OFFSET_DIV_DSTA = ((1 << SHIFT_DIV_DSTA) - 1);
1761 1 : nDstR = static_cast<GByte>((nDstR * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
1762 : SHIFT_DIV_DSTA);
1763 1 : nDstG = static_cast<GByte>((nDstG * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
1764 : SHIFT_DIV_DSTA);
1765 1 : nDstB = static_cast<GByte>((nDstB * nInvDstA + ROUND_OFFSET_DIV_DSTA) >>
1766 : SHIFT_DIV_DSTA);
1767 : if constexpr (bPixelSpaceIsOne)
1768 : {
1769 : pabyDst[i] = nDstR;
1770 : pabyDst[i + nBandSpace] = nDstG;
1771 : pabyDst[i + 2 * nBandSpace] = nDstB;
1772 : pabyDst[i + 3 * nBandSpace] = nDstA;
1773 : }
1774 : else
1775 : {
1776 1 : pabyDst[nDstOffset] = nDstR;
1777 1 : pabyDst[nDstOffset + nBandSpace] = nDstG;
1778 1 : pabyDst[nDstOffset + 2 * nBandSpace] = nDstB;
1779 1 : pabyDst[nDstOffset + 3 * nBandSpace] = nDstA;
1780 1 : nDstOffset += nPixelSpace;
1781 : }
1782 : }
1783 1 : }
1784 :
1785 : /************************************************************************/
1786 : /* BlendDataset::IRasterIO() */
1787 : /************************************************************************/
1788 :
1789 718 : CPLErr BlendDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
1790 : int nXSize, int nYSize, void *pData,
1791 : int nBufXSize, int nBufYSize,
1792 : GDALDataType eBufType, int nBandCount,
1793 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
1794 : GSpacing nLineSpace, GSpacing nBandSpace,
1795 : GDALRasterIOExtraArg *psExtraArg)
1796 : {
1797 : // Try to pass the request to the most appropriate overview dataset.
1798 718 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1799 : {
1800 2 : int bTried = FALSE;
1801 2 : const CPLErr eErr = TryOverviewRasterIO(
1802 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1803 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1804 : nBandSpace, psExtraArg, &bTried);
1805 2 : if (bTried)
1806 2 : return eErr;
1807 : }
1808 :
1809 716 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
1810 716 : const int nColorCount = m_oColorDS.GetRasterCount();
1811 716 : const int nOverlayCount = m_oOverlayDS.GetRasterCount();
1812 :
1813 : /************************************************************************/
1814 : /* HSV_VALUE */
1815 : /*************************************************************************/
1816 286 : if (nOverlayCount == 1 && m_opacity255Scale == 255 &&
1817 279 : m_operator == CompositionMode::HSV_VALUE && eRWFlag == GF_Read &&
1818 204 : eBufType == GDT_UInt8 && nBandCount == nBands &&
1819 1206 : IsAllBands(nBands, panBandMap) &&
1820 204 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
1821 : psExtraArg))
1822 : {
1823 198 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1824 198 : const GByte *pabyR = m_abyBuffer.data();
1825 198 : const GByte *pabyG = m_abyBuffer.data() + nPixelCount;
1826 198 : const GByte *pabyB = m_abyBuffer.data() + nPixelCount * 2;
1827 198 : const GByte *pabyValue = m_abyBuffer.data() + nPixelCount * nColorCount;
1828 198 : size_t nSrcIdx = 0;
1829 517 : for (int j = 0; j < nBufYSize; ++j)
1830 : {
1831 319 : auto nDstOffset = j * nLineSpace;
1832 319 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize &&
1833 317 : nBandSpace >= nLineSpace * nBufYSize)
1834 : {
1835 317 : patch_value_line(nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
1836 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
1837 317 : pabyDst + nDstOffset,
1838 317 : pabyDst + nDstOffset + nBandSpace,
1839 317 : pabyDst + nDstOffset + 2 * nBandSpace);
1840 317 : nSrcIdx += nBufXSize;
1841 : }
1842 : else
1843 : {
1844 262146 : for (int i = 0; i < nBufXSize;
1845 262144 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
1846 : {
1847 : float h, s;
1848 262144 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
1849 : &h, &s);
1850 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
1851 262144 : &pabyDst[nDstOffset + 0 * nBandSpace],
1852 262144 : &pabyDst[nDstOffset + 1 * nBandSpace],
1853 262144 : &pabyDst[nDstOffset + 2 * nBandSpace]);
1854 : }
1855 : }
1856 : }
1857 198 : if (nColorCount == 4)
1858 : {
1859 394 : for (int j = 0; j < nBufYSize; ++j)
1860 : {
1861 198 : auto nDstOffset = 3 * nBandSpace + j * nLineSpace;
1862 198 : const GByte *pabyA = m_abyBuffer.data() + nPixelCount * 3;
1863 198 : GDALCopyWords64(pabyA, GDT_UInt8, 1, pabyDst + nDstOffset,
1864 : GDT_UInt8, static_cast<int>(nPixelSpace),
1865 : nBufXSize);
1866 : }
1867 : }
1868 :
1869 198 : return CE_None;
1870 : }
1871 :
1872 : /************************************************************************/
1873 : /* SRC_OVER */
1874 : /************************************************************************/
1875 315 : else if (nOverlayCount == 4 && nColorCount == 4 &&
1876 314 : m_operator == CompositionMode::SRC_OVER && eRWFlag == GF_Read &&
1877 1 : eBufType == GDT_UInt8 && nBandCount == nBands &&
1878 834 : IsAllBands(nBands, panBandMap) &&
1879 1 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize,
1880 : nBufYSize, psExtraArg))
1881 : {
1882 1 : const GByte nOpacity = static_cast<GByte>(m_opacity255Scale);
1883 1 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1884 1 : const GByte *CPL_RESTRICT pabyR = m_abyBuffer.data();
1885 1 : const GByte *CPL_RESTRICT pabyG = m_abyBuffer.data() + nPixelCount;
1886 1 : const GByte *CPL_RESTRICT pabyB = m_abyBuffer.data() + nPixelCount * 2;
1887 1 : const GByte *CPL_RESTRICT pabyA = m_abyBuffer.data() + nPixelCount * 3;
1888 : const GByte *CPL_RESTRICT pabyOverlayR =
1889 1 : m_abyBuffer.data() + nPixelCount * nColorCount;
1890 : const GByte *CPL_RESTRICT pabyOverlayG =
1891 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
1892 : const GByte *CPL_RESTRICT pabyOverlayB =
1893 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 2);
1894 : const GByte *CPL_RESTRICT pabyOverlayA =
1895 1 : m_abyBuffer.data() + nPixelCount * (nColorCount + 3);
1896 1 : size_t nSrcIdx = 0;
1897 2 : for (int j = 0; j < nBufYSize; ++j)
1898 : {
1899 1 : auto nDstOffset = j * nLineSpace;
1900 1 : int i = 0;
1901 : #ifdef HAVE_SSE2
1902 1 : if (nPixelSpace == 1)
1903 : {
1904 2 : i = BlendSrcOverRGBA_SSE2(
1905 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
1906 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
1907 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
1908 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nBandSpace,
1909 : nBufXSize, nOpacity);
1910 1 : nSrcIdx += i;
1911 1 : nDstOffset += i;
1912 : }
1913 : #endif
1914 : #if !(defined(HAVE_SSE2) && defined(__GNUC__))
1915 : if (nPixelSpace == 1)
1916 : {
1917 : // Note: clang 20 does a rather good job at autovectorizing
1918 : // for SSE2, but BlendSrcOverRGBA_SSE2() performs better.
1919 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ true>(
1920 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
1921 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
1922 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
1923 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset,
1924 : /*nPixelSpace = */ 1, nBandSpace, i, nBufXSize, nOpacity);
1925 : }
1926 : else
1927 : #endif
1928 : {
1929 1 : BlendSrcOverRGBA_Generic</* bPixelSpaceIsOne = */ false>(
1930 : pabyR + nSrcIdx, pabyG + nSrcIdx, pabyB + nSrcIdx,
1931 : pabyA + nSrcIdx, pabyOverlayR + nSrcIdx,
1932 : pabyOverlayG + nSrcIdx, pabyOverlayB + nSrcIdx,
1933 1 : pabyOverlayA + nSrcIdx, pabyDst + nDstOffset, nPixelSpace,
1934 : nBandSpace, i, nBufXSize, nOpacity);
1935 : }
1936 1 : nSrcIdx += nBufXSize - i;
1937 : }
1938 1 : return CE_None;
1939 : }
1940 :
1941 : /************************************************************************/
1942 : /* OTHER OPERATORS */
1943 : /************************************************************************/
1944 :
1945 1335 : else if ((m_operator == CompositionMode::MULTIPLY ||
1946 301 : m_operator == CompositionMode::OVERLAY ||
1947 273 : m_operator == CompositionMode::SCREEN ||
1948 233 : m_operator == CompositionMode::HARD_LIGHT ||
1949 225 : m_operator == CompositionMode::DARKEN ||
1950 169 : m_operator == CompositionMode::LIGHTEN ||
1951 89 : m_operator == CompositionMode::COLOR_BURN ||
1952 517 : m_operator == CompositionMode::COLOR_DODGE) &&
1953 488 : eRWFlag == GF_Read && eBufType == GDT_UInt8 &&
1954 1522 : nBandCount == nBands && IsAllBands(nBands, panBandMap) &&
1955 488 : AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize, nBufXSize,
1956 : nBufYSize, psExtraArg))
1957 : {
1958 : // We should have optimized paths for 1, 2, 3 and 4 bands on input and overlay
1959 : // permutations but let's keep it simple for now.
1960 488 : const GByte nOpacity = static_cast<GByte>(m_opacity255Scale);
1961 488 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1962 488 : const GByte *CPL_RESTRICT pabyR = m_abyBuffer.data();
1963 488 : GByte *CPL_RESTRICT pabyG = nullptr;
1964 488 : GByte *CPL_RESTRICT pabyB = nullptr;
1965 488 : GByte *CPL_RESTRICT pabyA = nullptr;
1966 488 : switch (nColorCount)
1967 : {
1968 96 : case 2:
1969 96 : pabyA = m_abyBuffer.data() + nPixelCount;
1970 96 : break;
1971 0 : case 3:
1972 0 : pabyG = m_abyBuffer.data() + nPixelCount;
1973 0 : pabyB = m_abyBuffer.data() + nPixelCount * 2;
1974 0 : break;
1975 368 : case 4:
1976 368 : pabyG = m_abyBuffer.data() + nPixelCount;
1977 368 : pabyB = m_abyBuffer.data() + nPixelCount * 2;
1978 368 : pabyA = m_abyBuffer.data() + nPixelCount * 3;
1979 368 : break;
1980 : }
1981 :
1982 : const GByte *CPL_RESTRICT pabyOverlayR =
1983 488 : m_abyBuffer.data() + nPixelCount * nColorCount;
1984 488 : GByte *CPL_RESTRICT pabyOverlayG = nullptr;
1985 488 : GByte *CPL_RESTRICT pabyOverlayB = nullptr;
1986 488 : GByte *CPL_RESTRICT pabyOverlayA = nullptr;
1987 488 : switch (nOverlayCount)
1988 : {
1989 104 : case 2:
1990 104 : pabyOverlayA =
1991 104 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
1992 104 : break;
1993 0 : case 3:
1994 0 : pabyOverlayG =
1995 0 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
1996 0 : pabyOverlayB =
1997 0 : m_abyBuffer.data() + nPixelCount * (nColorCount + 2);
1998 0 : break;
1999 312 : case 4:
2000 312 : pabyOverlayG =
2001 312 : m_abyBuffer.data() + nPixelCount * (nColorCount + 1);
2002 312 : pabyOverlayB =
2003 312 : m_abyBuffer.data() + nPixelCount * (nColorCount + 2);
2004 312 : pabyOverlayA =
2005 312 : m_abyBuffer.data() + nPixelCount * (nColorCount + 3);
2006 312 : break;
2007 : }
2008 :
2009 488 : size_t nSrcIdx = 0;
2010 976 : for (int j = 0; j < nBufYSize; ++j)
2011 : {
2012 488 : auto nDstOffset = j * nLineSpace;
2013 488 : int i = 0;
2014 :
2015 488 : const GByte *CPL_RESTRICT pabyOverlayG_current =
2016 488 : pabyOverlayG ? pabyOverlayG + nSrcIdx : nullptr;
2017 488 : const GByte *CPL_RESTRICT pabyOverlayB_current =
2018 488 : pabyOverlayB ? pabyOverlayB + nSrcIdx : nullptr;
2019 488 : const GByte *CPL_RESTRICT pabyOverlayA_current =
2020 488 : pabyOverlayA ? pabyOverlayA + nSrcIdx : nullptr;
2021 :
2022 488 : const GByte *CPL_RESTRICT pabyG_current =
2023 488 : pabyG ? pabyG + nSrcIdx : nullptr;
2024 488 : const GByte *CPL_RESTRICT pabyB_current =
2025 488 : pabyB ? pabyB + nSrcIdx : nullptr;
2026 488 : const GByte *CPL_RESTRICT pabyA_current =
2027 488 : pabyA ? pabyA + nSrcIdx : nullptr;
2028 :
2029 : // Determine the number of bands
2030 488 : const int nInputBands{1 + (pabyG ? 2 : 0) + (pabyA ? 1 : 0)};
2031 488 : const int nOverlayBands{1 + (pabyOverlayG ? 2 : 0) +
2032 488 : (pabyOverlayA ? 1 : 0)};
2033 488 : const int nOutputBands = std::max(nInputBands, nOverlayBands);
2034 :
2035 488 : if (m_operator == CompositionMode::SCREEN)
2036 40 : BlendScreen_Generic(
2037 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2038 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2039 : pabyOverlayB_current, pabyOverlayA_current,
2040 40 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2041 40 : nOpacity, nOutputBands, m_bSwappedOpacity);
2042 448 : else if (m_operator == CompositionMode::MULTIPLY)
2043 216 : BlendMultiply_Generic(
2044 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2045 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2046 : pabyOverlayB_current, pabyOverlayA_current,
2047 216 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2048 216 : nOpacity, nOutputBands, m_bSwappedOpacity);
2049 232 : else if (m_operator == CompositionMode::HARD_LIGHT)
2050 8 : BlendHardLight_Generic(
2051 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2052 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2053 : pabyOverlayB_current, pabyOverlayA_current,
2054 8 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2055 8 : nOpacity, nOutputBands, m_bSwappedOpacity);
2056 224 : else if (m_operator == CompositionMode::OVERLAY)
2057 28 : BlendOverlay_Generic(
2058 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2059 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2060 : pabyOverlayB_current, pabyOverlayA_current,
2061 28 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2062 28 : nOpacity, nOutputBands, m_bSwappedOpacity);
2063 196 : else if (m_operator == CompositionMode::DARKEN)
2064 56 : BlendDarken_Generic(
2065 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2066 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2067 : pabyOverlayB_current, pabyOverlayA_current,
2068 56 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2069 56 : nOpacity, nOutputBands, m_bSwappedOpacity);
2070 140 : else if (m_operator == CompositionMode::LIGHTEN)
2071 80 : BlendLighten_Generic(
2072 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2073 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2074 : pabyOverlayB_current, pabyOverlayA_current,
2075 80 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2076 80 : nOpacity, nOutputBands, m_bSwappedOpacity);
2077 60 : else if (m_operator == CompositionMode::COLOR_BURN)
2078 32 : BlendColorBurn_Generic(
2079 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2080 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2081 : pabyOverlayB_current, pabyOverlayA_current,
2082 32 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2083 32 : nOpacity, nOutputBands, m_bSwappedOpacity);
2084 28 : else if (m_operator == CompositionMode::COLOR_DODGE)
2085 28 : BlendColorDodge_Generic(
2086 : pabyR + nSrcIdx, pabyG_current, pabyB_current,
2087 : pabyA_current, pabyOverlayR + nSrcIdx, pabyOverlayG_current,
2088 : pabyOverlayB_current, pabyOverlayA_current,
2089 28 : pabyDst + nDstOffset, nPixelSpace, nBandSpace, i, nBufXSize,
2090 28 : nOpacity, nOutputBands, m_bSwappedOpacity);
2091 : else
2092 : {
2093 0 : CPLAssert(false);
2094 : }
2095 :
2096 488 : nSrcIdx += nBufXSize - i;
2097 : }
2098 488 : return CE_None;
2099 : }
2100 :
2101 : /************************************************************************/
2102 : /* ERRORS */
2103 : /************************************************************************/
2104 29 : else if (m_ioError)
2105 : {
2106 6 : return CE_Failure;
2107 : }
2108 : else
2109 : {
2110 23 : const CPLErr eErr = GDALDataset::IRasterIO(
2111 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
2112 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
2113 : nBandSpace, psExtraArg);
2114 23 : m_ioError = eErr != CE_None;
2115 23 : return eErr;
2116 : }
2117 : }
2118 :
2119 : /************************************************************************/
2120 : /* SrcOverRGBOneComponent() */
2121 : /************************************************************************/
2122 :
2123 : // GCC and clang do a god job a auto vectorizing the below function
2124 : #if defined(__GNUC__) && !defined(__clang__)
2125 : __attribute__((optimize("tree-vectorize")))
2126 : #endif
2127 : static void
2128 12 : SrcOverRGB(const uint8_t *const __restrict pabyOverlay,
2129 : const uint8_t *const __restrict pabySrc,
2130 : uint8_t *const __restrict pabyDst, const size_t N,
2131 : const uint8_t nOpacity)
2132 : {
2133 915 : for (size_t i = 0; i < N; ++i)
2134 : {
2135 903 : const uint8_t nOverlay = pabyOverlay[i];
2136 903 : const uint8_t nSrc = pabySrc[i];
2137 903 : pabyDst[i] = static_cast<uint8_t>(
2138 903 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) / 256);
2139 : }
2140 12 : }
2141 :
2142 : /************************************************************************/
2143 : /* BlendBand::IRasterIO() */
2144 : /************************************************************************/
2145 :
2146 566 : CPLErr BlendBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
2147 : int nXSize, int nYSize, void *pData, int nBufXSize,
2148 : int nBufYSize, GDALDataType eBufType,
2149 : GSpacing nPixelSpace, GSpacing nLineSpace,
2150 : GDALRasterIOExtraArg *psExtraArg)
2151 : {
2152 : // Try to pass the request to the most appropriate overview dataset.
2153 566 : if (nBufXSize < nXSize && nBufYSize < nYSize)
2154 : {
2155 1 : int bTried = FALSE;
2156 1 : const CPLErr eErr = TryOverviewRasterIO(
2157 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
2158 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
2159 1 : if (bTried)
2160 1 : return eErr;
2161 : }
2162 :
2163 565 : const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
2164 565 : const int nColorCount = m_oBlendDataset.m_oColorDS.GetRasterCount();
2165 565 : const int nOverlayCount = m_oBlendDataset.m_oOverlayDS.GetRasterCount();
2166 565 : if (nBand == 4 && m_oBlendDataset.m_operator == CompositionMode::HSV_VALUE)
2167 : {
2168 9 : if (nColorCount == 3)
2169 : {
2170 0 : GByte ch = 255;
2171 0 : for (int iY = 0; iY < nBufYSize; ++iY)
2172 : {
2173 0 : GDALCopyWords64(&ch, GDT_UInt8, 0,
2174 0 : static_cast<GByte *>(pData) + iY * nLineSpace,
2175 : eBufType, static_cast<int>(nPixelSpace),
2176 : nBufXSize);
2177 : }
2178 0 : return CE_None;
2179 : }
2180 : else
2181 : {
2182 9 : CPLAssert(nColorCount == 4);
2183 9 : return m_oBlendDataset.m_oColorDS.GetRasterBand(4)->RasterIO(
2184 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
2185 9 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg);
2186 : }
2187 : }
2188 25 : else if (nOverlayCount == 3 && nColorCount == 3 &&
2189 18 : m_oBlendDataset.m_operator == CompositionMode::SRC_OVER &&
2190 593 : eRWFlag == GF_Read && eBufType == GDT_UInt8 &&
2191 12 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
2192 : nBufXSize, nBufYSize,
2193 : psExtraArg))
2194 : {
2195 12 : const int nOpacity = m_oBlendDataset.m_opacity255Scale;
2196 : const GByte *const CPL_RESTRICT pabySrc =
2197 12 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * (nBand - 1);
2198 : const GByte *const CPL_RESTRICT pabyOverlay =
2199 12 : m_oBlendDataset.m_abyBuffer.data() +
2200 12 : nPixelCount * (nColorCount + nBand - 1);
2201 12 : GByte *const CPL_RESTRICT pabyDst = static_cast<GByte *>(pData);
2202 12 : size_t nSrcIdx = 0;
2203 24 : for (int j = 0; j < nBufYSize; ++j)
2204 : {
2205 12 : auto nDstOffset = j * nLineSpace;
2206 12 : if (nPixelSpace == 1)
2207 : {
2208 12 : SrcOverRGB(pabyOverlay + nSrcIdx, pabySrc + nSrcIdx,
2209 12 : pabyDst + nDstOffset, nBufXSize,
2210 : static_cast<uint8_t>(nOpacity));
2211 12 : nSrcIdx += nBufXSize;
2212 : }
2213 : else
2214 : {
2215 0 : for (int i = 0; i < nBufXSize;
2216 0 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
2217 : {
2218 0 : const int nOverlay = pabyOverlay[nSrcIdx];
2219 0 : const int nSrc = pabySrc[nSrcIdx];
2220 0 : pabyDst[nDstOffset] = static_cast<GByte>(
2221 0 : (nOverlay * nOpacity + nSrc * (255 - nOpacity) + 255) /
2222 : 256);
2223 : }
2224 : }
2225 : }
2226 12 : return CE_None;
2227 : }
2228 1020 : else if (eRWFlag == GF_Read && eBufType == GDT_UInt8 &&
2229 476 : m_oBlendDataset.AcquireSourcePixels(nXOff, nYOff, nXSize, nYSize,
2230 : nBufXSize, nBufYSize,
2231 : psExtraArg))
2232 : {
2233 476 : GByte *pabyDst = static_cast<GByte *>(pData);
2234 :
2235 476 : if (m_oBlendDataset.m_operator == CompositionMode::MULTIPLY ||
2236 364 : m_oBlendDataset.m_operator == CompositionMode::SCREEN ||
2237 324 : m_oBlendDataset.m_operator == CompositionMode::HARD_LIGHT ||
2238 316 : m_oBlendDataset.m_operator == CompositionMode::OVERLAY ||
2239 288 : m_oBlendDataset.m_operator == CompositionMode::DARKEN ||
2240 256 : m_oBlendDataset.m_operator == CompositionMode::LIGHTEN ||
2241 200 : m_oBlendDataset.m_operator == CompositionMode::COLOR_BURN ||
2242 168 : m_oBlendDataset.m_operator == CompositionMode::COLOR_DODGE)
2243 : {
2244 336 : CPLAssert(nBand <= 4);
2245 336 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
2246 : const GByte *pabyG =
2247 336 : nColorCount >= 3
2248 336 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
2249 336 : : nullptr;
2250 : const GByte *pabyB =
2251 336 : nColorCount >= 3
2252 336 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2
2253 336 : : nullptr;
2254 :
2255 336 : GByte *pabyA = nullptr;
2256 :
2257 336 : if (nColorCount == 2)
2258 : {
2259 0 : pabyA = m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
2260 : }
2261 336 : else if (nColorCount == 4)
2262 : {
2263 336 : pabyA = m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 3;
2264 : }
2265 :
2266 : // Retrieve single band value as R
2267 : const GByte *pabyOverlayR =
2268 336 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
2269 :
2270 : const GByte *pabyOverlayG =
2271 336 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
2272 312 : nPixelCount * (nColorCount + 1)
2273 336 : : nullptr;
2274 : const GByte *pabyOverlayB =
2275 336 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
2276 312 : nPixelCount * (nColorCount + 2)
2277 336 : : nullptr;
2278 : const GByte *pabyOverlayA =
2279 336 : (nOverlayCount == 2 || nOverlayCount == 4)
2280 672 : ? m_oBlendDataset.m_abyBuffer.data() +
2281 312 : nPixelCount * (nColorCount + nOverlayCount - 1)
2282 336 : : nullptr;
2283 :
2284 : // Determine the number of bands
2285 336 : const int nInputBands{1 + (pabyG ? 2 : 0) + (pabyA ? 1 : 0)};
2286 336 : const int nOverlayBands{1 + (pabyOverlayG ? 2 : 0) +
2287 336 : (pabyOverlayA ? 1 : 0)};
2288 336 : const int nOutputBands = std::max(nInputBands, nOverlayBands);
2289 336 : const bool bSwappedOpacity{m_oBlendDataset.m_bSwappedOpacity};
2290 :
2291 336 : size_t nSrcIdx = 0;
2292 672 : for (int j = 0; j < nBufYSize; ++j)
2293 : {
2294 336 : auto nDstOffset = j * nLineSpace;
2295 672 : for (int i = 0; i < nBufXSize;
2296 336 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
2297 : {
2298 : // TODO: This need to be optimized for requesting a single band
2299 336 : std::vector<GByte> byteBuffer;
2300 336 : byteBuffer.resize(std::max(nColorCount, nOverlayCount));
2301 336 : if (m_oBlendDataset.m_operator == CompositionMode::SCREEN)
2302 40 : BlendScreen_Generic(
2303 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2304 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2305 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2306 : static_cast<GByte>(
2307 40 : m_oBlendDataset.m_opacity255Scale),
2308 : nOutputBands, bSwappedOpacity);
2309 296 : else if (m_oBlendDataset.m_operator ==
2310 : CompositionMode::HARD_LIGHT)
2311 8 : BlendHardLight_Generic(
2312 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2313 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2314 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2315 : static_cast<GByte>(
2316 8 : m_oBlendDataset.m_opacity255Scale),
2317 : nOutputBands, bSwappedOpacity);
2318 288 : else if (m_oBlendDataset.m_operator ==
2319 : CompositionMode::OVERLAY)
2320 28 : BlendOverlay_Generic(
2321 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2322 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2323 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2324 : static_cast<GByte>(
2325 28 : m_oBlendDataset.m_opacity255Scale),
2326 : nOutputBands, bSwappedOpacity);
2327 260 : else if (m_oBlendDataset.m_operator ==
2328 : CompositionMode::MULTIPLY)
2329 112 : BlendMultiply_Generic(
2330 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2331 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2332 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2333 : static_cast<GByte>(
2334 112 : m_oBlendDataset.m_opacity255Scale),
2335 : nOutputBands, bSwappedOpacity);
2336 148 : else if (m_oBlendDataset.m_operator ==
2337 : CompositionMode::DARKEN)
2338 32 : BlendDarken_Generic(
2339 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2340 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2341 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2342 : static_cast<GByte>(
2343 32 : m_oBlendDataset.m_opacity255Scale),
2344 : nOutputBands, bSwappedOpacity);
2345 116 : else if (m_oBlendDataset.m_operator ==
2346 : CompositionMode::LIGHTEN)
2347 56 : BlendLighten_Generic(
2348 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2349 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2350 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2351 : static_cast<GByte>(
2352 56 : m_oBlendDataset.m_opacity255Scale),
2353 : nOutputBands, bSwappedOpacity);
2354 60 : else if (m_oBlendDataset.m_operator ==
2355 : CompositionMode::COLOR_BURN)
2356 32 : BlendColorBurn_Generic(
2357 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2358 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2359 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2360 : static_cast<GByte>(
2361 32 : m_oBlendDataset.m_opacity255Scale),
2362 : nOutputBands, bSwappedOpacity);
2363 28 : else if (m_oBlendDataset.m_operator ==
2364 : CompositionMode::COLOR_DODGE)
2365 28 : BlendColorDodge_Generic(
2366 : pabyR, pabyG, pabyB, pabyA, pabyOverlayR,
2367 : pabyOverlayG, pabyOverlayB, pabyOverlayA,
2368 : byteBuffer.data(), nPixelSpace, 1, nSrcIdx, 1,
2369 : static_cast<GByte>(
2370 28 : m_oBlendDataset.m_opacity255Scale),
2371 : nOutputBands, bSwappedOpacity);
2372 : else
2373 : {
2374 0 : CPLAssert(false);
2375 : }
2376 :
2377 336 : pabyDst[nDstOffset] = byteBuffer[nBand - 1];
2378 : }
2379 336 : }
2380 : }
2381 140 : else if (m_oBlendDataset.m_operator == CompositionMode::SRC_OVER)
2382 : {
2383 3 : const auto RGBToGrayScale = [](int R, int G, int B)
2384 : {
2385 : // Equivalent to R * 0.299 + G * 0.587 + B * 0.114
2386 : // but using faster computation
2387 3 : return (R * 306 + G * 601 + B * 117) / 1024;
2388 : };
2389 :
2390 : const GByte *paby =
2391 83 : (nBand <= nColorCount) ? m_oBlendDataset.m_abyBuffer.data() +
2392 83 : nPixelCount * (nBand - 1)
2393 0 : : (nBand == 4 && nColorCount == 2)
2394 0 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
2395 83 : : nullptr;
2396 : const GByte *pabyA =
2397 83 : (nColorCount == 4)
2398 106 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 3
2399 23 : : nColorCount == 2
2400 4 : ? m_oBlendDataset.m_abyBuffer.data() + nPixelCount
2401 83 : : nullptr;
2402 : const GByte *pabyOverlay =
2403 83 : (nBand <= nOverlayCount)
2404 91 : ? m_oBlendDataset.m_abyBuffer.data() +
2405 75 : nPixelCount * (nColorCount + nBand - 1)
2406 8 : : (nBand <= 3) ? m_oBlendDataset.m_abyBuffer.data() +
2407 7 : nPixelCount * nColorCount
2408 83 : : nullptr;
2409 : const GByte *pabyOverlayA =
2410 80 : (nOverlayCount == 2 || nOverlayCount == 4)
2411 163 : ? m_oBlendDataset.m_abyBuffer.data() +
2412 62 : nPixelCount * (nColorCount + nOverlayCount - 1)
2413 83 : : nullptr;
2414 : const GByte *pabyOverlayR =
2415 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
2416 149 : ? m_oBlendDataset.m_abyBuffer.data() +
2417 3 : nPixelCount * nColorCount
2418 83 : : nullptr;
2419 : const GByte *pabyOverlayG =
2420 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
2421 149 : ? m_oBlendDataset.m_abyBuffer.data() +
2422 3 : nPixelCount * (nColorCount + 1)
2423 83 : : nullptr;
2424 : const GByte *pabyOverlayB =
2425 66 : (nOverlayCount >= 3 && nColorCount < 3 && nBand <= 3)
2426 149 : ? m_oBlendDataset.m_abyBuffer.data() +
2427 3 : nPixelCount * (nColorCount + 2)
2428 83 : : nullptr;
2429 :
2430 83 : size_t nSrcIdx = 0;
2431 714 : for (int j = 0; j < nBufYSize; ++j)
2432 : {
2433 631 : auto nDstOffset = j * nLineSpace;
2434 98753 : for (int i = 0; i < nBufXSize;
2435 98122 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
2436 : {
2437 : // Corrected to take into account m_opacity255Scale
2438 98122 : const int nOverlayA =
2439 98122 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
2440 97606 : m_oBlendDataset.m_opacity255Scale +
2441 : 255) /
2442 : 256)
2443 516 : : m_oBlendDataset.m_opacity255Scale;
2444 :
2445 98122 : const int nSrcA = pabyA ? pabyA[nSrcIdx] : 255;
2446 :
2447 98122 : const int nSrcAMul255MinusOverlayA =
2448 98122 : (nSrcA * (255 - nOverlayA) + 255) / 256;
2449 98122 : const int nDstA = nOverlayA + nSrcAMul255MinusOverlayA;
2450 98122 : if (nBand != 4)
2451 : {
2452 : const int nOverlay =
2453 3 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
2454 147439 : ? RGBToGrayScale(pabyOverlayR[nSrcIdx],
2455 3 : pabyOverlayG[nSrcIdx],
2456 3 : pabyOverlayB[nSrcIdx])
2457 73718 : : pabyOverlay ? pabyOverlay[nSrcIdx]
2458 73721 : : 255;
2459 :
2460 73721 : const int nSrc = paby ? paby[nSrcIdx] : 255;
2461 73721 : int nDst = (nOverlay * nOverlayA +
2462 73721 : nSrc * nSrcAMul255MinusOverlayA + 255) /
2463 : 256;
2464 73721 : if (nDstA != 0 && nDstA != 255)
2465 19672 : nDst = (nDst * 255 + (nDstA / 2)) / nDstA;
2466 73721 : pabyDst[nDstOffset] =
2467 73721 : static_cast<GByte>(std::min(nDst, 255));
2468 : }
2469 : else
2470 : {
2471 24401 : pabyDst[nDstOffset] = static_cast<GByte>(nDstA);
2472 : }
2473 : }
2474 : }
2475 : }
2476 57 : else if (nOverlayCount == 1 && m_oBlendDataset.m_opacity255Scale == 255)
2477 : {
2478 27 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
2479 : const GByte *pabyG =
2480 27 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
2481 : const GByte *pabyB =
2482 27 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
2483 27 : CPLAssert(m_oBlendDataset.m_operator == CompositionMode::HSV_VALUE);
2484 27 : size_t nSrcIdx = 0;
2485 : const GByte *pabyValue =
2486 27 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
2487 411 : for (int j = 0; j < nBufYSize; ++j)
2488 : {
2489 384 : auto nDstOffset = j * nLineSpace;
2490 384 : if (nPixelSpace == 1 && nLineSpace >= nPixelSpace * nBufXSize)
2491 : {
2492 884 : patch_value_line(
2493 : nBufXSize, pabyR + nSrcIdx, pabyG + nSrcIdx,
2494 : pabyB + nSrcIdx, pabyValue + nSrcIdx,
2495 378 : nBand == 1 ? pabyDst + nDstOffset : nullptr,
2496 378 : nBand == 2 ? pabyDst + nDstOffset : nullptr,
2497 378 : nBand == 3 ? pabyDst + nDstOffset : nullptr);
2498 378 : nSrcIdx += nBufXSize;
2499 : }
2500 : else
2501 : {
2502 786438 : for (int i = 0; i < nBufXSize;
2503 786432 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
2504 : {
2505 : float h, s;
2506 786432 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx],
2507 786432 : pabyB[nSrcIdx], &h, &s);
2508 786432 : if (nBand == 1)
2509 : {
2510 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx],
2511 262144 : &pabyDst[nDstOffset], nullptr, nullptr);
2512 : }
2513 524288 : else if (nBand == 2)
2514 : {
2515 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
2516 262144 : &pabyDst[nDstOffset], nullptr);
2517 : }
2518 : else
2519 : {
2520 262144 : CPLAssert(nBand == 3);
2521 262144 : hsv_to_rgb(h, s, pabyValue[nSrcIdx], nullptr,
2522 262144 : nullptr, &pabyDst[nDstOffset]);
2523 : }
2524 : }
2525 : }
2526 27 : }
2527 : }
2528 : else
2529 : {
2530 30 : CPLAssert(m_oBlendDataset.m_operator == CompositionMode::HSV_VALUE);
2531 30 : CPLAssert(nBand <= 3);
2532 30 : const GByte *pabyR = m_oBlendDataset.m_abyBuffer.data();
2533 : const GByte *pabyG =
2534 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount;
2535 : const GByte *pabyB =
2536 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * 2;
2537 : const GByte *pabyValue =
2538 30 : m_oBlendDataset.m_abyBuffer.data() + nPixelCount * nColorCount;
2539 : const GByte *pabyOverlayR =
2540 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
2541 12 : nPixelCount * nColorCount
2542 30 : : nullptr;
2543 : const GByte *pabyOverlayG =
2544 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
2545 12 : nPixelCount * (nColorCount + 1)
2546 30 : : nullptr;
2547 : const GByte *pabyOverlayB =
2548 30 : nOverlayCount >= 3 ? m_oBlendDataset.m_abyBuffer.data() +
2549 12 : nPixelCount * (nColorCount + 2)
2550 30 : : nullptr;
2551 : const GByte *pabyOverlayA =
2552 24 : (nOverlayCount == 2 || nOverlayCount == 4)
2553 54 : ? m_oBlendDataset.m_abyBuffer.data() +
2554 12 : nPixelCount * (nColorCount + nOverlayCount - 1)
2555 30 : : nullptr;
2556 :
2557 30 : size_t nSrcIdx = 0;
2558 60 : for (int j = 0; j < nBufYSize; ++j)
2559 : {
2560 30 : auto nDstOffset = j * nLineSpace;
2561 3932190 : for (int i = 0; i < nBufXSize;
2562 3932160 : ++i, ++nSrcIdx, nDstOffset += nPixelSpace)
2563 : {
2564 3932160 : const int nColorR = pabyR[nSrcIdx];
2565 3932160 : const int nColorG = pabyG[nSrcIdx];
2566 3932160 : const int nColorB = pabyB[nSrcIdx];
2567 : const int nOverlayV =
2568 1572860 : (pabyOverlayR && pabyOverlayG && pabyOverlayB)
2569 5505020 : ? std::max({pabyOverlayR[nSrcIdx],
2570 1572860 : pabyOverlayG[nSrcIdx],
2571 1572860 : pabyOverlayB[nSrcIdx]})
2572 2359300 : : pabyValue[nSrcIdx];
2573 3932160 : const int nOverlayA =
2574 3932160 : pabyOverlayA ? ((pabyOverlayA[nSrcIdx] *
2575 1572860 : m_oBlendDataset.m_opacity255Scale +
2576 : 255) /
2577 : 256)
2578 2359300 : : m_oBlendDataset.m_opacity255Scale;
2579 : const int nColorValue =
2580 3932160 : std::max({nColorR, nColorG, nColorB});
2581 :
2582 : float h, s;
2583 3932160 : rgb_to_hs(pabyR[nSrcIdx], pabyG[nSrcIdx], pabyB[nSrcIdx],
2584 : &h, &s);
2585 :
2586 3932160 : const GByte nTargetValue = static_cast<GByte>(
2587 3932160 : (nOverlayV * nOverlayA +
2588 3932160 : nColorValue * (255 - nOverlayA) + 255) /
2589 : 256);
2590 :
2591 3932160 : if (nBand == 1)
2592 : {
2593 1310720 : hsv_to_rgb(h, s, nTargetValue, &pabyDst[nDstOffset],
2594 : nullptr, nullptr);
2595 : }
2596 2621440 : else if (nBand == 2)
2597 : {
2598 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr,
2599 1310720 : &pabyDst[nDstOffset], nullptr);
2600 : }
2601 : else
2602 : {
2603 1310720 : CPLAssert(nBand == 3);
2604 1310720 : hsv_to_rgb(h, s, nTargetValue, nullptr, nullptr,
2605 1310720 : &pabyDst[nDstOffset]);
2606 : }
2607 : }
2608 : }
2609 : }
2610 :
2611 476 : return CE_None;
2612 : }
2613 68 : else if (m_oBlendDataset.m_ioError)
2614 : {
2615 0 : return CE_Failure;
2616 : }
2617 : else
2618 : {
2619 68 : const CPLErr eErr = GDALRasterBand::IRasterIO(
2620 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
2621 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
2622 68 : m_oBlendDataset.m_ioError = eErr != CE_None;
2623 68 : return eErr;
2624 : }
2625 : }
2626 :
2627 : } // namespace
2628 :
2629 : /************************************************************************/
2630 : /* GDALRasterBlendAlgorithm::ValidateGlobal() */
2631 : /************************************************************************/
2632 :
2633 340 : bool GDALRasterBlendAlgorithm::ValidateGlobal()
2634 : {
2635 : auto poSrcDS =
2636 340 : m_inputDataset.empty() ? nullptr : m_inputDataset[0].GetDatasetRef();
2637 340 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
2638 340 : if (poSrcDS)
2639 : {
2640 674 : if (poSrcDS->GetRasterCount() == 0 || poSrcDS->GetRasterCount() > 4 ||
2641 337 : poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_UInt8)
2642 : {
2643 1 : ReportError(CE_Failure, CPLE_NotSupported,
2644 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
2645 : "supported as input");
2646 1 : return false;
2647 : }
2648 : }
2649 339 : if (poOverlayDS)
2650 : {
2651 339 : if (poOverlayDS->GetRasterCount() == 0 ||
2652 678 : poOverlayDS->GetRasterCount() > 4 ||
2653 339 : poOverlayDS->GetRasterBand(1)->GetRasterDataType() != GDT_UInt8)
2654 : {
2655 1 : ReportError(CE_Failure, CPLE_NotSupported,
2656 : "Only 1-band, 2-band, 3-band or 4-band Byte dataset "
2657 : "supported as overlay");
2658 1 : return false;
2659 : }
2660 : }
2661 :
2662 338 : if (poSrcDS && poOverlayDS)
2663 : {
2664 669 : if (poSrcDS->GetRasterXSize() != poOverlayDS->GetRasterXSize() ||
2665 334 : poSrcDS->GetRasterYSize() != poOverlayDS->GetRasterYSize())
2666 : {
2667 2 : ReportError(CE_Failure, CPLE_IllegalArg,
2668 : "Input dataset and overlay dataset must have "
2669 : "the same dimensions");
2670 2 : return false;
2671 : }
2672 :
2673 333 : if (!BandCountIsCompatibleWithCompositionMode(poSrcDS->GetRasterCount(),
2674 : m_operator))
2675 : {
2676 : const int minRequiredBands{
2677 1 : MinBandCountForCompositionMode(m_operator)};
2678 : const int maxRequiredBands{
2679 1 : MaxBandCountForCompositionMode(m_operator)};
2680 1 : if (minRequiredBands != maxRequiredBands)
2681 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2682 : "Input dataset has %d band(s), but operator %s "
2683 : "requires between %d and %d bands",
2684 : poSrcDS->GetRasterCount(),
2685 2 : CompositionModeToString(m_operator).c_str(),
2686 : minRequiredBands, maxRequiredBands);
2687 : else
2688 0 : ReportError(CE_Failure, CPLE_IllegalArg,
2689 : "Input dataset has %d band(s), but operator %s "
2690 : "requires %d bands",
2691 : poSrcDS->GetRasterCount(),
2692 0 : CompositionModeToString(m_operator).c_str(),
2693 : minRequiredBands);
2694 1 : return false;
2695 : }
2696 : }
2697 :
2698 : // Check that for LIGHTEN and DARKEN, the source dataset and destination dataset
2699 : // have the same number of color bands (do not consider alpha)
2700 335 : if (poSrcDS && poOverlayDS &&
2701 332 : (m_operator == CompositionMode::LIGHTEN ||
2702 284 : m_operator == CompositionMode::DARKEN))
2703 : {
2704 : const int nSrcColorBands =
2705 140 : (poSrcDS->GetRasterCount() == 2 || poSrcDS->GetRasterCount() == 4)
2706 132 : ? poSrcDS->GetRasterCount() - 1
2707 8 : : poSrcDS->GetRasterCount();
2708 84 : const int nOverlayColorBands = (poOverlayDS->GetRasterCount() == 2 ||
2709 56 : poOverlayDS->GetRasterCount() == 4)
2710 132 : ? poOverlayDS->GetRasterCount() - 1
2711 8 : : poOverlayDS->GetRasterCount();
2712 84 : if (nSrcColorBands != nOverlayColorBands)
2713 : {
2714 16 : ReportError(
2715 : CE_Failure, CPLE_IllegalArg,
2716 : "For LIGHTEN and DARKEN operators, the source dataset "
2717 : "and overlay dataset must have the same number of "
2718 : "bands (without considering alpha). They have %d and %d "
2719 : "bands respectively",
2720 : nSrcColorBands, nOverlayColorBands);
2721 16 : return false;
2722 : }
2723 : }
2724 :
2725 319 : return true;
2726 : }
2727 :
2728 : /************************************************************************/
2729 : /* GDALRasterBlendAlgorithm::RunStep() */
2730 : /************************************************************************/
2731 :
2732 158 : bool GDALRasterBlendAlgorithm::RunStep(GDALPipelineStepRunContext &)
2733 : {
2734 158 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
2735 158 : CPLAssert(poSrcDS);
2736 :
2737 158 : auto poOverlayDS = m_overlayDataset.GetDatasetRef();
2738 158 : CPLAssert(poOverlayDS);
2739 :
2740 : // If any of the dataset single band has a color table implicitly convert it to RGBA by calling
2741 : // GDALTranslate with -expand RGBA
2742 : auto convertToRGBAifNeeded =
2743 316 : [](GDALDataset *&poDS,
2744 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser> &ds)
2745 : -> bool
2746 : {
2747 374 : if (poDS->GetRasterCount() == 1 &&
2748 58 : poDS->GetRasterBand(1)->GetColorTable() != nullptr)
2749 : {
2750 12 : CPLStringList aosOptions;
2751 6 : aosOptions.AddString("-of");
2752 6 : aosOptions.AddString("VRT");
2753 6 : aosOptions.AddString("-expand");
2754 6 : aosOptions.AddString("RGBA");
2755 : GDALTranslateOptions *translateOptions =
2756 6 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
2757 :
2758 6 : ds.reset(GDALDataset::FromHandle(GDALTranslate(
2759 : "", GDALDataset::ToHandle(poDS), translateOptions, nullptr)));
2760 :
2761 6 : GDALTranslateOptionsFree(translateOptions);
2762 :
2763 6 : if (ds != nullptr)
2764 : {
2765 6 : poDS = ds.get();
2766 6 : return true;
2767 : }
2768 : else
2769 : {
2770 0 : return false;
2771 : }
2772 : }
2773 310 : return true;
2774 : };
2775 :
2776 158 : if (!convertToRGBAifNeeded(poSrcDS, m_poTmpSrcDS))
2777 : {
2778 0 : ReportError(CE_Failure, CPLE_AppDefined,
2779 : "Conversion of source dataset color table to RGBA failed");
2780 0 : return false;
2781 : }
2782 :
2783 158 : if (!convertToRGBAifNeeded(poOverlayDS, m_poTmpOverlayDS))
2784 : {
2785 0 : ReportError(CE_Failure, CPLE_AppDefined,
2786 : "Conversion of overlay dataset color table to RGBA failed");
2787 0 : return false;
2788 : }
2789 :
2790 158 : if (!ValidateGlobal())
2791 0 : return false;
2792 :
2793 158 : const int nOpacity255Scale =
2794 158 : (m_opacity * 255 + OPACITY_INPUT_RANGE / 2) / OPACITY_INPUT_RANGE;
2795 :
2796 158 : bool bSwappedOpacity = false;
2797 : // Many algorithms are commutative regarding the two inputs but BlendDataset assume
2798 : // RGB(A) is in the source (and not in the overlay).
2799 420 : if ((m_operator == CompositionMode::MULTIPLY ||
2800 104 : m_operator == CompositionMode::SCREEN ||
2801 94 : m_operator == CompositionMode::HARD_LIGHT ||
2802 323 : m_operator == CompositionMode::OVERLAY) &&
2803 73 : (poSrcDS->GetRasterCount() < poOverlayDS->GetRasterCount()))
2804 : {
2805 10 : bSwappedOpacity = true;
2806 10 : std::swap(poSrcDS, poOverlayDS);
2807 : }
2808 :
2809 158 : m_outputDataset.Set(std::make_unique<BlendDataset>(
2810 158 : *poSrcDS, *poOverlayDS, m_operator, nOpacity255Scale, bSwappedOpacity));
2811 :
2812 158 : return true;
2813 : }
2814 :
2815 : GDALRasterBlendAlgorithmStandalone::~GDALRasterBlendAlgorithmStandalone() =
2816 : default;
2817 :
2818 : //! @endcond
|
