From 2073a5ae4b13510d8d73e30a565c3543bcc69bc2 Mon Sep 17 00:00:00 2001
From: Yermalayeu Ihar <ermig@tut.by>
Date: Wed, 5 Feb 2025 17:43:15 +0300
Subject: [PATCH] +add Base implementation of function BgrToLab (part 4:
 refactoring).

---
 src/Simd/SimdBaseBgrToLab.cpp | 237 +++++++---------------------------
 src/Simd/SimdBgrToLab.h       | 106 +++++++--------
 2 files changed, 101 insertions(+), 242 deletions(-)

diff --git a/src/Simd/SimdBaseBgrToLab.cpp b/src/Simd/SimdBaseBgrToLab.cpp
index 89caa26df7..420bbd546c 100644
--- a/src/Simd/SimdBaseBgrToLab.cpp
+++ b/src/Simd/SimdBaseBgrToLab.cpp
@@ -26,217 +26,74 @@
 #include "Simd/SimdMath.h"
 #include "Simd/SimdLog.h"
 
+#define SIMD_BGR_TO_LAB_OPENCV_COMPATIBILITY
+
 namespace Simd
 {
     namespace Base
     {
-#define SIMD_BGR_TO_LAB_OPENCV_COMPATIBILITY
+        uint16_t LabGammaTab[LabGammaTabSize];
+        uint16_t LabCbrtTab[LabCbrtTabSize];
+        uint32_t LabCoeffsTab[9];
 
-        namespace LabV1
+        SIMD_INLINE float ApplyGamma(float x)
         {
-            inline double FromRaw(uint64_t raw)
-            {
-                return *((double*)&raw);
-            }
-
-            SIMD_INLINE float mullAdd(float a, float b, float c)
-            {
-                return a * b + c;
-            }
-
-#define  CV_DESCALE(x,n)     (((x) + (1 << ((n)-1))) >> (n))
-            const int xyz_shift = 12;
-
-            //const double D65[3] = { 0.950456, 1., 1.088754 };
-            const double D65[3] = { FromRaw(0x3fee6a22b3892ee8),
-                                 1.0 ,
-                                 FromRaw(0x3ff16b8950763a19) };
-
-            //static const double sRGB2XYZ_D65[] =
-            //{
-            //    0.412453, 0.357580, 0.180423,
-            //    0.212671, 0.715160, 0.072169,
-            //    0.019334, 0.119193, 0.950227
-            //};
-
-            const double sRGB2XYZ_D65[] =
-            {
-                FromRaw(0x3fda65a14488c60d),
-                FromRaw(0x3fd6e297396d0918),
-                FromRaw(0x3fc71819d2391d58),
-                FromRaw(0x3fcb38cda6e75ff6),
-                FromRaw(0x3fe6e297396d0918),
-                FromRaw(0x3fb279aae6c8f755),
-                FromRaw(0x3f93cc4ac6cdaf4b),
-                FromRaw(0x3fbe836eb4e98138),
-                FromRaw(0x3fee68427418d691)
-            };
-
-            enum { LAB_CBRT_TAB_SIZE = 1024, GAMMA_TAB_SIZE = 1024 };
-
-            static const float GammaTabScale((int)GAMMA_TAB_SIZE);
-
-            static uint16_t sRGBGammaTab_b[256];
-
-#undef lab_shift
-#define lab_shift xyz_shift
-#define gamma_shift 3
-#define lab_shift2 (lab_shift + gamma_shift)
-#define LAB_CBRT_TAB_SIZE_B (256*3/2*(1<<gamma_shift))
-            static uint16_t LabCbrtTab_b[LAB_CBRT_TAB_SIZE_B];
-
-//all constants should be presented through integers to keep bit-exactness
-            static const double gammaThreshold = double(809) / double(20000);    //  0.04045
-            static const double gammaInvThreshold = double(7827) / double(2500000); //  0.0031308
-            static const double gammaLowScale = double(323) / double(25);       // 12.92
-            static const double gammaPower = double(12) / double(5);         //  2.4
-            static const double gammaXshift = double(11) / double(200);       //  0.055
+            return x <= 0.04045f ? x * (1.f / 12.92f) : (float)std::pow((double)(x + 0.055) * (1. / 1.055), 2.4);
+        }
 
-            static const float lthresh = float(216) / float(24389); // 0.008856f = (6/29)^3
-            static const float lscale = float(841) / float(108); // 7.787f = (29/3)^3/(29*4)
-            static const float lbias = float(16) / float(116);
+        SIMD_INLINE float ApplyCbrt(float x)
+        {
+            static const float threshold = float(216) / float(24389);
+            static const float scale = float(841) / float(108);
+            static const float bias = float(16) / float(116);
+            return x < threshold ? (x * scale + bias) : cbrt(x);
+        }
 
-            static inline float applyGamma(float x)
+        static bool LabTabsInited = false;
+        void LabInitTabs()
+        {
+            if (LabTabsInited)
+                return;
+            const float intScale(255 * (1 << LAB_GAMMA_SHIFT));
+            for (int i = 0; i < LabGammaTabSize; i++)
             {
-                //return x <= 0.04045f ? x*(1.f/12.92f) : (float)std::pow((double)(x + 0.055)*(1./1.055), 2.4);
-
-                double xd = x;
-                return (xd <= gammaThreshold ?
-                    xd / gammaLowScale :
-                    pow((xd + gammaXshift) / (1.0 + gammaXshift), gammaPower));
+                float x = float(i) / 255.0f;
+                LabGammaTab[i] = (uint16_t)(Round(intScale * ApplyGamma(x)));
             }
-
-            static void initLabTabs()
+            const float tabScale(1.0 / (255.0 * (1 << LAB_GAMMA_SHIFT)));
+            const float shift2(1 << LAB_SHIFT2);
+            for (int i = 0; i < LabCbrtTabSize; i++)
             {
-                static bool initialized = false;
-                if (!initialized)
-                {
-                    float scale = 1.0f / float(GammaTabScale);
-                    static const float intScale(255 * (1 << gamma_shift));
-                    for (int i = 0; i < 256; i++)
-                    {
-                        float x = float(i) / 255.0f;
-                        sRGBGammaTab_b[i] = (uint16_t)(Round(intScale * applyGamma(x)));
-                    }
-
-                    static const float cbTabScale(1.0 / (255.0 * (1 << gamma_shift)));
-                    static const float lshift2(1 << lab_shift2);
-                    for (int i = 0; i < LAB_CBRT_TAB_SIZE_B; i++)
-                    {
-                        float x = cbTabScale * float(i);
-                        LabCbrtTab_b[i] = (uint16_t)(Round(lshift2 * (x < lthresh ? mullAdd(x,  lscale,  lbias) : cbrt(x))));
-                    }
-#if defined(SIMD_BGR_TO_LAB_OPENCV_COMPATIBILITY)
-                    if(LabCbrtTab_b[324] == 17746)
-                        LabCbrtTab_b[324] = 17745;
-                    if (LabCbrtTab_b[49] == 9455)
-                        LabCbrtTab_b[49] = 9454;
-#endif
-                    initialized = true;
-                }
+                float x = tabScale * float(i);
+                LabCbrtTab[i] = (uint16_t)Round(shift2 * ApplyCbrt(x));
             }
-
-            struct RGB2Lab_b
-            {
-                typedef uint8_t channel_type;
-
-                RGB2Lab_b(int _srccn, int blueIdx)
-                    : srccn(_srccn)
-                {
-                    initLabTabs();
-
-                    double whitePt[3];
-                    for (int i = 0; i < 3; i++)
-                        whitePt[i] = D65[i];
-
-                    const double lshift(1 << lab_shift);
-                    for (int i = 0; i < 3; i++)
-                    {
-                        double c[3];
-                        for (int j = 0; j < 3; j++)
-                            c[j] = sRGB2XYZ_D65[i * 3 + j];
-                        coeffs[i * 3 + (blueIdx ^ 2)] = Round(lshift * c[0] / whitePt[i]);
-                        coeffs[i * 3 + 1] = Round(lshift * c[1] / whitePt[i]);
-                        coeffs[i * 3 + blueIdx] = Round(lshift * c[2] / whitePt[i]);
-
-                        assert(coeffs[i * 3] >= 0 && coeffs[i * 3 + 1] >= 0 && coeffs[i * 3 + 2] >= 0 &&
-                            coeffs[i * 3] + coeffs[i * 3 + 1] + coeffs[i * 3 + 2] < 2 * (1 << lab_shift));
-                    }
-                }
-
-                void operator()(const uint8_t* src, uint8_t* dst, int n) const
-                {
-                    const int Lscale = (116 * 255 + 50) / 100;
-                    const int Lshift = -((16 * 255 * (1 << lab_shift2) + 50) / 100);
-                    const uint16_t* tab = sRGBGammaTab_b;
-                    int i, scn = srccn;
-                    int C0 = coeffs[0], C1 = coeffs[1], C2 = coeffs[2],
-                        C3 = coeffs[3], C4 = coeffs[4], C5 = coeffs[5],
-                        C6 = coeffs[6], C7 = coeffs[7], C8 = coeffs[8];
-
-                    for (int i = 0; i < n; i++, src += scn, dst += 3)
-                    {
-                        int R = tab[src[0]], G = tab[src[1]], B = tab[src[2]];
-                        int fX = LabCbrtTab_b[CV_DESCALE(R * C0 + G * C1 + B * C2, lab_shift)];
-                        int fY = LabCbrtTab_b[CV_DESCALE(R * C3 + G * C4 + B * C5, lab_shift)];
-                        int fZ = LabCbrtTab_b[CV_DESCALE(R * C6 + G * C7 + B * C8, lab_shift)];
-
-                        int L = CV_DESCALE(Lscale * fY + Lshift, lab_shift2);
-                        int a = CV_DESCALE(500 * (fX - fY) + 128 * (1 << lab_shift2), lab_shift2);
-                        int b = CV_DESCALE(200 * (fY - fZ) + 128 * (1 << lab_shift2), lab_shift2);
-
-#if defined(_WIN32) && 0
-                        if (L == 45 && a == 165)
-                        {
-                            std::cout << i << " old Lab = {" << L << ", " << a << ", " << b << "}";
-
-                            int R = tab[src[0]], G = tab[src[1]], B = tab[src[2]];
-                            int fX = LabCbrtTab_b[CV_DESCALE(R * C0 + G * C1 + B * C2, lab_shift)];
-                            int fY = LabCbrtTab_b[CV_DESCALE(R * C3 + G * C4 + B * C5, lab_shift)] - 1;
-                            int fZ = LabCbrtTab_b[CV_DESCALE(R * C6 + G * C7 + B * C8, lab_shift)];
-
-                            int L = CV_DESCALE(Lscale * fY + Lshift, lab_shift2);
-                            int a = CV_DESCALE(500 * (fX - fY) + 128 * (1 << lab_shift2), lab_shift2);
-                            int b = CV_DESCALE(200 * (fY - fZ) + 128 * (1 << lab_shift2), lab_shift2);
-
-                            std::cout << " new Lab = {" << L << ", " << a << ", " << b << "}";
-                            std::cout << " index " << CV_DESCALE(R * C3 + G * C4 + B * C5, lab_shift);
-                            std::cout << " value " << LabCbrtTab_b[CV_DESCALE(R * C3 + G * C4 + B * C5, lab_shift)];
-                            std::cout << std::endl;
-                            exit(0);
-                        }
+#if defined(SIMD_BGR_TO_LAB_OPENCV_COMPATIBILITY)
+            if (LabCbrtTab[324] == 17746)
+                LabCbrtTab[324] = 17745;
+            if (LabCbrtTab[49] == 9455)
+                LabCbrtTab[49] = 9454;
 #endif
-                        dst[0] = Base::RestrictRange(L);
-                        dst[1] = Base::RestrictRange(a);
-                        dst[2] = Base::RestrictRange(b);
-                    }
-                }
-
-                int srccn;
-                int coeffs[9];
-            };
+            const double D65[3] = { 0.950456, 1., 1.088754 };
+            const double sRGB2XYZ_D65[9] = { 0.412453, 0.357580, 0.180423, 0.212671, 0.715160, 0.072169, 0.019334, 0.119193, 0.950227 };
+            const double shift(1 << LAB_SHIFT);
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    LabCoeffsTab[i * 3 + j] = Round(shift * sRGB2XYZ_D65[i * 3 + j] / D65[i]);
+            LabTabsInited = true;
         }
 
+        //--------------------------------------------------------------------------------------------------
+
         void BgrToLab(const uint8_t* bgr, size_t bgrStride, size_t width, size_t height, uint8_t* lab, size_t labStride)
         {
-#if 1
-            LabV1::RGB2Lab_b _lab(3, 0);
+            LabInitTabs();
             for (size_t row = 0; row < height; ++row)
             {
-                const uint8_t* pBgr = bgr + row * bgrStride;
+                const uint8_t* pBgr = bgr + row * bgrStride, * pEnd = pBgr + width * 3;
                 uint8_t* pLab = lab + row * labStride;
-                _lab(pBgr, pLab, width);
+                for (; pBgr < pEnd; pBgr += 3, pLab += 3)
+                    RgbToLab(pBgr[2], pBgr[1], pBgr[0], pLab);
             }
-
-#else
-            for (size_t row = 0; row < height; ++row)
-            {
-                const uint8_t* pBgr = bgr + row * bgrStride;
-                uint8_t* pLab = lab + row * labStride;
-                for (const uint8_t* pBgrEnd = pBgr + width * 3; pBgr < pBgrEnd; pBgr += 3, pLab += 3)
-                    BgrToLab(pBgr[0], pBgr[1], pBgr[2], pLab);
-            }
-#endif
         }
     }
 }
diff --git a/src/Simd/SimdBgrToLab.h b/src/Simd/SimdBgrToLab.h
index 7d2f913078..75572ac60d 100644
--- a/src/Simd/SimdBgrToLab.h
+++ b/src/Simd/SimdBgrToLab.h
@@ -24,67 +24,69 @@
 #ifndef __SimdBgrToLab_h__
 #define __SimdBgrToLab_h__
 
-#include "Simd/SimdDefs.h"
+#include "Simd/SimdMath.h"
 
 namespace Simd
 {
     namespace Base
     {
-        namespace LabV0
+        const int LAB_GAMMA_SHIFT = 3;
+
+        const int LabGammaTabSize = 256;
+        extern uint16_t LabGammaTab[LabGammaTabSize];
+
+        const int LabCbrtTabSize = 256 * 3 / 2 * (1 << LAB_GAMMA_SHIFT);
+        extern uint16_t LabCbrtTab[LabCbrtTabSize];
+
+        const int LabCoeffsTabSize = 9;
+        extern uint32_t LabCoeffsTab[LabCoeffsTabSize];
+
+        void LabInitTabs();
+
+        const int LAB_SHIFT = 12;
+        const int LAB_ROUND = 1 << (LAB_SHIFT - 1);
+
+        const int LAB_SHIFT2 = LAB_SHIFT + LAB_GAMMA_SHIFT;
+        const int LAB_ROUND2 = 1 << (LAB_SHIFT2 - 1);
+
+        const int LAB_L_SCALE = (116 * 255 + 50) / 100;
+        const int LAB_L_SHIFT = -((16 * 255 * (1 << LAB_SHIFT2) + 50) / 100);
+
+        const int LAB_A_SCALE = 500;
+        const int LAB_B_SCALE = 200;
+        const int LAB_AB_SHIFT = 128 * (1 << LAB_SHIFT2);
+
+        SIMD_INLINE int LabDescale(int value)
+        {
+            return (value + LAB_ROUND) >> LAB_SHIFT;
+        }
+
+        SIMD_INLINE int LabDescale2(int value)
         {
-#define X(xyz) xyz.data[0]
-#define Y(xyz) xyz.data[1]
-#define Z(xyz) xyz.data[2]
-
-            typedef union
-            {
-                double data[3];
-                struct { double r, g, b; };
-                struct { double L, A, B; };
-                struct { double l, c, h; };
-            } DoubleTriplet;
-
-            const double eps = (6 * 6 * 6) / (29.0 * 29.0 * 29.0), kap = (29 * 29 * 29) / (3.0 * 3.0 * 3.0);
-
-            const DoubleTriplet xyzReferenceValues = { {0.95047, 1.0, 1.08883} };
-
-            SIMD_INLINE DoubleTriplet xyzFromRgb(DoubleTriplet rgb)
-            {
-                for (int i = 0; i < 3; ++i)
-                {
-                    double v = rgb.data[i];
-                    rgb.data[i] = (v > 0.04045) ? pow(((v + 0.055) / 1.055), 2.4) : (v / 12.92);
-                }
-                DoubleTriplet temp = { {
-                    rgb.r * 0.4124564 + rgb.g * 0.3575761 + rgb.b * 0.1804375,
-                    rgb.r * 0.2126729 + rgb.g * 0.7151522 + rgb.b * 0.0721750,
-                    rgb.r * 0.0193339 + rgb.g * 0.1191920 + rgb.b * 0.9503041
-                    } };
-                return temp;
-            }
-
-            SIMD_INLINE DoubleTriplet labFromXyz(DoubleTriplet xyz)
-            {
-                for (int i = 0; i < 3; ++i)
-                {
-                    xyz.data[i] /= xyzReferenceValues.data[i];
-                    double v = xyz.data[i];
-                    xyz.data[i] = (v > eps) ? pow(v, (1 / 3.0)) : ((kap * v + 16) / 116.0);
-                }
-                DoubleTriplet temp = { {(116 * Y(xyz)) - 16, 500 * (X(xyz) - Y(xyz)), 200 * (Y(xyz) - Z(xyz))} };
-                return temp;
-            }
-
-            DoubleTriplet labFromRgb(DoubleTriplet rgb) { return labFromXyz(xyzFromRgb(rgb)); }
+            return (value + LAB_ROUND2) >> LAB_SHIFT2;
         }
 
-        SIMD_INLINE void BgrToLab(int blue, int green, int red, uint8_t* lab)
+        SIMD_INLINE void RgbToLab(int red, int green, int blue, uint8_t *lab)
         {
-            LabV0::DoubleTriplet _rgb{red / 255.0 , green / 255.0, blue / 255.0 };
-            LabV0::DoubleTriplet _lab = labFromRgb(_rgb);
-            lab[0] = int(_lab.data[0]);
-            lab[1] = int(_lab.data[1]);
-            lab[2] = int(_lab.data[2]);
+            int R = LabGammaTab[red];
+            int G = LabGammaTab[green];
+            int B = LabGammaTab[blue];
+
+            int iX = LabDescale(R * LabCoeffsTab[0] + G * LabCoeffsTab[1] + B * LabCoeffsTab[2]);
+            int iY = LabDescale(R * LabCoeffsTab[3] + G * LabCoeffsTab[4] + B * LabCoeffsTab[5]);
+            int iZ = LabDescale(R * LabCoeffsTab[6] + G * LabCoeffsTab[7] + B * LabCoeffsTab[8]);
+
+            int fX = LabCbrtTab[iX];
+            int fY = LabCbrtTab[iY];
+            int fZ = LabCbrtTab[iZ];
+
+            int L = LabDescale2(LAB_L_SCALE * fY + LAB_L_SHIFT);
+            int a = LabDescale2(LAB_A_SCALE * (fX - fY) + LAB_AB_SHIFT);
+            int b = LabDescale2(LAB_B_SCALE * (fY - fZ) + LAB_AB_SHIFT);
+
+            lab[0] = Base::RestrictRange(L);
+            lab[1] = Base::RestrictRange(a);
+            lab[2] = Base::RestrictRange(b);
         }
     }
 }