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