skinpeel

icbc.h (127319B)

---

 // icbc.h v1.05
 // A High Quality SIMD BC1 Encoder by Ignacio Castano <castano@gmail.com>.
 //
 // LICENSE:
 //  MIT license at the end of this file.
 
 #ifndef ICBC_H
 #define ICBC_H
 
 namespace icbc {
 
     enum Decoder {
         Decoder_D3D10 = 0,
         Decoder_NVIDIA = 1,
         Decoder_AMD = 2,    // Apple's M1 decoder appears to match AMD's
         Decoder_Intel = 3
     };
 
     void init(Decoder decoder = Decoder_D3D10);
 
     enum Quality {
         Quality_Level0,  // Box fit only.
         Quality_Level1,  // Box fit + least squares fit.
         Quality_Level2,  // Cluster fit 4, threshold = 24.
         Quality_Level3,  // Cluster fit 4, threshold = 32.
         Quality_Level4,  // Cluster fit 4, threshold = 48.
         Quality_Level5,  // Cluster fit 4, threshold = 64.
         Quality_Level6,  // Cluster fit 4, threshold = 96.
         Quality_Level7,  // Cluster fit 4, threshold = 128.
         Quality_Level8,  // Cluster fit 4+3, threshold = 256.
         Quality_Level9,  // Cluster fit 4+3, threshold = 256 + Refinement.
 
         Quality_Fast = Quality_Level1,
         Quality_Default = Quality_Level8,
         Quality_Max = Quality_Level9,
     };
 
     void decode_bc1(const void * block, unsigned char rgba_block[16 * 4], Decoder decoder = Decoder_D3D10);
     void decode_bc3(const void * block, unsigned char rgba_block[16 * 4], Decoder decoder = Decoder_D3D10);
     //void decode_bc4(const void * block, bool snorm, unsigned char rgba_block[16 * 4], Decoder decoder = Decoder_D3D10);
     //void decode_bc5(const void * block, bool snorm, unsigned char rgba_block[16 * 4], Decoder decoder = Decoder_D3D10);
 
     float evaluate_bc1_error(const unsigned char rgba_block[16 * 4], const void * block, Decoder decoder = Decoder_D3D10);
     float evaluate_bc3_error(const unsigned char rgba_block[16 * 4], const void * block, bool alpha_blend, Decoder decoder = Decoder_D3D10);
     //float evaluate_bc4_error(const float rgba_block[16 * 4], const void * block, bool snorm, Decoder decoder = Decoder_D3D10);
     //float evaluate_bc5_error(const float rgba_block[16 * 4], const void * block, bool snorm, Decoder decoder = Decoder_D3D10);
 
     float compress_bc1(Quality level, const float * input_colors, const float * input_weights, const float color_weights[3], bool three_color_mode, bool three_color_black, void * output);
     float compress_bc3(Quality level, const float * input_colors, const float * input_weights, const float color_weights[3], bool six_alpha_mode, void * output);
     //float compress_bc4(Quality level, const float * input_colors, const float * input_weights, bool snorm, bool six_alpha_mode, void * output);
     //float compress_bc5(Quality level, const float * input_colors, const float * input_weights, bool snorm, bool six_alpha_mode, void * output);
 }
 
 #endif // ICBC_H
 
 #ifdef ICBC_IMPLEMENTATION
 
 // Instruction level support must be chosen at compile time setting ICBC_SIMD to one of these values:
 #define ICBC_SCALAR 0
 #define ICBC_SSE2   1
 #define ICBC_SSE41  2
 #define ICBC_AVX1   3
 #define ICBC_AVX2   4
 #define ICBC_AVX512 5
 #define ICBC_NEON   -1
 #define ICBC_VMX    -2
 
 #if defined(__i386__) || defined(_M_IX86) || defined(__x86_64__) || defined(_M_X64)
     #define ICBC_X86 1
 #endif
 
 #if defined(__x86_64__) || defined(_M_X64)
     #define ICBC_X64 1
 #endif
 
 #if defined(__arm__) || defined(_M_ARM) || defined(__aarch64__) || defined(_M_ARM64)
     #define ICBC_ARM 1
 #endif
 
 #if defined(__aarch64__) || defined(_M_ARM64)
     #define ICBC_ARM64 1
 #endif
 
 #if (defined(__PPC__) || defined(_M_PPC))
     #define ICBC_PPC 1
 #endif
 
 // SIMD version.
 #ifndef ICBC_SIMD
     #if ICBC_X86
         #if __AVX512F__
             #define ICBC_SIMD ICBC_AVX512
         #elif __AVX2__
             #define ICBC_SIMD ICBC_AVX2
         #elif __AVX__
             #define ICBC_SIMD ICBC_AVX1
         #elif __SSE4_1__
             #define ICBC_SIMD ICBC_SSE41
         #elif __SSE2__ || ICBC_X64
             #define ICBC_SIMD ICBC_SSE2
         #else
             #define ICBC_SIMD ICBC_SCALAR
         #endif
     #endif
 
     #if ICBC_ARM
         #if __ARM_NEON__
             #define ICBC_SIMD ICBC_NEON
         #else
             #define ICBC_SIMD ICBC_SCALAR
         #endif
     #endif
 
     #if ICBC_PPC
         #define ICBC_SIMD ICBC_VMX
     #endif
 #endif
 
 // AVX1 does not require FMA, and depending on whether it's Intel or AMD you may have FMA3 or FMA4. What a mess.
 #ifndef ICBC_USE_FMA
 //#define ICBC_USE_FMA 3
 //#define ICBC_USE_FMA 4
 #endif
 
 #if ICBC_SIMD >= ICBC_AVX2
 #define ICBC_BMI2 1
 #endif
 
 // Apparently rcp is not deterministic (different precision on Intel and AMD), enable if you don't care about that for a small performance boost.
 //#define ICBC_USE_RCP 1
 
 #if ICBC_SIMD == ICBC_AVX2
 #define ICBC_USE_AVX2_PERMUTE2 1    // Using permutevar8x32 and bitops.
 #endif
 
 #if ICBC_SIMD == ICBC_AVX512
 #define ICBC_USE_AVX512_PERMUTE 1
 #endif
 
 #if ICBC_SIMD == ICBC_NEON
 #define ICBC_USE_NEON_VTL 1         // Enable this for a 0.89x performance improvement.
 #endif
 
 #if ICBC_SIMD == ICBC_SSE41
 #define ICBC_USE_SSSE3_SHUFFLEB 1   // Enable this for 0.88x performance improvement.
 #endif
 
 
 // Some experimental knobs:
 #define ICBC_PERFECT_ROUND 0        // Enable perfect rounding to compute cluster fit residual.
 
 
 #if ICBC_SIMD >= ICBC_SSE2
 #include <emmintrin.h>
 #endif
 
 #if ICBC_SIMD >= ICBC_SSE41
 #include <smmintrin.h>
 #endif
 
 #if ICBC_SIMD >= ICBC_AVX1
 #include <immintrin.h>
 #endif
 
 #if ICBC_SIMD >= ICBC_AVX512 && _MSC_VER
 #include <zmmintrin.h>
 #endif
 
 #if ICBC_SIMD == ICBC_NEON
 #include <arm_neon.h>
 #endif
 
 #if ICBC_SIMD == ICBC_VMX
 #include <altivec.h>
 #endif
 
 #if _MSC_VER
 #include <intrin.h> // _BitScanReverse
 #endif
 
 #include <stdint.h>
 #include <stdlib.h> // abs
 #include <string.h> // memset
 #include <math.h>   // fabsf
 #include <float.h>  // FLT_MAX
 
 #ifndef ICBC_ASSERT
 #if _DEBUG
 #define ICBC_ASSERT assert
 #include <assert.h>
 #else
 #define ICBC_ASSERT(x)
 #endif
 #endif
 
 namespace icbc {
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Basic Templates
 
 template <typename T> inline void swap(T & a, T & b) {
     T temp(a);
     a = b;
     b = temp;
 }
 
 template <typename T> inline T max(const T & a, const T & b) {
     return (b < a) ? a : b;
 }
 
 template <typename T> inline T min(const T & a, const T & b) {
     return (a < b) ? a : b;
 }
 
 template <typename T> inline T clamp(const T & x, const T & a, const T & b) {
     return min(max(x, a), b);
 }
 
 template <typename T> inline T square(const T & a) {
     return a * a;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Basic Types
 
 typedef uint8_t uint8;
 typedef int8_t int8;
 typedef uint16_t uint16;
 typedef uint32_t uint32;
 typedef uint32_t uint;
 
 
 struct Color16 {
     union {
         struct {
             uint16 b : 5;
             uint16 g : 6;
             uint16 r : 5;
         };
         uint16 u;
     };
 };
 
 struct Color32 {
     union {
         struct {
             uint8 b, g, r, a;
         };
         uint32 u;
     };
 };
 
 struct BlockBC1 {
     Color16 col0;
     Color16 col1;
     uint32 indices;
 };
 
 struct BlockBC4 {
     uint8 alpha0;
     uint8 alpha1;
     uint8 indices[6];
 };
 
 struct BlockBC3 {
     BlockBC4 alpha;
     BlockBC1 rgb;
 };
 
 struct BlockBC5 {
     BlockBC4 x;
     BlockBC1 y;
 };
 
 
 
 struct Vector3 {
     float x;
     float y;
     float z;
 
     inline void operator+=(Vector3 v) {
         x += v.x; y += v.y; z += v.z;
     }
     inline void operator*=(Vector3 v) {
         x *= v.x; y *= v.y; z *= v.z;
     }
     inline void operator*=(float s) {
         x *= s; y *= s; z *= s;
     }
 };
 
 struct Vector4 {
     union {
         struct {
             float x, y, z, w;
         };
         Vector3 xyz;
     };
 };
 
 
 inline Vector3 operator*(Vector3 v, float s) {
     return { v.x * s, v.y * s, v.z * s };
 }
 
 inline Vector3 operator*(float s, Vector3 v) {
     return { v.x * s, v.y * s, v.z * s };
 }
 
 inline Vector3 operator*(Vector3 a, Vector3 b) {
     return { a.x * b.x, a.y * b.y, a.z * b.z };
 }
 
 inline float dot(Vector3 a, Vector3 b) {
     return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
 inline Vector3 operator+(Vector3 a, Vector3 b) {
     return { a.x + b.x, a.y + b.y, a.z + b.z };
 }
 
 inline Vector3 operator-(Vector3 a, Vector3 b) {
     return { a.x - b.x, a.y - b.y, a.z - b.z };
 }
 
 inline Vector3 operator/(Vector3 v, float s) {
     return { v.x / s, v.y / s, v.z / s };
 }
 
 inline float saturate(float x) {
     return clamp(x, 0.0f, 1.0f);
 }
 
 inline Vector3 saturate(Vector3 v) {
     return { saturate(v.x), saturate(v.y), saturate(v.z) };
 }
 
 inline Vector3 min(Vector3 a, Vector3 b) {
     return { min(a.x, b.x), min(a.y, b.y), min(a.z, b.z) };
 }
 
 inline Vector3 max(Vector3 a, Vector3 b) {
     return { max(a.x, b.x), max(a.y, b.y), max(a.z, b.z) };
 }
 
 inline bool operator==(const Vector3 & a, const Vector3 & b) {
     return a.x == b.x && a.y == b.y && a.z == b.z;
 }
 
 inline Vector3 scalar_to_vector3(float f) {
     return {f, f, f};
 }
 
 inline float lengthSquared(Vector3 v) {
     return dot(v, v);
 }
 
 inline bool equal(float a, float b, float epsilon = 0.0001) {
     // http://realtimecollisiondetection.net/blog/?p=89
     //return fabsf(a - b) < epsilon * max(1.0f, max(fabsf(a), fabsf(b)));
     return fabsf(a - b) < epsilon;
 }
 
 inline bool equal(Vector3 a, Vector3 b, float epsilon) {
     return equal(a.x, b.x, epsilon) && equal(a.y, b.y, epsilon) && equal(a.z, b.z, epsilon);
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // SIMD
 
 #ifndef ICBC_ALIGN
 #if __GNUC__
 #   define ICBC_ALIGN __attribute__ ((__aligned__ (icbc::VEC_SIZE*4)))
 #else // _MSC_VER
 #   define ICBC_ALIGN __declspec(align(icbc::VEC_SIZE*4))
 #endif
 #endif
 
 #if __GNUC__
 #define ICBC_FORCEINLINE inline __attribute__((always_inline))
 #else
 #define ICBC_FORCEINLINE __forceinline
 #endif
 
 
 // Count trailing zeros (BSR).
 ICBC_FORCEINLINE int ctz(uint mask) {
 #if __GNUC__
     return __builtin_ctz(mask);
 #else
     unsigned long index;
     _BitScanReverse(&index, mask);
     return (int)index;
 #endif
 }
 
 
 #if ICBC_SIMD == ICBC_SCALAR  // Purely scalar version.
 
 constexpr int VEC_SIZE = 1;
 
 using VFloat = float;
 using VMask = bool;
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) { return v; }
 ICBC_FORCEINLINE VFloat vzero() { return 0.0f; }
 ICBC_FORCEINLINE VFloat vbroadcast(float x) { return x; }
 ICBC_FORCEINLINE VFloat vload(const float * ptr) { return *ptr; }
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) { return 1.0f / a; }
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) { return a * b + c; }
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) { return a * b - c; }
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) { return a * b - c * d; }
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) { return min(max(a, 0.0f), 1.0f); }
 ICBC_FORCEINLINE VFloat vround01(VFloat a) { return float(int(a + 0.5f)); }
 ICBC_FORCEINLINE VFloat lane_id() { return 0; }
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) { return mask ? b : a; }
 ICBC_FORCEINLINE VMask vbroadcast(bool b) { return b; }
 ICBC_FORCEINLINE bool all(VMask m) { return m; }
 ICBC_FORCEINLINE bool any(VMask m) { return m; }
 ICBC_FORCEINLINE uint mask(VMask m) { return (uint)m; }
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) { return 0; }
 ICBC_FORCEINLINE void vtranspose4(VFloat & a, VFloat & b, VFloat & c, VFloat & d) {}
 
 #elif ICBC_SIMD == ICBC_SSE2 || ICBC_SIMD == ICBC_SSE41
 
 constexpr int VEC_SIZE = 4;
 
 #if __GNUC__
 // GCC needs a struct so that we can overload operators.
 union VFloat {
     __m128 v;
     float m128_f32[VEC_SIZE];
 
     VFloat() {}
     VFloat(__m128 v) : v(v) {}
     operator __m128 & () { return v; }
 };
 union VMask {
     __m128 m;
 
     VMask() {}
     VMask(__m128 m) : m(m) {}
     operator __m128 & () { return m; }
 };
 #else
 using VFloat = __m128;
 using VMask = __m128;
 #endif
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) {
     return v.m128_f32[i];
 }
 
 ICBC_FORCEINLINE VFloat vzero() {
     return _mm_setzero_ps();
 }
 
 ICBC_FORCEINLINE VFloat vbroadcast(float x) {
     return _mm_set1_ps(x);
 }
 
 ICBC_FORCEINLINE VFloat vload(const float * ptr) {
     return _mm_load_ps(ptr);
 }
 
 ICBC_FORCEINLINE VFloat vgather(const float * base, VFloat index) {
     VFloat v;
     for (int i = 0; i < VEC_SIZE; i++) {
         lane(v, i) = base[int(lane(index, i))];
     }
     return v;
 }
 
 ICBC_FORCEINLINE VFloat operator+(VFloat a, VFloat b) {
     return _mm_add_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator-(VFloat a, VFloat b) {
     return _mm_sub_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator*(VFloat a, VFloat b) {
     return _mm_mul_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) {
 #if ICBC_USE_RCP
     VFloat r = _mm_rcp_ps(a);
     return _mm_mul_ps(r, _mm_sub_ps(vbroadcast(2.0f), _mm_mul_ps(r, a)));   // r * (2 - r * a)
 #else
     return _mm_div_ps(vbroadcast(1.0f), a);
 #endif
 }
 
 // a*b+c
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) {
     return a * b + c;
 }
 
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) {
     return a * b - c;
 }
 
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) {
     return a * b - c * d;
 }
 
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) {
     auto zero = _mm_setzero_ps();
     auto one = _mm_set1_ps(1.0f);
     return _mm_min_ps(_mm_max_ps(a, zero), one);
 }
 
 // Assumes a is in [0, 1] range.
 ICBC_FORCEINLINE VFloat vround01(VFloat a) {
 #if ICBC_SIMD == ICBC_SSE41
     return _mm_round_ps(a, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC);
 #else
     return _mm_cvtepi32_ps(_mm_cvttps_epi32(a + vbroadcast(0.5f)));
 #endif
 }
 
 ICBC_FORCEINLINE VFloat vtruncate(VFloat a) {
 #if ICBC_SIMD == ICBC_SSE41
     return _mm_round_ps(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
 #else
     return _mm_cvtepi32_ps(_mm_cvttps_epi32(a));
 #endif
 }
 
 ICBC_FORCEINLINE VFloat lane_id() {
     return _mm_set_ps(3, 2, 1, 0);
 }
 
 ICBC_FORCEINLINE VMask operator> (VFloat A, VFloat B) { return _mm_cmpgt_ps(A, B); }
 ICBC_FORCEINLINE VMask operator>=(VFloat A, VFloat B) { return _mm_cmpge_ps(A, B); }
 ICBC_FORCEINLINE VMask operator< (VFloat A, VFloat B) { return _mm_cmplt_ps(A, B); }
 ICBC_FORCEINLINE VMask operator<=(VFloat A, VFloat B) { return _mm_cmple_ps(A, B); }
 
 ICBC_FORCEINLINE VMask operator| (VMask A, VMask B) { return _mm_or_ps(A, B); }
 ICBC_FORCEINLINE VMask operator& (VMask A, VMask B) { return _mm_and_ps(A, B); }
 ICBC_FORCEINLINE VMask operator^ (VMask A, VMask B) { return _mm_xor_ps(A, B); }
 
 // mask ? b : a
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) {
 #if ICBC_SIMD == ICBC_SSE41
     return _mm_blendv_ps(a, b, mask);
 #else
     return _mm_or_ps(_mm_andnot_ps(mask, a), _mm_and_ps(mask, b));
 #endif
 }
 
 ICBC_FORCEINLINE VMask vbroadcast(bool b) { 
     return _mm_castsi128_ps(_mm_set1_epi32(-int32_t(b)));
 }
 
 ICBC_FORCEINLINE bool all(VMask m) {
     int value = _mm_movemask_ps(m);
     return value == 0x7;
 }
 
 ICBC_FORCEINLINE bool any(VMask m) {
     int value = _mm_movemask_ps(m);
     return value != 0;
 }
 
 ICBC_FORCEINLINE uint mask(VMask m) {
     return (uint)_mm_movemask_ps(m);
 }
 
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) {
 
     // First do an horizontal reduction.                            // v = [ D C | B A ]
     VFloat shuf = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 0, 1));    //     [ C D | A B ]
     VFloat mins = _mm_min_ps(v, shuf);                              // mins = [ D+C C+D | B+A A+B ]
     shuf        = _mm_movehl_ps(shuf, mins);                        //        [   C   D | D+C C+D ]  // let the compiler avoid a mov by reusing shuf
     mins        = _mm_min_ss(mins, shuf);
     mins =      _mm_shuffle_ps(mins, mins, _MM_SHUFFLE(0, 0, 0, 0));
 
     // Then find the index.
     uint mask = _mm_movemask_ps(v <= mins);
     return ctz(mask);
 }
 
 // https://gcc.gnu.org/legacy-ml/gcc-patches/2005-10/msg00324.html
 ICBC_FORCEINLINE void vtranspose4(VFloat & r0, VFloat & r1, VFloat & r2, VFloat & r3) {
     VFloat t0 = _mm_unpacklo_ps(r0, r1);
     VFloat t1 = _mm_unpacklo_ps(r2, r3);
     VFloat t2 = _mm_unpackhi_ps(r0, r1);
     VFloat t3 = _mm_unpackhi_ps(r2, r3);
     r0 = _mm_movelh_ps(t0, t1);
     r1 = _mm_movehl_ps(t1, t0);
     r2 = _mm_movelh_ps(t2, t3);
     r3 = _mm_movehl_ps(t3, t2);
 }
 
 
 #if ICBC_SIMD == ICBC_SSE41
 
 ICBC_FORCEINLINE VFloat vpermute2(VFloat tab0, VFloat tab1, __m128i idx) {
 
     // @@ Precompute this:
     tab1 = _mm_xor_ps(tab1, tab0);
 
     __m128i result = _mm_shuffle_epi8(_mm_castps_si128(tab0), idx);
     idx = _mm_sub_epi8(idx, _mm_set1_epi8(16));
 
     result = _mm_xor_si128(result, _mm_shuffle_epi8(_mm_castps_si128(tab1), idx));
 
     return _mm_castsi128_ps(result);
 }
 
 ICBC_FORCEINLINE VFloat vpermute4(VFloat tab0, VFloat tab1, VFloat tab2, VFloat tab3, __m128i idx) {
 
     // @@ Precompute this:
     tab3 = _mm_xor_ps(tab3, tab2);
     tab2 = _mm_xor_ps(tab2, tab1);
     tab1 = _mm_xor_ps(tab1, tab0);
 
     __m128i result = _mm_shuffle_epi8(_mm_castps_si128(tab0), idx);
     idx = _mm_sub_epi8(idx, _mm_set1_epi8(16));
 
     result = _mm_xor_si128(result, _mm_shuffle_epi8(_mm_castps_si128(tab1), idx));
     idx = _mm_sub_epi8(idx, _mm_set1_epi8(16));
 
     result = _mm_xor_si128(result, _mm_shuffle_epi8(_mm_castps_si128(tab2), idx));
     idx = _mm_sub_epi8(idx, _mm_set1_epi8(16));
 
     result = _mm_xor_si128(result, _mm_shuffle_epi8(_mm_castps_si128(tab3), idx));
 
     return _mm_castsi128_ps(result);
 }
 
 #endif
 
 
 #elif ICBC_SIMD == ICBC_AVX1 || ICBC_SIMD == ICBC_AVX2
 
 constexpr int VEC_SIZE = 8;
 
 #if __GNUC__
 union VFloat {
     __m256 v;
     float m256_f32[VEC_SIZE];
 
     VFloat() {}
     VFloat(__m256 v) : v(v) {}
     operator __m256 & () { return v; }
 };
 union VInt {
     __m256i v;
     int m256_i32[VEC_SIZE];
 
     VInt() {}
     VInt(__m256i v) : v(v) {}
     operator __m256i & () { return v; }
 };
 union VMask {
     __m256 m;
 
     VMask() {}
     VMask(__m256 m) : m(m) {}
     operator __m256 & () { return m; }
 };
 #else
 using VFloat = __m256;
 using VInt = __m256i;
 using VMask = __m256;   // Emulate mask vector using packed float.
 #endif
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) {
     return v.m256_f32[i];
 }
 
 ICBC_FORCEINLINE VFloat vzero() {
     return _mm256_setzero_ps();
 }
 
 ICBC_FORCEINLINE VFloat vbroadcast(float a) {
     return _mm256_set1_ps(a);
 }
 
 ICBC_FORCEINLINE VFloat vload(const float * ptr) {
     return _mm256_load_ps(ptr);
 }
 
 ICBC_FORCEINLINE VFloat vgather(const float * base, VFloat index) {
 #if ICBC_SIMD == ICBC_AVX2
     return _mm256_i32gather_ps(base, _mm256_cvtps_epi32(index), 4);
 #else
     VFloat v;
     for (int i = 0; i < VEC_SIZE; i++) {
         lane(v, i) = base[int(lane(index, i))];
     }
     return v;
 #endif
 }
 
 ICBC_FORCEINLINE VFloat operator+(VFloat a, VFloat b) {
     return _mm256_add_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator-(VFloat a, VFloat b) {
     return _mm256_sub_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator*(VFloat a, VFloat b) {
     return _mm256_mul_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) {
 #if ICBC_USE_RCP
     #if ICBC_SIMD == ICBC_AVX512
         VFloat r = _mm256_rcp14_ps(a);
     #else
         VFloat r = _mm256_rcp_ps(a);
     #endif
 
     // r = r * (2 - r * a)
     #if ICBC_USE_FMA == 3 || ICBC_AVX2
         return _mm256_mul_ps(r, _mm256_fnmadd_ps(r, a, vbroadcast(2.0f)));
     #else
         return _mm256_mul_ps(r, _mm256_sub_ps(vbroadcast(2.0f), _mm256_mul_ps(r, a)));
     #endif
 #else
     return _mm256_div_ps(vbroadcast(1.0f), a);
 #endif
 }
 
 // a*b+c
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) {
 #if ICBC_USE_FMA == 3 || ICBC_SIMD == ICBC_AVX2
     return _mm256_fmadd_ps(a, b, c);
 #elif ICBC_USE_FMA == 4
     return _mm256_macc_ps(a, b, c);
 #else
     return ((a * b) + c);
 #endif
 }
 
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) {
 #if ICBC_USE_FMA == 3 || ICBC_SIMD == ICBC_AVX2
     return _mm256_fmsub_ps(a, b, c);
 #elif ICBC_USE_FMA == 4
     return _mm256_msub_ps(a, b, c);
 #else
     return ((a * b) - c);
 #endif
 }
 
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) {
     return vmsub(a, b, c * d);
 }
 
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) {
     __m256 zero = _mm256_setzero_ps();
     __m256 one = _mm256_set1_ps(1.0f);
     return _mm256_min_ps(_mm256_max_ps(a, zero), one);
 }
 
 ICBC_FORCEINLINE VFloat vround01(VFloat a) {
     return _mm256_round_ps(a, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC);
 }
 
 ICBC_FORCEINLINE VFloat vtruncate(VFloat a) {
     return _mm256_round_ps(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
 }
 
 ICBC_FORCEINLINE VFloat lane_id() {
     return _mm256_set_ps(7, 6, 5, 4, 3, 2, 1, 0);
 }
 
 ICBC_FORCEINLINE VMask operator> (VFloat A, VFloat B) { return _mm256_cmp_ps(A, B, _CMP_GT_OQ); }
 ICBC_FORCEINLINE VMask operator>=(VFloat A, VFloat B) { return _mm256_cmp_ps(A, B, _CMP_GE_OQ); }
 ICBC_FORCEINLINE VMask operator< (VFloat A, VFloat B) { return _mm256_cmp_ps(A, B, _CMP_LT_OQ); }
 ICBC_FORCEINLINE VMask operator<=(VFloat A, VFloat B) { return _mm256_cmp_ps(A, B, _CMP_LE_OQ); }
 
 ICBC_FORCEINLINE VMask operator| (VMask A, VMask B) { return _mm256_or_ps(A, B); }
 ICBC_FORCEINLINE VMask operator& (VMask A, VMask B) { return _mm256_and_ps(A, B); }
 ICBC_FORCEINLINE VMask operator^ (VMask A, VMask B) { return _mm256_xor_ps(A, B); }
 
 // mask ? b : a
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) {
     return _mm256_blendv_ps(a, b, mask);
 }
 
 ICBC_FORCEINLINE VMask vbroadcast(bool b) { 
     return _mm256_castsi256_ps(_mm256_set1_epi32(-int32_t(b)));
 }
 
 ICBC_FORCEINLINE bool all(VMask m) {
     __m256 zero = _mm256_setzero_ps();
     return _mm256_testc_ps(_mm256_cmp_ps(zero, zero, _CMP_EQ_UQ), m) == 0;
 }
 
 ICBC_FORCEINLINE bool any(VMask m) {
     return _mm256_testz_ps(m, m) == 0;
 }
 
 ICBC_FORCEINLINE uint mask(VMask m) {
     return (uint)_mm256_movemask_ps(m);
 }
 
 // This is missing on some GCC versions.
 #if !defined _mm256_set_m128
 #define _mm256_set_m128(hi, lo) _mm256_insertf128_ps(_mm256_castps128_ps256(lo), (hi), 0x1)
 #endif
 
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) {
 
     __m128 vlow  = _mm256_castps256_ps128(v);
     __m128 vhigh = _mm256_extractf128_ps(v, 1);
            vlow  = _mm_min_ps(vlow, vhigh);
 
     // First do an horizontal reduction.                                // v = [ D C | B A ]
     __m128 shuf = _mm_shuffle_ps(vlow, vlow, _MM_SHUFFLE(2, 3, 0, 1));  //     [ C D | A B ]
     __m128 mins = _mm_min_ps(vlow, shuf);                            // mins = [ D+C C+D | B+A A+B ]
     shuf        = _mm_movehl_ps(shuf, mins);                         //        [   C   D | D+C C+D ]
     mins        = _mm_min_ss(mins, shuf);
 
     VFloat vmin = _mm256_permute_ps(_mm256_set_m128(mins, mins), 0); // _MM256_PERMUTE(0, 0, 0, 0, 0, 0, 0, 0)
 
     // Then find the index.
     uint mask = _mm256_movemask_ps(v <= vmin);
     return ctz(mask);
 }
 
 // AoS to SoA
 ICBC_FORCEINLINE void vtranspose4(VFloat & a, VFloat & b, VFloat & c, VFloat & d) {
     VFloat r0 = _mm256_unpacklo_ps(a, b);
     VFloat r1 = _mm256_unpacklo_ps(c, d);
     VFloat r2 = _mm256_permute2f128_ps(r0, r1, 0x20);
     VFloat r3 = _mm256_permute2f128_ps(r0, r1, 0x31);
     r0 = _mm256_unpackhi_ps(a, b);
     r1 = _mm256_unpackhi_ps(c, d);
     a = _mm256_unpacklo_ps(r2, r3);
     b = _mm256_unpackhi_ps(r2, r3);
     r2 = _mm256_permute2f128_ps(r0, r1, 0x20);
     r3 = _mm256_permute2f128_ps(r0, r1, 0x31);
     c = _mm256_unpacklo_ps(r2, r3);
     d = _mm256_unpackhi_ps(r2, r3);
 }
 
 ICBC_FORCEINLINE VInt vzeroi() {
     return _mm256_setzero_si256();
 }
 
 ICBC_FORCEINLINE VInt vbroadcast(int a) {
     return _mm256_set1_epi32(a);
 }
 
 ICBC_FORCEINLINE VInt vload(const int * ptr) {
     return _mm256_load_si256((const __m256i*)ptr);
 }
 
 ICBC_FORCEINLINE VInt operator- (VInt A, int b) { return _mm256_sub_epi32(A, _mm256_set1_epi32(b)); }
 ICBC_FORCEINLINE VInt operator& (VInt A, int b) { return _mm256_and_si256(A, _mm256_set1_epi32(b)); }
 ICBC_FORCEINLINE VInt operator>> (VInt A, int b) { return _mm256_srli_epi32(A, b); }
 
 ICBC_FORCEINLINE VMask operator> (VInt A, int b) { return _mm256_castsi256_ps(_mm256_cmpgt_epi32(A, _mm256_set1_epi32(b))); }
 ICBC_FORCEINLINE VMask operator>= (VInt A, int b) { return _mm256_castsi256_ps(_mm256_cmpgt_epi32(A, _mm256_set1_epi32(b-1))); }
 ICBC_FORCEINLINE VMask operator== (VInt A, int b) { return _mm256_castsi256_ps(_mm256_cmpeq_epi32(A, _mm256_set1_epi32(b))); }
 
 // mask ? v[idx] : 0
 ICBC_FORCEINLINE VFloat vpermuteif(VMask mask, VFloat v, VInt idx) {
     return _mm256_and_ps(_mm256_permutevar8x32_ps(v, idx), mask);
 }
 
 // mask ? (idx > 8 ? vhi[idx] : vlo[idx]) : 0
 ICBC_FORCEINLINE VFloat vpermute2if(VMask mask, VFloat vlo, VFloat vhi, VInt idx) {
 #if 0
     VMask mhi = idx > 7;
     vlo = _mm256_permutevar8x32_ps(vlo, idx);
     vhi = _mm256_permutevar8x32_ps(vhi, idx);
     VFloat v = _mm256_blendv_ps(vlo, vhi, mhi);
     return _mm256_and_ps(v, mask);
 #else
     // Fabian Giesen says not to mix _mm256_blendv_ps and _mm256_permutevar8x32_ps since they contend for the same gates and instead suggests the following:
     vhi = _mm256_xor_ps(vhi, vlo);
     VFloat v = _mm256_permutevar8x32_ps(vlo, idx);
     VMask mhi = idx > 7;
     v = _mm256_xor_ps(v, _mm256_and_ps(_mm256_permutevar8x32_ps(vhi, idx), mhi));
     return _mm256_and_ps(v, mask);
 #endif
 }
 
 
 #elif ICBC_SIMD == ICBC_AVX512
 
 constexpr int VEC_SIZE = 16;
 
 #if __GNUC__
 union VFloat {
     __m512 v;
     float m512_f32[VEC_SIZE];
 
     VFloat() {}
     VFloat(__m512 v) : v(v) {}
     operator __m512 & () { return v; }
 };
 union VInt {
     __m512i v;
     int m512i_i32[VEC_SIZE];
 
     VInt() {}
     VInt(__m512i v) : v(v) {}
     operator __m512i & () { return v; }
 };
 #else
 using VFloat = __m512;
 using VInt = __m512i;
 #endif
 struct VMask { __mmask16 m; };
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) {
     return v.m512_f32[i];
 }
 
 ICBC_FORCEINLINE VFloat vzero() {
     return _mm512_setzero_ps();
 }
 
 ICBC_FORCEINLINE VFloat vbroadcast(float a) {
     return _mm512_set1_ps(a);
 }
 
 ICBC_FORCEINLINE VFloat vload(const float * ptr) {
     return _mm512_load_ps(ptr);
 }
 
 ICBC_FORCEINLINE VFloat vload(VMask mask, const float * ptr) {
     return _mm512_mask_load_ps(_mm512_undefined(), mask.m, ptr);
 }
 
 ICBC_FORCEINLINE VFloat vload(VMask mask, const float * ptr, float fallback) {
     return _mm512_mask_load_ps(_mm512_set1_ps(fallback), mask.m, ptr);
 }
 
 ICBC_FORCEINLINE VFloat vgather(const float * base, VFloat index) {
     return _mm512_i32gather_ps(_mm512_cvtps_epi32(index), base, 4);
 }
 
 ICBC_FORCEINLINE VFloat operator+(VFloat a, VFloat b) {
     return _mm512_add_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator-(VFloat a, VFloat b) {
     return _mm512_sub_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator*(VFloat a, VFloat b) {
     return _mm512_mul_ps(a, b);
 }
 
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) {
 #if ICBC_USE_RCP
     VFloat r = _mm512_rcp14_ps(a);
 
     // r = r * (2 - r * a)
     return _mm512_mul_ps(r, _mm512_fnmadd_ps(r, a, vbroadcast(2.0f)));
 #else
     return _mm512_div_ps(vbroadcast(1.0f), a);
 #endif
 }
 
 // a*b+c
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) {
     return _mm512_fmadd_ps(a, b, c);
 }
 
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) {
     return _mm512_fmsub_ps(a, b, c);
 }
 
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) {
     return vmsub(a, b, c * d);
 }
 
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) {
     auto zero = _mm512_setzero_ps();
     auto one = _mm512_set1_ps(1.0f);
     return _mm512_min_ps(_mm512_max_ps(a, zero), one);
 }
 
 ICBC_FORCEINLINE VFloat vround01(VFloat a) {
     return _mm512_roundscale_ps(a, _MM_FROUND_TO_NEAREST_INT);
 }
 
 ICBC_FORCEINLINE VFloat vtruncate(VFloat a) {
     return _mm512_roundscale_ps(a, _MM_FROUND_TO_ZERO);
 }
 
 ICBC_FORCEINLINE VFloat lane_id() {
     return _mm512_set_ps(15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0);
 }
 
 ICBC_FORCEINLINE VMask operator> (VFloat A, VFloat B) { return { _mm512_cmp_ps_mask(A, B, _CMP_GT_OQ) }; }
 ICBC_FORCEINLINE VMask operator>=(VFloat A, VFloat B) { return { _mm512_cmp_ps_mask(A, B, _CMP_GE_OQ) }; }
 ICBC_FORCEINLINE VMask operator< (VFloat A, VFloat B) { return { _mm512_cmp_ps_mask(A, B, _CMP_LT_OQ) }; }
 ICBC_FORCEINLINE VMask operator<=(VFloat A, VFloat B) { return { _mm512_cmp_ps_mask(A, B, _CMP_LE_OQ) }; }
 
 ICBC_FORCEINLINE VMask operator! (VMask A) { return { _mm512_knot(A.m) }; }
 ICBC_FORCEINLINE VMask operator| (VMask A, VMask B) { return { _mm512_kor(A.m, B.m) }; }
 ICBC_FORCEINLINE VMask operator& (VMask A, VMask B) { return { _mm512_kand(A.m, B.m) }; }
 ICBC_FORCEINLINE VMask operator^ (VMask A, VMask B) { return { _mm512_kxor(A.m, B.m) }; }
 
 // mask ? b : a
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) {
     return _mm512_mask_blend_ps(mask.m, a, b);
 }
 
 ICBC_FORCEINLINE VMask vbroadcast(bool b) { 
     return { __mmask16(-int16_t(b)) };
 }
 
 ICBC_FORCEINLINE bool all(VMask mask) {
     return mask.m == 0xFFFFFFFF;
 }
 
 ICBC_FORCEINLINE bool any(VMask mask) {
     return mask.m != 0;
 }
 
 ICBC_FORCEINLINE uint mask(VMask mask) {
     return mask.m;
 }
 
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) {
 
     // First do an horizontal reduction.
     VFloat vmin = vbroadcast(_mm512_reduce_min_ps(v));
 
     // Then find the index.
     VMask mask = (v <= vmin);
     return ctz(mask.m);
 }
 
 //ICBC_FORCEINLINE void vtranspose4(VFloat & a, VFloat & b, VFloat & c, VFloat & d); // @@
 
 ICBC_FORCEINLINE int lane(VInt v, int i) {
     //return _mm256_extract_epi32(v, i);
     return v.m512i_i32[i];
 }
 
 ICBC_FORCEINLINE VInt vzeroi() {
     return _mm512_setzero_epi32();
 }
 
 ICBC_FORCEINLINE VInt vbroadcast(int a) {
     return _mm512_set1_epi32(a);
 }
 
 ICBC_FORCEINLINE VInt vload(const int * ptr) {
     return _mm512_load_epi32(ptr);
 }
 
 ICBC_FORCEINLINE VInt operator- (VInt A, int b) { return _mm512_sub_epi32(A, vbroadcast(b)); }
 ICBC_FORCEINLINE VInt operator& (VInt A, int b) { return _mm512_and_epi32(A, vbroadcast(b)); }
 ICBC_FORCEINLINE VInt operator>> (VInt A, int b) { return _mm512_srli_epi32(A, b); }
 
 ICBC_FORCEINLINE VMask operator> (VInt A, int b) { return { _mm512_cmpgt_epi32_mask(A, vbroadcast(b)) }; }
 ICBC_FORCEINLINE VMask operator>=(VInt A, int b) { return { _mm512_cmpge_epi32_mask(A, vbroadcast(b)) }; }
 ICBC_FORCEINLINE VMask operator== (VInt A, int b) { return { _mm512_cmpeq_epi32_mask(A, vbroadcast(b)) }; }
 
 // mask ? v[idx] : 0
 ICBC_FORCEINLINE VFloat vpermuteif(VMask mask, VFloat v, VInt idx) {
     return _mm512_maskz_permutexvar_ps(mask.m, idx, v);
 }
 
 
 #elif ICBC_SIMD == ICBC_NEON
 
 constexpr int VEC_SIZE = 4;
 
 #if __GNUC__
 union VFloat {
     float32x4_t v;
     float e[4];
     VFloat() {}
     VFloat(float32x4_t v) : v(v) {}
     operator float32x4_t & () { return v; }
 };
 struct VMask {
     uint32x4_t v;
     VMask() {}
     VMask(uint32x4_t v) : v(v) {}
     operator uint32x4_t & () { return v; }
 };
 #else
 using VFloat = float32x4_t;
 using VMask = uint32x4_t;
 #endif
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) {
 #if defined(__clang__)
     return v.e[i];
 #elif defined(_MSC_VER)
     return v.n128_f32[i];
 #else
     return v.v[i];
 #endif
 }
 
 ICBC_FORCEINLINE VFloat vzero() {
     return vdupq_n_f32(0.0f);
 }
 
 ICBC_FORCEINLINE VFloat vbroadcast(float a) {
     return vdupq_n_f32(a);
 }
 
 ICBC_FORCEINLINE VFloat vload(const float * ptr) {
     return vld1q_f32(ptr);
 }
 
 ICBC_FORCEINLINE VFloat operator+(VFloat a, VFloat b) {
     return vaddq_f32(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator-(VFloat a, VFloat b) {
     return vsubq_f32(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator*(VFloat a, VFloat b) {
     return vmulq_f32(a, b);
 }
 
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) {
 #if ICBC_USE_RCP
     return vrecpeq_f32(a);
 #else
     //VFloat rcp = vrecpeq_f32(a);
     //return vmulq_f32(vrecpsq_f32(a, rcp), rcp);
     return vdivq_f32(vbroadcast(1.0f), a);
 #endif
 }
 
 // a*b+c
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) {
     return vmlaq_f32(c, a, b);
 }
 
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) {
     return a * b - c;
 }
 
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) {
     // vmlsq_f32(a, b, c) == a - (b * c)
     return vmlsq_f32(a * b, c, d);
 }
 
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) {
     return vminq_f32(vmaxq_f32(a, vzero()), vbroadcast(1));
 }
 
 ICBC_FORCEINLINE VFloat vround01(VFloat a) {
     return vrndnq_f32(a);   // Round to integral (to nearest, ties to even)
     //return vcvtq_f32_s32(vcvtq_s32_f32(a + vbroadcast(0.5)));
 }
 
 ICBC_FORCEINLINE VFloat vtruncate(VFloat a) {
     return vrndq_f32(a);
     //return vcvtq_f32_s32(vcvtq_s32_f32(a));
 }
 
 ICBC_FORCEINLINE VFloat lane_id() {
     ICBC_ALIGN float data[4] = { 0, 1, 2, 3 };
 	return vld1q_f32(data);
 }
 
 ICBC_FORCEINLINE VMask operator> (VFloat A, VFloat B) { return { vcgtq_f32(A, B) }; }
 ICBC_FORCEINLINE VMask operator>=(VFloat A, VFloat B) { return { vcgeq_f32(A, B) }; }
 ICBC_FORCEINLINE VMask operator< (VFloat A, VFloat B) { return { vcltq_f32(A, B) }; }
 ICBC_FORCEINLINE VMask operator<=(VFloat A, VFloat B) { return { vcleq_f32(A, B) }; }
 
 ICBC_FORCEINLINE VMask operator| (VMask A, VMask B) { return { vorrq_u32(A, B) }; }
 ICBC_FORCEINLINE VMask operator& (VMask A, VMask B) { return { vandq_u32(A, B) }; }
 ICBC_FORCEINLINE VMask operator^ (VMask A, VMask B) { return { veorq_u32(A, B) }; }
 
 // mask ? b : a
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) {
     return vbslq_f32(mask, b, a);
 }
 
 ICBC_FORCEINLINE bool all(VMask mask) {
     return vminvq_u32(mask) != 0;
     // uint32x2_t m2 = vpmin_u32(vget_low_u32(mask), vget_high_u32(mask));
     // uint32x2_t m1 = vpmin_u32(m2, m2);
     // return vget_lane_u32(m1, 0) != 0;
 }
 
 ICBC_FORCEINLINE bool any(VMask mask) {
     return vmaxvq_u32(mask) != 0;
     // uint32x2_t m2 = vpmax_u32(vget_low_u32(mask), vget_high_u32(mask));
     // uint32x2_t m1 = vpmax_u32(m2, m2);
     // return vget_lane_u32(m1, 0) != 0;
 }
 
 ICBC_FORCEINLINE uint mask(VMask b) {
     const uint32x4_t movemask = { 1, 2, 4, 8 };
     return vaddvq_u32(vandq_u32(b, movemask));
 }
 
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) {
     // float32x2_t m2 = vpmin_f32(vget_low_u32(v), vget_high_u32(v));
     // float32x2_t m1 = vpmin_f32(m2, m2);
     // VFloat vmin = vdupq_lane_f32(m1, 0);
     VFloat vmin = vdupq_n_f32(vminvq_f32(v));
     return ctz(mask(v <= vmin));
 }
 
 // https://github.com/Maratyszcza/NNPACK/blob/master/src/neon/transpose.h
 ICBC_FORCEINLINE void vtranspose4(VFloat & a, VFloat & b, VFloat & c, VFloat & d) {
     // row0 = ( x00 x01 x02 x03 )
     // row1 = ( x10 x11 x12 x13 )
     // row2 = ( x20 x21 x22 x23 )
     // row3 = ( x30 x31 x32 x33 )
 
     // row01 = ( x00 x10 x02 x12 ), ( x01 x11 x03, x13 )
     // row23 = ( x20 x30 x22 x32 ), ( x21 x31 x23, x33 )
     float32x4x2_t row01 = vtrnq_f32(a, b);
     float32x4x2_t row23 = vtrnq_f32(c, d);
 
     // row0 = ( x00 x10 x20 x30 )
     // row1 = ( x01 x11 x21 x31 )
     // row2 = ( x02 x12 x22 x32 )
     // row3 = ( x03 x13 x23 x33 )
     a = vcombine_f32(vget_low_f32(row01.val[0]), vget_low_f32(row23.val[0]));
     b = vcombine_f32(vget_low_f32(row01.val[1]), vget_low_f32(row23.val[1]));
     c = vcombine_f32(vget_high_f32(row01.val[0]), vget_high_f32(row23.val[0]));
     d = vcombine_f32(vget_high_f32(row01.val[1]), vget_high_f32(row23.val[1]));
 }
 
 #elif ICBC_SIMD == ICBC_VMX
 
 constexpr int VEC_SIZE = 4;
 
 union VFloat {
     vectro float v;
     float e[4];
     VFloat() {}
     VFloat(vector float v) : v(v) {}
     operator vector float & () { return v; }
 };
 struct VMask {
     vector unsigned int v;
     VMask() {}
     VMask(vector unsigned int v) : v(v) {}
     operator vector unsigned int & () { return v; }
 };
 
 ICBC_FORCEINLINE float & lane(VFloat & v, int i) {
     return v.e[i];
 }
 
 ICBC_FORCEINLINE VFloat vzero() {
     return vec_splats(0.0f);
 }
 
 ICBC_FORCEINLINE VFloat vbroadcast(float a) {
     return vec_splats(a);
 }
 
 ICBC_FORCEINLINE VFloat vload(const float * ptr) {
     return vec_ld(ptr)
 }
 
 ICBC_FORCEINLINE VFloat operator+(VFloat a, VFloat b) {
     return vec_add(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator-(VFloat a, VFloat b) {
     return vec_sub(a, b);
 }
 
 ICBC_FORCEINLINE VFloat operator*(VFloat a, VFloat b) {
     return vec_madd(a, b, vec_splats(-0.0f));
 }
 
 ICBC_FORCEINLINE VFloat vrcp(VFloat a) {
 #if ICBC_USE_RCP
     // get the reciprocal estimate
     vector float estimate = vec_re( v.vec );
 
     // one round of Newton-Rhaphson refinement
     vector float diff = vec_nmsub( estimate, v.vec, vec_splats( 1.0f ) );
     return vec_madd(diff, estimate, estimate );
 #else
     return vec_div(vec_splats(1),a);
 #endif
 }
 
 // a*b+c
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, VFloat b, VFloat c) {
     return vec_madd(a, b, c);
 }
 
 ICBC_FORCEINLINE VFloat vmsub(VFloat a, VFloat b, VFloat c) {
     return vec_msub(a, b, c); // @@ Is this right?
 }
 
 ICBC_FORCEINLINE VFloat vm2sub(VFloat a, VFloat b, VFloat c, VFloat d) {
     return vmsub(a, b, c * d);
 }
 
 ICBC_FORCEINLINE VFloat vsaturate(VFloat a) {
     return vec_min(vec_max(a, vzero()), vbroadcast(1));
 }
 
 ICBC_FORCEINLINE VFloat vround01(VFloat a) {
     // @@ Assumes a is positive and ~small
     return vec_trunc(a + vbroadcast(0.5));
 }
 
 ICBC_FORCEINLINE VFloat vtruncate(VFloat a) {
     return vec_trunc(a);
 }
 
 ICBC_FORCEINLINE VFloat lane_id() {
     return (VFloat){ 0, 1, 2, 3 };
 }
 
 ICBC_FORCEINLINE VMask operator> (VFloat A, VFloat B) { return { vec_cmpgt(A, B) }; }
 ICBC_FORCEINLINE VMask operator>=(VFloat A, VFloat B) { return { vec_cmpge(A, B) }; }
 ICBC_FORCEINLINE VMask operator< (VFloat A, VFloat B) { return { vec_cmplt(A, B) }; }
 ICBC_FORCEINLINE VMask operator<=(VFloat A, VFloat B) { return { vec_cmple(A, B) }; }
 
 ICBC_FORCEINLINE VMask operator| (VMask A, VMask B) { return { vec_or(A, B) }; }
 ICBC_FORCEINLINE VMask operator& (VMask A, VMask B) { return { vec_and(A, B) }; }
 ICBC_FORCEINLINE VMask operator^ (VMask A, VMask B) { return { vec_xor(A, B) }; }
 
 // mask ? b : a
 ICBC_FORCEINLINE VFloat vselect(VMask mask, VFloat a, VFloat b) {
     return vec_sel(a, b, mask);
 }
 
 ICBC_FORCEINLINE int reduce_min_index(VFloat v) {
 
     //VFloat vmin =  //@@ Horizontal min?
     //return vec_cmpeq_idx(v, vmin);
 
     // @@ Is there a better way to do this reduction?
     int min_idx = 0;
     float min_value = lane(v, 0);
 
     for (int i = 1; i < VEC_SIZE; i++) {
         float value = lane(v, i);
         if (value < min_value) {
             min_value = value;
             min_idx = i;
         }
     }
 
     return min_idx;
 }
 
 ICBC_FORCEINLINE void vtranspose4(VFloat & a, VFloat & b, VFloat & c, VFloat & d) {
     VFloat t1 = vec_mergeh(a, c);
     VFloat t2 = vec_mergel(a, c);
     VFloat t3 = vec_mergeh(b, d);
     VFloat t4 = vec_mergel(b, d);
     a = vec_mergeh(t1, t3);
     b = vec_mergel(t1, t3);
     c = vec_mergeh(t2, t4);
     d = vec_mergel(t2, t4);
 }
 
 #endif // ICBC_SIMD == *
 
 #if ICBC_SIMD != ICBC_SCALAR
 ICBC_FORCEINLINE VFloat vmadd(VFloat a, float b, VFloat c) {
     VFloat vb = vbroadcast(b);
     return vmadd(a, vb, c);
 }
 #endif
 
 struct VVector3 {
     VFloat x;
     VFloat y;
     VFloat z;
 };
 
 ICBC_FORCEINLINE VVector3 vbroadcast(Vector3 v) {
     VVector3 v8;
 
     v8.x = vbroadcast(v.x);
     v8.y = vbroadcast(v.y);
     v8.z = vbroadcast(v.z);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 vbroadcast(float x, float y, float z) {
     VVector3 v8;
     v8.x = vbroadcast(x);
     v8.y = vbroadcast(y);
     v8.z = vbroadcast(z);
     return v8;
 }
 
 /*ICBC_FORCEINLINE VVector3 vload(const Vector3 * v) {
     // @@ See this for a 8x3 transpose: https://software.intel.com/content/www/us/en/develop/articles/3d-vector-normalization-using-256-bit-intel-advanced-vector-extensions-intel-avx.html
 }*/
 
 ICBC_FORCEINLINE VVector3 vload(const Vector4 * ptr) {
 #if ICBC_SIMD == ICBC_AVX512
 
     // @@ AVX512 transpose not implemented.
     __m512i vindex = _mm512_set_epi32(4 * 15, 4 * 14, 4 * 13, 4 * 12, 4 * 11, 4 * 10, 4 * 9, 4 * 8, 4 * 7, 4 * 6, 4 * 5, 4 * 4, 4 * 3, 4 * 2, 4 * 1, 0);
 
     VVector3 v;
     v.x = _mm512_i32gather_ps(vindex, &ptr->x, 4);
     v.y = _mm512_i32gather_ps(vindex, &ptr->y, 4);
     v.z = _mm512_i32gather_ps(vindex, &ptr->z, 4);
     return v;
 
 #else
 
     VVector3 v;
     v.x = vload(&ptr->x + 0 * VEC_SIZE);
     v.y = vload(&ptr->x + 1 * VEC_SIZE);
     v.z = vload(&ptr->x + 2 * VEC_SIZE);
     VFloat tmp = vload(&ptr->x + 3 * VEC_SIZE);
 
     vtranspose4(v.x, v.y, v.z, tmp);
 
     return v;
 #endif
 }
 
 
 ICBC_FORCEINLINE VVector3 operator+(VVector3 a, VVector3 b) {
     VVector3 v;
     v.x = (a.x + b.x);
     v.y = (a.y + b.y);
     v.z = (a.z + b.z);
     return v;
 }
 
 ICBC_FORCEINLINE VVector3 operator-(VVector3 a, VVector3 b) {
     VVector3 v8;
     v8.x = (a.x - b.x);
     v8.y = (a.y - b.y);
     v8.z = (a.z - b.z);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 operator*(VVector3 a, VVector3 b) {
     VVector3 v8;
     v8.x = (a.x * b.x);
     v8.y = (a.y * b.y);
     v8.z = (a.z * b.z);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 operator*(VVector3 a, VFloat b) {
     VVector3 v8;
     v8.x = (a.x * b);
     v8.y = (a.y * b);
     v8.z = (a.z * b);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 vmadd(VVector3 a, VVector3 b, VVector3 c) {
     VVector3 v8;
     v8.x = vmadd(a.x, b.x, c.x);
     v8.y = vmadd(a.y, b.y, c.y);
     v8.z = vmadd(a.z, b.z, c.z);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 vmadd(VVector3 a, VFloat b, VVector3 c) {
     VVector3 v8;
     v8.x = vmadd(a.x, b, c.x);
     v8.y = vmadd(a.y, b, c.y);
     v8.z = vmadd(a.z, b, c.z);
     return v8;
 }
 
 #if ICBC_SIMD != ICBC_SCALAR
 ICBC_FORCEINLINE VVector3 vmadd(VVector3 a, float b, VVector3 c) {
     VVector3 v8;
     VFloat vb = vbroadcast(b);
     v8.x = vmadd(a.x, vb, c.x);
     v8.y = vmadd(a.y, vb, c.y);
     v8.z = vmadd(a.z, vb, c.z);
     return v8;
 }
 #endif
 
 ICBC_FORCEINLINE VVector3 vmsub(VVector3 a, VFloat b, VFloat c) {
     VVector3 v8;
     v8.x = vmsub(a.x, b, c);
     v8.y = vmsub(a.y, b, c);
     v8.z = vmsub(a.z, b, c);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 vmsub(VVector3 a, VFloat b, VVector3 c) {
     VVector3 v8;
     v8.x = vmsub(a.x, b, c.x);
     v8.y = vmsub(a.y, b, c.y);
     v8.z = vmsub(a.z, b, c.z);
     return v8;
 }
 
 ICBC_FORCEINLINE VVector3 vm2sub(VVector3 a, VFloat b, VVector3 c, VFloat d) {
     VVector3 v;
     v.x = vm2sub(a.x, b, c.x, d);
     v.y = vm2sub(a.y, b, c.y, d);
     v.z = vm2sub(a.z, b, c.z, d);
     return v;
 }
 
 ICBC_FORCEINLINE VVector3 vm2sub(VVector3 a, VVector3 b, VVector3 c, VFloat d) {
     VVector3 v;
     v.x = vm2sub(a.x, b.x, c.x, d);
     v.y = vm2sub(a.y, b.y, c.y, d);
     v.z = vm2sub(a.z, b.z, c.z, d);
     return v;
 }
 
 ICBC_FORCEINLINE VVector3 vm2sub(VVector3 a, VVector3 b, VVector3 c, VVector3 d) {
     VVector3 v;
     v.x = vm2sub(a.x, b.x, c.x, d.x);
     v.y = vm2sub(a.y, b.y, c.y, d.y);
     v.z = vm2sub(a.z, b.z, c.z, d.z);
     return v;
 }
 
 ICBC_FORCEINLINE VFloat vdot(VVector3 a, VVector3 b) {
     VFloat r;
     r = a.x * b.x + a.y * b.y + a.z * b.z;
     return r;
 }
 
 // Length squared.
 ICBC_FORCEINLINE VFloat vlen2(VVector3 v) {
     return vdot(v, v);
 }
 
 
 // mask ? b : a
 ICBC_FORCEINLINE VVector3 vselect(VMask mask, VVector3 a, VVector3 b) {
     VVector3 r;
     r.x = vselect(mask, a.x, b.x);
     r.y = vselect(mask, a.y, b.y);
     r.z = vselect(mask, a.z, b.z);
     return r;
 }
 
 
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Color conversion functions.
 
 static const float midpoints5[32] = {
     0.015686f, 0.047059f, 0.078431f, 0.111765f, 0.145098f, 0.176471f, 0.207843f, 0.241176f, 0.274510f, 0.305882f, 0.337255f, 0.370588f, 0.403922f, 0.435294f, 0.466667f, 0.5f,
     0.533333f, 0.564706f, 0.596078f, 0.629412f, 0.662745f, 0.694118f, 0.725490f, 0.758824f, 0.792157f, 0.823529f, 0.854902f, 0.888235f, 0.921569f, 0.952941f, 0.984314f, FLT_MAX
 };
 
 static const float midpoints6[64] = {
     0.007843f, 0.023529f, 0.039216f, 0.054902f, 0.070588f, 0.086275f, 0.101961f, 0.117647f, 0.133333f, 0.149020f, 0.164706f, 0.180392f, 0.196078f, 0.211765f, 0.227451f, 0.245098f,
     0.262745f, 0.278431f, 0.294118f, 0.309804f, 0.325490f, 0.341176f, 0.356863f, 0.372549f, 0.388235f, 0.403922f, 0.419608f, 0.435294f, 0.450980f, 0.466667f, 0.482353f, 0.500000f,
     0.517647f, 0.533333f, 0.549020f, 0.564706f, 0.580392f, 0.596078f, 0.611765f, 0.627451f, 0.643137f, 0.658824f, 0.674510f, 0.690196f, 0.705882f, 0.721569f, 0.737255f, 0.754902f,
     0.772549f, 0.788235f, 0.803922f, 0.819608f, 0.835294f, 0.850980f, 0.866667f, 0.882353f, 0.898039f, 0.913725f, 0.929412f, 0.945098f, 0.960784f, 0.976471f, 0.992157f, FLT_MAX
 };
 
 /*void init_tables() {
     for (int i = 0; i < 31; i++) {
         float f0 = float(((i+0) << 3) | ((i+0) >> 2)) / 255.0f;
         float f1 = float(((i+1) << 3) | ((i+1) >> 2)) / 255.0f;
         midpoints5[i] = (f0 + f1) * 0.5;
     }
     midpoints5[31] = FLT_MAX;
 
     for (int i = 0; i < 63; i++) {
         float f0 = float(((i+0) << 2) | ((i+0) >> 4)) / 255.0f;
         float f1 = float(((i+1) << 2) | ((i+1) >> 4)) / 255.0f;
         midpoints6[i] = (f0 + f1) * 0.5;
     }
     midpoints6[63] = FLT_MAX;
 }*/
 
 ICBC_FORCEINLINE VFloat vround5(VFloat x) {
     const VFloat rb_scale = vbroadcast(31.0f);
     const VFloat rb_inv_scale = vbroadcast(1.0f / 31.0f);
 
 #if ICBC_PERFECT_ROUND
     VFloat q = vtruncate(x * rb_scale);
     VFloat mp = vgather(midpoints5, q);
     //return (q + (vbroadcast(1.0f) & (x > mp))) * rb_inv_scale;
     return (q + vselect(x > mp, vzero(), vbroadcast(1))) * rb_inv_scale;
 #else
     return vround01(x * rb_scale) * rb_inv_scale;
 #endif
 }
 
 ICBC_FORCEINLINE VFloat vround6(VFloat x) {
     const VFloat g_scale = vbroadcast(63.0f);
     const VFloat g_inv_scale = vbroadcast(1.0f / 63.0f);
 
 #if ICBC_PERFECT_ROUND
     VFloat q = vtruncate(x * g_scale);
     VFloat mp = vgather(midpoints6, q);
     //return (q + (vbroadcast(1) & (x > mp))) * g_inv_scale;
     return (q + vselect(x > mp, vzero(), vbroadcast(1))) * g_inv_scale;
 #else
     return vround01(x * g_scale) * g_inv_scale;
 #endif
 }
 
 ICBC_FORCEINLINE VVector3 vround_ept(VVector3 v) {
     VVector3 r;
     r.x = vround5(vsaturate(v.x));
     r.y = vround6(vsaturate(v.y));
     r.z = vround5(vsaturate(v.z));
     return r;
 }
 
 static Color16 vector3_to_color16(const Vector3 & v) {
 
     // Truncate.
     uint r = uint(clamp(v.x * 31.0f, 0.0f, 31.0f));
 	uint g = uint(clamp(v.y * 63.0f, 0.0f, 63.0f));
 	uint b = uint(clamp(v.z * 31.0f, 0.0f, 31.0f));
 
     // Round exactly according to 565 bit-expansion.
     r += (v.x > midpoints5[r]);
     g += (v.y > midpoints6[g]);
     b += (v.z > midpoints5[b]);
 
     Color16 c;
     c.u = uint16((r << 11) | (g << 5) | b);
     return c;
 }
 
 static Color32 bitexpand_color16_to_color32(Color16 c16) {
     Color32 c32;
     //c32.b = (c16.b << 3) | (c16.b >> 2);
     //c32.g = (c16.g << 2) | (c16.g >> 4);
     //c32.r = (c16.r << 3) | (c16.r >> 2);
     //c32.a = 0xFF;
 
     c32.u = ((c16.u << 3) & 0xf8) | ((c16.u << 5) & 0xfc00) | ((c16.u << 8) & 0xf80000);
     c32.u |= (c32.u >> 5) & 0x070007;
     c32.u |= (c32.u >> 6) & 0x000300;
     c32.u |= 0xff000000;	
 
     return c32;
 }
 
 inline Vector3 color_to_vector3(Color32 c) {
     return { c.r / 255.0f, c.g / 255.0f, c.b / 255.0f };
 }
 
 inline Color32 vector3_to_color32(Vector3 v) {
     Color32 color;
     color.r = uint8(saturate(v.x) * 255 + 0.5f);
     color.g = uint8(saturate(v.y) * 255 + 0.5f);
     color.b = uint8(saturate(v.z) * 255 + 0.5f);
     color.a = 255;
     return color;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Input block processing.
 
 inline bool is_black(Vector3 c) {
     // This large threshold seems to improve compression. This is not forcing these texels to be black, just
     // causes them to be ignored during PCA.
     // @@ We may want to adjust this based on the quality level, since this increases the number of blocks that try cluster-3 fitting.
     //return c.x < midpoints5[0] && c.y < midpoints6[0] && c.z < midpoints5[0];
     //return c.x < 1.0f / 32 && c.y < 1.0f / 32 && c.z < 1.0f / 32;
     return c.x < 1.0f / 8 && c.y < 1.0f / 8 && c.z < 1.0f / 8;
 }
 
 // Find similar colors and combine them together.
 static int reduce_colors(const Vector4 * input_colors, const float * input_weights, int count, float threshold, Vector3 * colors, float * weights, bool * any_black)
 {
 #if 0
     for (int i = 0; i < 16; i++) {
         colors[i] = input_colors[i].xyz;
         weights[i] = input_weights[i];
     }
     return 16;
 #else
     *any_black = false;
 
     int n = 0;
     for (int i = 0; i < count; i++)
     {
         Vector3 ci = input_colors[i].xyz;
         float wi = input_weights[i];
 
         if (wi > 0) {
             // Find matching color.
             int j;
             for (j = 0; j < n; j++) {
                 if (equal(colors[j], ci, threshold)) {
                     colors[j] = (colors[j] * weights[j] + ci * wi) / (weights[j] + wi);
                     weights[j] += wi;
                     break;
                 }
             }
 
             // No match found. Add new color.
             if (j == n) {
                 colors[n] = ci;
                 weights[n] = wi;
                 n++;
             }
 
             if (is_black(ci)) {
                 *any_black = true;
             }
         }
     }
 
     ICBC_ASSERT(n <= count);
 
     return n;
 #endif
 }
 
 static int skip_blacks(const Vector3 * input_colors, const float * input_weights, int count, Vector3 * colors, float * weights)
 {
     int n = 0;
     for (int i = 0; i < count; i++)
     {
         Vector3 ci = input_colors[i];
         float wi = input_weights[i];
 
         if (is_black(ci)) {
             continue;
         }
 
         colors[n] = ci;
         weights[n] = wi;
         n += 1;
     }
 
     return n;
 }
 
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // PCA
 
 static Vector3 computeCentroid(int n, const Vector3 *__restrict points, const float *__restrict weights)
 {
     Vector3 centroid = { 0 };
     float total = 0.0f;
 
     for (int i = 0; i < n; i++)
     {
         total += weights[i];
         centroid += weights[i] * points[i];
     }
     centroid *= (1.0f / total);
 
     return centroid;
 }
 
 static Vector3 computeCovariance(int n, const Vector3 *__restrict points, const float *__restrict weights, float *__restrict covariance)
 {
     // compute the centroid
     Vector3 centroid = computeCentroid(n, points, weights);
 
     // compute covariance matrix
     for (int i = 0; i < 6; i++)
     {
         covariance[i] = 0.0f;
     }
 
     for (int i = 0; i < n; i++)
     {
         Vector3 a = (points[i] - centroid);    // @@ I think weight should be squared, but that seems to increase the error slightly.
         Vector3 b = weights[i] * a;
 
         covariance[0] += a.x * b.x;
         covariance[1] += a.x * b.y;
         covariance[2] += a.x * b.z;
         covariance[3] += a.y * b.y;
         covariance[4] += a.y * b.z;
         covariance[5] += a.z * b.z;
     }
 
     return centroid;
 }
 
 // @@ We should be able to do something cheaper...
 static Vector3 estimatePrincipalComponent(const float * __restrict matrix)
 {
     const Vector3 row0 = { matrix[0], matrix[1], matrix[2] };
     const Vector3 row1 = { matrix[1], matrix[3], matrix[4] };
     const Vector3 row2 = { matrix[2], matrix[4], matrix[5] };
 
     float r0 = lengthSquared(row0);
     float r1 = lengthSquared(row1);
     float r2 = lengthSquared(row2);
 
     if (r0 > r1 && r0 > r2) return row0;
     if (r1 > r2) return row1;
     return row2;
 }
 
 static inline Vector3 firstEigenVector_PowerMethod(const float *__restrict matrix)
 {
     if (matrix[0] == 0 && matrix[3] == 0 && matrix[5] == 0)
     {
         return {0};
     }
 
     Vector3 v = estimatePrincipalComponent(matrix);
 
     const int NUM = 6;
     for (int i = 0; i < NUM; i++)
     {
         float x = v.x * matrix[0] + v.y * matrix[1] + v.z * matrix[2];
         float y = v.x * matrix[1] + v.y * matrix[3] + v.z * matrix[4];
         float z = v.x * matrix[2] + v.y * matrix[4] + v.z * matrix[5];
 
         float norm = max(max(x, y), z);
 
         v = { x, y, z };
         v *= (1.0f / norm);
     }
 
     return v;
 }
 
 static Vector3 computePrincipalComponent_PowerMethod(int n, const Vector3 *__restrict points, const float *__restrict weights)
 {
     float matrix[6];
     computeCovariance(n, points, weights, matrix);
 
     return firstEigenVector_PowerMethod(matrix);
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // SAT
 
 struct SummedAreaTable {
     ICBC_ALIGN float r[16];
     ICBC_ALIGN float g[16];
     ICBC_ALIGN float b[16];
     ICBC_ALIGN float w[16];
 };
 
 int compute_sat(const Vector3 * colors, const float * weights, int count, SummedAreaTable * sat)
 {
     // I've tried using a lower quality approximation of the principal direction, but the best fit line seems to produce best results.
     Vector3 principal = computePrincipalComponent_PowerMethod(count, colors, weights);
 
     // build the list of values
     int order[16];
     float dps[16];
     for (int i = 0; i < count; ++i)
     {
         order[i] = i;
         dps[i] = dot(colors[i], principal);
     }
 
     // stable sort
     for (int i = 0; i < count; ++i)
     {
         for (int j = i; j > 0 && dps[j] < dps[j - 1]; --j)
         {
             swap(dps[j], dps[j - 1]);
             swap(order[j], order[j - 1]);
         }
     }
 
 #if 0
     u8 byte_order[16*4];
     for (int i = 0; i < count; ++i)
     {
     }
     
     // @@ Load colors and weights.
     VFloat table_r[16/VEC_SIZE];
     VFloat table_g[16/VEC_SIZE];
     VFloat table_b[16/VEC_SIZE];
     VFloat table_w[16/VEC_SIZE];
 
     VFloat sat_r[16/VEC_SIZE];
     VFloat sat_g[16/VEC_SIZE];
     VFloat sat_b[16/VEC_SIZE];
     VFloat sat_w[16/VEC_SIZE];
 
     for (int i = 0; i < 16; i += VEC_SIZE) {
         VFloat w = table_lookup(table_w, order);
         sat_r[i] = table_lookup(table_r, order) * w;
         sat_g[i] = table_lookup(table_g, order) * w;
         sat_b[i] = table_lookup(table_b, order) * w;
         sat_w[i] = w;
     }
 
     // Somewhat brute force way to compute sum tables. Is it worth it?
     for (int i = 0; i < 16; i += 1) {
         float r = sat_r[i/VEC_SIZE][i%VEC_SIZE];
 
         if (i > lane_id)
             sat_r[i] += sat_r[lane_id];
     }
 #endif
 
     float w = weights[order[0]];
     sat->r[0] = colors[order[0]].x * w;
     sat->g[0] = colors[order[0]].y * w;
     sat->b[0] = colors[order[0]].z * w;
     sat->w[0] = w;
 
     for (int i = 1; i < count; i++) {
         w = weights[order[i]];
         sat->r[i] = sat->r[i - 1] + colors[order[i]].x * w;
         sat->g[i] = sat->g[i - 1] + colors[order[i]].y * w;
         sat->b[i] = sat->b[i - 1] + colors[order[i]].z * w;
         sat->w[i] = sat->w[i - 1] + w;
     }
 
     for (int i = count; i < 16; i++) {
         sat->r[i] = FLT_MAX;
         sat->g[i] = FLT_MAX;
         sat->b[i] = FLT_MAX;
         sat->w[i] = FLT_MAX;
     }
 
     return count;
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Cluster Fit
 
 struct Combinations {
     uint8 c0, c1, c2, pad;
 };
 
 static ICBC_ALIGN int s_fourClusterTotal[16];
 static ICBC_ALIGN int s_threeClusterTotal[16];
 static ICBC_ALIGN Combinations s_fourCluster[968 + 8];
 static ICBC_ALIGN Combinations s_threeCluster[152 + 8];
 
 #if ICBC_USE_NEON_VTL || ICBC_USE_SSSE3_SHUFFLEB
 static uint8 s_byte_vtl_index0_4[4 * 968];
 static uint8 s_byte_vtl_index1_4[4 * 968];
 static uint8 s_byte_vtl_index2_4[4 * 968];
 
 static uint8 s_byte_vtl_index0_3[4 * 152];
 static uint8 s_byte_vtl_index1_3[4 * 152];
 #endif
 
 static float w1_3, w2_3;
 static float w1_9, w4_9, w2_9;
 
 
 static void init_cluster_tables() {
 
     for (int t = 1, i = 0; t <= 16; t++) {
         for (int c0 = 0; c0 <= t; c0++) {
             for (int c1 = 0; c1 <= t - c0; c1++) {
                 for (int c2 = 0; c2 <= t - c0 - c1; c2++) {
 
                     // Skip this cluster so that the total is a multiple of 8
                     if (c0 == 0 && c1 == 0 && c2 == 0) continue;
 
                     bool found = false;
                     if (t > 1) {
                         for (int j = 0; j < s_fourClusterTotal[t-2]; j++) {
                             if (s_fourCluster[j].c0 == c0 && s_fourCluster[j].c1 == c0+c1 && s_fourCluster[j].c2 == c0+c1+c2) {
                                 found = true;
                                 break;
                             }
                         }
                     }
 
                     if (!found) {
                         s_fourCluster[i].c0 = uint8(c0);
                         s_fourCluster[i].c1 = uint8(c0+c1);
                         s_fourCluster[i].c2 = uint8(c0+c1+c2);
                         i++;
                     }
                 }
             }
         }
 
         s_fourClusterTotal[t - 1] = i;
     }
 
     // Replicate last entry.
     for (int i = 0; i < 8; i++) {
         s_fourCluster[968 + i] = s_fourCluster[968-1];
     }
 
     for (int t = 1, i = 0; t <= 16; t++) {
         for (int c0 = 0; c0 <= t; c0++) {
             for (int c1 = 0; c1 <= t - c0; c1++) {
 
                 // Skip this cluster so that the total is a multiple of 8
                 if (c0 == 0 && c1 == 0) continue;
 
                 bool found = false;
                 if (t > 1) {
                     for (int j = 0; j < s_threeClusterTotal[t - 2]; j++) {
                         if (s_threeCluster[j].c0 == c0 && s_threeCluster[j].c1 == c0+c1) {
                             found = true;
                             break;
                         }
                     }
                 }
 
                 if (!found) {
                     s_threeCluster[i].c0 = uint8(c0);
                     s_threeCluster[i].c1 = uint8(c0 + c1);
                     i++;
                 }
             }
         }
 
         s_threeClusterTotal[t - 1] = i;
     }
 
     // Replicate last entry.
     for (int i = 0; i < 8; i++) {
         s_threeCluster[152 + i] = s_threeCluster[152 - 1];
     }
 
 #if ICBC_USE_NEON_VTL || ICBC_USE_SSSE3_SHUFFLEB
     for (int i = 0; i < 968; i++) {
         int c0 = (s_fourCluster[i].c0) - 1;
         s_byte_vtl_index0_4[4 * i + 0] = c0 >= 0 ? uint8(c0 * 4 + 0) : 255;
         s_byte_vtl_index0_4[4 * i + 1] = c0 >= 0 ? uint8(c0 * 4 + 1) : 255;
         s_byte_vtl_index0_4[4 * i + 2] = c0 >= 0 ? uint8(c0 * 4 + 2) : 255;
         s_byte_vtl_index0_4[4 * i + 3] = c0 >= 0 ? uint8(c0 * 4 + 3) : 255;
 
         int c1 = (s_fourCluster[i].c1) - 1;
         s_byte_vtl_index1_4[4 * i + 0] = c1 >= 0 ? uint8(c1 * 4 + 0) : 255;
         s_byte_vtl_index1_4[4 * i + 1] = c1 >= 0 ? uint8(c1 * 4 + 1) : 255;
         s_byte_vtl_index1_4[4 * i + 2] = c1 >= 0 ? uint8(c1 * 4 + 2) : 255;
         s_byte_vtl_index1_4[4 * i + 3] = c1 >= 0 ? uint8(c1 * 4 + 3) : 255;
 
         int c2 = (s_fourCluster[i].c2) - 1;
         s_byte_vtl_index2_4[4 * i + 0] = c2 >= 0 ? uint8(c2 * 4 + 0) : 255;
         s_byte_vtl_index2_4[4 * i + 1] = c2 >= 0 ? uint8(c2 * 4 + 1) : 255;
         s_byte_vtl_index2_4[4 * i + 2] = c2 >= 0 ? uint8(c2 * 4 + 2) : 255;
         s_byte_vtl_index2_4[4 * i + 3] = c2 >= 0 ? uint8(c2 * 4 + 3) : 255;
     }
 
     for (int i = 0; i < 152; i++) {
         int c0 = (s_threeCluster[i].c0) - 1;
         s_byte_vtl_index0_3[4 * i + 0] = c0 >= 0 ? uint8(c0 * 4 + 0) : 255;
         s_byte_vtl_index0_3[4 * i + 1] = c0 >= 0 ? uint8(c0 * 4 + 1) : 255;
         s_byte_vtl_index0_3[4 * i + 2] = c0 >= 0 ? uint8(c0 * 4 + 2) : 255;
         s_byte_vtl_index0_3[4 * i + 3] = c0 >= 0 ? uint8(c0 * 4 + 3) : 255;
 
         int c1 = (s_threeCluster[i].c1) - 1;
         s_byte_vtl_index1_3[4 * i + 0] = c1 >= 0 ? uint8(c1 * 4 + 0) : 255;
         s_byte_vtl_index1_3[4 * i + 1] = c1 >= 0 ? uint8(c1 * 4 + 1) : 255;
         s_byte_vtl_index1_3[4 * i + 2] = c1 >= 0 ? uint8(c1 * 4 + 2) : 255;
         s_byte_vtl_index1_3[4 * i + 3] = c1 >= 0 ? uint8(c1 * 4 + 3) : 255;
     }
 #endif
 }
 
 
 
 static void cluster_fit_three(const SummedAreaTable & sat, int count, Vector3 metric_sqr, Vector3 * start, Vector3 * end)
 {
     const float r_sum = sat.r[count-1];
     const float g_sum = sat.g[count-1];
     const float b_sum = sat.b[count-1];
     const float w_sum = sat.w[count-1];
 
     VFloat vbesterror = vbroadcast(FLT_MAX);
     VVector3 vbeststart = { vzero(), vzero(), vzero() };
     VVector3 vbestend = { vzero(), vzero(), vzero() };
 
     // check all possible clusters for this total order
     const int total_order_count = s_threeClusterTotal[count - 1];
 
     for (int i = 0; i < total_order_count; i += VEC_SIZE)
     {
         VVector3 x0, x1;
         VFloat w0, w1;
 
 #if ICBC_USE_AVX512_PERMUTE
 
         auto loadmask = lane_id() < vbroadcast(float(count));
 
         // Load sat in one register:
         VFloat vrsat = vload(loadmask, sat.r, FLT_MAX);
         VFloat vgsat = vload(loadmask, sat.g, FLT_MAX);
         VFloat vbsat = vload(loadmask, sat.b, FLT_MAX);
         VFloat vwsat = vload(loadmask, sat.w, FLT_MAX);
 
         // Load 4 uint8 per lane.
         VInt packedClusterIndex = vload((int *)&s_threeCluster[i]);
 
         VInt c0 = (packedClusterIndex & 0xFF) - 1;
         VInt c1 = (packedClusterIndex >> 8) - 1;
 
         x0.x = vpermuteif(c0 >= 0, vrsat, c0);
         x0.y = vpermuteif(c0 >= 0, vgsat, c0);
         x0.z = vpermuteif(c0 >= 0, vbsat, c0);
         w0   = vpermuteif(c0 >= 0, vwsat, c0);
 
         x1.x = vpermuteif(c1 >= 0, vrsat, c1);
         x1.y = vpermuteif(c1 >= 0, vgsat, c1);
         x1.z = vpermuteif(c1 >= 0, vbsat, c1);
         w1   = vpermuteif(c1 >= 0, vwsat, c1);
 
 #elif ICBC_USE_AVX2_PERMUTE2
 
         // Load 4 uint8 per lane. @@ Ideally I should pack this better and load only 2.
         VInt packedClusterIndex = vload((int *)&s_threeCluster[i]);
 
         VInt c0 = (packedClusterIndex & 0xFF) - 1;
         VInt c1 = (packedClusterIndex >> 8) - 1; // No need for & 0xFF
 
         if (count <= 8) {
             // Load sat.r in one register:
             VFloat rLo = vload(sat.r);
             VFloat gLo = vload(sat.g);
             VFloat bLo = vload(sat.b);
             VFloat wLo = vload(sat.w);
 
             x0.x = vpermuteif(c0 >= 0, rLo, c0);
             x0.y = vpermuteif(c0 >= 0, gLo, c0);
             x0.z = vpermuteif(c0 >= 0, bLo, c0);
             w0   = vpermuteif(c0 >= 0, wLo, c0);
 
             x1.x = vpermuteif(c1 >= 0, rLo, c1);
             x1.y = vpermuteif(c1 >= 0, gLo, c1);
             x1.z = vpermuteif(c1 >= 0, bLo, c1);
             w1   = vpermuteif(c1 >= 0, wLo, c1);
         }
         else {
             // Load sat.r in two registers:
             VFloat rLo = vload(sat.r); VFloat rHi = vload(sat.r + 8);
             VFloat gLo = vload(sat.g); VFloat gHi = vload(sat.g + 8);
             VFloat bLo = vload(sat.b); VFloat bHi = vload(sat.b + 8);
             VFloat wLo = vload(sat.w); VFloat wHi = vload(sat.w + 8);
 
             x0.x = vpermute2if(c0 >= 0, rLo, rHi, c0);
             x0.y = vpermute2if(c0 >= 0, gLo, gHi, c0);
             x0.z = vpermute2if(c0 >= 0, bLo, bHi, c0);
             w0   = vpermute2if(c0 >= 0, wLo, wHi, c0);
 
             x1.x = vpermute2if(c1 >= 0, rLo, rHi, c1);
             x1.y = vpermute2if(c1 >= 0, gLo, gHi, c1);
             x1.z = vpermute2if(c1 >= 0, bLo, bHi, c1);
             w1   = vpermute2if(c1 >= 0, wLo, wHi, c1);
         }
 
 #elif ICBC_USE_SSSE3_SHUFFLEB
 
         if (count <= 8) {
             VFloat rsat[2] = { vload(sat.r), vload(sat.r+4) };
             VFloat gsat[2] = { vload(sat.g), vload(sat.g+4) };
             VFloat bsat[2] = { vload(sat.b), vload(sat.b+4) };
             VFloat wsat[2] = { vload(sat.w), vload(sat.w+4) };
 
             __m128i idx0 = (__m128i &)s_byte_vtl_index0_3[4*i];
             x0.x = vpermute2(rsat[0], rsat[1], idx0);
             x0.y = vpermute2(gsat[0], gsat[1], idx0);
             x0.z = vpermute2(bsat[0], bsat[1], idx0);
             w0   = vpermute2(wsat[0], wsat[1], idx0);
 
             __m128i idx1 = (__m128i &)s_byte_vtl_index1_3[4*i];
             x1.x = vpermute2(rsat[0], rsat[1], idx1);
             x1.y = vpermute2(gsat[0], gsat[1], idx1);
             x1.z = vpermute2(bsat[0], bsat[1], idx1);
             w1   = vpermute2(wsat[0], wsat[1], idx1);
         }
         else {
             VFloat rsat[4] = { vload(sat.r), vload(sat.r+4), vload(sat.r+8), vload(sat.r+12) };
             VFloat gsat[4] = { vload(sat.g), vload(sat.g+4), vload(sat.g+8), vload(sat.g+12) };
             VFloat bsat[4] = { vload(sat.b), vload(sat.b+4), vload(sat.b+8), vload(sat.b+12) };
             VFloat wsat[4] = { vload(sat.w), vload(sat.w+4), vload(sat.w+8), vload(sat.w+12) };
 
             __m128i idx0 = (__m128i &)s_byte_vtl_index0_3[4*i];
             x0.x = vpermute4(rsat[0], rsat[1], rsat[2], rsat[3], idx0);
             x0.y = vpermute4(gsat[0], gsat[1], gsat[2], gsat[3], idx0);
             x0.z = vpermute4(bsat[0], bsat[1], bsat[2], bsat[3], idx0);
             w0   = vpermute4(wsat[0], wsat[1], wsat[2], wsat[3], idx0);
 
             __m128i idx1 = (__m128i &)s_byte_vtl_index1_3[4*i];
             x1.x = vpermute4(rsat[0], rsat[1], rsat[2], rsat[3], idx1);
             x1.y = vpermute4(gsat[0], gsat[1], gsat[2], gsat[3], idx1);
             x1.z = vpermute4(bsat[0], bsat[1], bsat[2], bsat[3], idx1);
             w1   = vpermute4(wsat[0], wsat[1], wsat[2], wsat[3], idx1);
         }
 
 #elif ICBC_USE_NEON_VTL
 
         uint8x16_t idx0 = (uint8x16_t &)s_byte_vtl_index0_3[4*i];
         uint8x16_t idx1 = (uint8x16_t &)s_byte_vtl_index1_3[4*i];
 
         /*if (count <= 4) {
             uint8x16_t rsat1 = (uint8x16_t&)sat.r;
             uint8x16_t gsat1 = (uint8x16_t&)sat.g;
             uint8x16_t bsat1 = (uint8x16_t&)sat.b;
             uint8x16_t wsat1 = (uint8x16_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl1q_u8(rsat1, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl1q_u8(gsat1, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl1q_u8(bsat1, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl1q_u8(wsat1, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl1q_u8(rsat1, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl1q_u8(gsat1, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl1q_u8(bsat1, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl1q_u8(wsat1, idx1));
         }
         else*/
         if (count <= 8) {
             uint8x16x2_t rsat2 = (uint8x16x2_t&)sat.r;
             uint8x16x2_t gsat2 = (uint8x16x2_t&)sat.g;
             uint8x16x2_t bsat2 = (uint8x16x2_t&)sat.b;
             uint8x16x2_t wsat2 = (uint8x16x2_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl2q_u8(rsat2, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl2q_u8(gsat2, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl2q_u8(bsat2, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl2q_u8(wsat2, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl2q_u8(rsat2, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl2q_u8(gsat2, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl2q_u8(bsat2, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl2q_u8(wsat2, idx1));
         }
         else if (count <= 12) {
             uint8x16x3_t rsat3 = (uint8x16x3_t&)sat.r;
             uint8x16x3_t gsat3 = (uint8x16x3_t&)sat.g;
             uint8x16x3_t bsat3 = (uint8x16x3_t&)sat.b;
             uint8x16x3_t wsat3 = (uint8x16x3_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl3q_u8(rsat3, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl3q_u8(gsat3, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl3q_u8(bsat3, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl3q_u8(wsat3, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl3q_u8(rsat3, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl3q_u8(gsat3, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl3q_u8(bsat3, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl3q_u8(wsat3, idx1));
         }
         else {
             uint8x16x4_t rsat4 = (uint8x16x4_t&)sat.r;
             uint8x16x4_t gsat4 = (uint8x16x4_t&)sat.g;
             uint8x16x4_t bsat4 = (uint8x16x4_t&)sat.b;
             uint8x16x4_t wsat4 = (uint8x16x4_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl4q_u8(rsat4, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl4q_u8(gsat4, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl4q_u8(bsat4, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl4q_u8(wsat4, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl4q_u8(rsat4, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl4q_u8(gsat4, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl4q_u8(bsat4, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl4q_u8(wsat4, idx1));
         }
 
 #else
 
         // Scalar path
         x0.x = vzero(); x0.y = vzero(); x0.z = vzero(); w0 = vzero();
         x1.x = vzero(); x1.y = vzero(); x1.z = vzero(); w1 = vzero();
 
         for (int l = 0; l < VEC_SIZE; l++) {
             int c0 = s_threeCluster[i + l].c0 - 1;
             if (c0 >= 0) {
                 lane(x0.x, l) = sat.r[c0];
                 lane(x0.y, l) = sat.g[c0];
                 lane(x0.z, l) = sat.b[c0];
                 lane(w0, l) = sat.w[c0];
             }
 
             int c1 = s_threeCluster[i + l].c1 - 1;
             if (c1 >= 0) {
                 lane(x1.x, l) = sat.r[c1];
                 lane(x1.y, l) = sat.g[c1];
                 lane(x1.z, l) = sat.b[c1];
                 lane(w1, l) = sat.w[c1];
             }
         }
 
 #endif
 
         VFloat w2 = vbroadcast(w_sum) - w1;
         x1 = x1 - x0;
         w1 = w1 - w0;
 
         VFloat alphabeta_sum = w1 * vbroadcast(0.25f);
         VFloat alpha2_sum = w0 + alphabeta_sum;
         VFloat beta2_sum = w2 + alphabeta_sum;
         VFloat factor = vrcp(vm2sub(alpha2_sum, beta2_sum, alphabeta_sum, alphabeta_sum));
 
         VVector3 alphax_sum = x0 + x1 * vbroadcast(0.5f);
         VVector3 betax_sum = vbroadcast(r_sum, g_sum, b_sum) - alphax_sum;
 
         VVector3 a = vm2sub(alphax_sum, beta2_sum, betax_sum, alphabeta_sum) * factor;
         VVector3 b = vm2sub(betax_sum, alpha2_sum, alphax_sum, alphabeta_sum) * factor;
 
         // snap to the grid
         a = vround_ept(a);
         b = vround_ept(b);
 
         // compute the error
         VVector3 e2 = vm2sub(a, vmsub(b, alphabeta_sum, alphax_sum), b, betax_sum) * vbroadcast(2.0f);
         VVector3 e1 = vmadd(a * a, alpha2_sum, vmadd(b * b, beta2_sum, e2));
 
         // apply the metric to the error term
         VFloat error = vdot(e1, vbroadcast(metric_sqr));
 
 
         // keep the solution if it wins
         auto mask = (error < vbesterror);
 
         // I could mask the unused lanes here, but instead I set the invalid SAT entries to FLT_MAX.
         //mask = (mask & (vbroadcast(total_order_count) >= tid8(i))); // This doesn't seem to help. Is it OK to consider elements out of bounds?
 
         vbesterror = vselect(mask, vbesterror, error);
         vbeststart = vselect(mask, vbeststart, a);
         vbestend = vselect(mask, vbestend, b);
     }
 
     int bestindex = reduce_min_index(vbesterror);
 
     start->x = lane(vbeststart.x, bestindex);
     start->y = lane(vbeststart.y, bestindex);
     start->z = lane(vbeststart.z, bestindex);
     end->x = lane(vbestend.x, bestindex);
     end->y = lane(vbestend.y, bestindex);
     end->z = lane(vbestend.z, bestindex);
 }
 
 
 static void cluster_fit_four(const SummedAreaTable & sat, int count, Vector3 metric_sqr, Vector3 * start, Vector3 * end)
 {
     const float r_sum = sat.r[count-1];
     const float g_sum = sat.g[count-1];
     const float b_sum = sat.b[count-1];
     const float w_sum = sat.w[count-1];
 
     VFloat vbesterror = vbroadcast(FLT_MAX);
     VVector3 vbeststart = { vzero(), vzero(), vzero() };
     VVector3 vbestend = { vzero(), vzero(), vzero() };
 
     // check all possible clusters for this total order
     const int total_order_count = s_fourClusterTotal[count - 1];
 
     for (int i = 0; i < total_order_count; i += VEC_SIZE)
     {
         VVector3 x0, x1, x2;
         VFloat w0, w1, w2;
 
         /*
         // Another approach would be to load and broadcast one color at a time like I do in my old CUDA implementation.
         uint akku = 0;
 
         // Compute alpha & beta for this permutation.
         #pragma unroll
         for (int i = 0; i < 16; i++)
         {
             const uint bits = permutation >> (2*i);
 
             alphax_sum += alphaTable4[bits & 3] * colors[i];
             akku += prods4[bits & 3];
         }
 
         float alpha2_sum = float(akku >> 16);
         float beta2_sum = float((akku >> 8) & 0xff);
         float alphabeta_sum = float(akku & 0xff);
         float3 betax_sum = 9.0f * color_sum - alphax_sum;
         */
 
 #if ICBC_USE_AVX512_PERMUTE
 
         auto loadmask = lane_id() < vbroadcast(float(count));
 
         // Load sat in one register:
         VFloat vrsat = vload(loadmask, sat.r, FLT_MAX);
         VFloat vgsat = vload(loadmask, sat.g, FLT_MAX);
         VFloat vbsat = vload(loadmask, sat.b, FLT_MAX);
         VFloat vwsat = vload(loadmask, sat.w, FLT_MAX);
 
         // Load 4 uint8 per lane.
         VInt packedClusterIndex = vload((int *)&s_fourCluster[i]);
 
         VInt c0 = (packedClusterIndex & 0xFF) - 1;
         VInt c1 = ((packedClusterIndex >> 8) & 0xFF) - 1;
         VInt c2 = ((packedClusterIndex >> 16)) - 1; // No need for &
 
         x0.x = vpermuteif(c0 >= 0, vrsat, c0);
         x0.y = vpermuteif(c0 >= 0, vgsat, c0);
         x0.z = vpermuteif(c0 >= 0, vbsat, c0);
         w0   = vpermuteif(c0 >= 0, vwsat, c0);
 
         x1.x = vpermuteif(c1 >= 0, vrsat, c1);
         x1.y = vpermuteif(c1 >= 0, vgsat, c1);
         x1.z = vpermuteif(c1 >= 0, vbsat, c1);
         w1   = vpermuteif(c1 >= 0, vwsat, c1);
 
         x2.x = vpermuteif(c2 >= 0, vrsat, c2);
         x2.y = vpermuteif(c2 >= 0, vgsat, c2);
         x2.z = vpermuteif(c2 >= 0, vbsat, c2);
         w2   = vpermuteif(c2 >= 0, vwsat, c2);
 
 #elif ICBC_USE_AVX2_PERMUTE2
 
         // Load 4 uint8 per lane.
         VInt packedClusterIndex = vload((int *)&s_fourCluster[i]);
 
         VInt c0 = (packedClusterIndex & 0xFF) - 1;
         VInt c1 = ((packedClusterIndex >> 8) & 0xFF) - 1;
         VInt c2 = ((packedClusterIndex >> 16)) - 1; // No need for &
 
         if (count <= 8) {
             // Load sat.r in one register:
             VFloat rLo = vload(sat.r);
             VFloat gLo = vload(sat.g);
             VFloat bLo = vload(sat.b);
             VFloat wLo = vload(sat.w);
 
             x0.x = vpermuteif(c0 >= 0, rLo, c0);
             x0.y = vpermuteif(c0 >= 0, gLo, c0);
             x0.z = vpermuteif(c0 >= 0, bLo, c0);
             w0   = vpermuteif(c0 >= 0, wLo, c0);
 
             x1.x = vpermuteif(c1 >= 0, rLo, c1);
             x1.y = vpermuteif(c1 >= 0, gLo, c1);
             x1.z = vpermuteif(c1 >= 0, bLo, c1);
             w1   = vpermuteif(c1 >= 0, wLo, c1);
 
             x2.x = vpermuteif(c2 >= 0, rLo, c2);
             x2.y = vpermuteif(c2 >= 0, gLo, c2);
             x2.z = vpermuteif(c2 >= 0, bLo, c2);
             w2   = vpermuteif(c2 >= 0, wLo, c2);
         }
         else {
             // Load sat.r in two registers:
             VFloat rLo = vload(sat.r); VFloat rHi = vload(sat.r + 8);
             VFloat gLo = vload(sat.g); VFloat gHi = vload(sat.g + 8);
             VFloat bLo = vload(sat.b); VFloat bHi = vload(sat.b + 8);
             VFloat wLo = vload(sat.w); VFloat wHi = vload(sat.w + 8);
 
             x0.x = vpermute2if(c0 >= 0, rLo, rHi, c0);
             x0.y = vpermute2if(c0 >= 0, gLo, gHi, c0);
             x0.z = vpermute2if(c0 >= 0, bLo, bHi, c0);
             w0   = vpermute2if(c0 >= 0, wLo, wHi, c0);
 
             x1.x = vpermute2if(c1 >= 0, rLo, rHi, c1);
             x1.y = vpermute2if(c1 >= 0, gLo, gHi, c1);
             x1.z = vpermute2if(c1 >= 0, bLo, bHi, c1);
             w1   = vpermute2if(c1 >= 0, wLo, wHi, c1);
 
             x2.x = vpermute2if(c2 >= 0, rLo, rHi, c2);
             x2.y = vpermute2if(c2 >= 0, gLo, gHi, c2);
             x2.z = vpermute2if(c2 >= 0, bLo, bHi, c2);
             w2   = vpermute2if(c2 >= 0, wLo, wHi, c2);
         }
 
 #elif ICBC_USE_SSSE3_SHUFFLEB
 
         if (count <= 8) {
             VFloat rsat[2] = { vload(sat.r), vload(sat.r+4) };
             VFloat gsat[2] = { vload(sat.g), vload(sat.g+4) };
             VFloat bsat[2] = { vload(sat.b), vload(sat.b+4) };
             VFloat wsat[2] = { vload(sat.w), vload(sat.w+4) };
 
             __m128i idx0 = (__m128i &)s_byte_vtl_index0_4[4*i];
             x0.x = vpermute2(rsat[0], rsat[1], idx0);
             x0.y = vpermute2(gsat[0], gsat[1], idx0);
             x0.z = vpermute2(bsat[0], bsat[1], idx0);
             w0   = vpermute2(wsat[0], wsat[1], idx0);
 
             __m128i idx1 = (__m128i &)s_byte_vtl_index1_4[4*i];
             x1.x = vpermute2(rsat[0], rsat[1], idx1);
             x1.y = vpermute2(gsat[0], gsat[1], idx1);
             x1.z = vpermute2(bsat[0], bsat[1], idx1);
             w1   = vpermute2(wsat[0], wsat[1], idx1);
 
             __m128i idx2 = (__m128i &)s_byte_vtl_index2_4[4*i];
             x2.x = vpermute2(rsat[0], rsat[1], idx2);
             x2.y = vpermute2(gsat[0], gsat[1], idx2);
             x2.z = vpermute2(bsat[0], bsat[1], idx2);
             w2   = vpermute2(wsat[0], wsat[1], idx2);
         }
         else {
             VFloat rsat[4] = { vload(sat.r), vload(sat.r+4), vload(sat.r+8), vload(sat.r+12) };
             VFloat gsat[4] = { vload(sat.g), vload(sat.g+4), vload(sat.g+8), vload(sat.g+12) };
             VFloat bsat[4] = { vload(sat.b), vload(sat.b+4), vload(sat.b+8), vload(sat.b+12) };
             VFloat wsat[4] = { vload(sat.w), vload(sat.w+4), vload(sat.w+8), vload(sat.w+12) };
 
             __m128i idx0 = (__m128i &)s_byte_vtl_index0_4[4*i];
             x0.x = vpermute4(rsat[0], rsat[1], rsat[2], rsat[3], idx0);
             x0.y = vpermute4(gsat[0], gsat[1], gsat[2], gsat[3], idx0);
             x0.z = vpermute4(bsat[0], bsat[1], bsat[2], bsat[3], idx0);
             w0   = vpermute4(wsat[0], wsat[1], wsat[2], wsat[3], idx0);
 
             __m128i idx1 = (__m128i &)s_byte_vtl_index1_4[4*i];
             x1.x = vpermute4(rsat[0], rsat[1], rsat[2], rsat[3], idx1);
             x1.y = vpermute4(gsat[0], gsat[1], gsat[2], gsat[3], idx1);
             x1.z = vpermute4(bsat[0], bsat[1], bsat[2], bsat[3], idx1);
             w1   = vpermute4(wsat[0], wsat[1], wsat[2], wsat[3], idx1);
 
             __m128i idx2 = (__m128i &)s_byte_vtl_index2_4[4*i];
             x2.x = vpermute4(rsat[0], rsat[1], rsat[2], rsat[3], idx2);
             x2.y = vpermute4(gsat[0], gsat[1], gsat[2], gsat[3], idx2);
             x2.z = vpermute4(bsat[0], bsat[1], bsat[2], bsat[3], idx2);
             w2   = vpermute4(wsat[0], wsat[1], wsat[2], wsat[3], idx2);
         }
 
 #elif ICBC_USE_NEON_VTL
 
         /*
         // c0 = { pq[12],pq[12],pq[12],pq[12], pq[8],pq[8],pq[8],pq[8], pq[4],pq[4],pq[4],pq[4], pq[0],pq[0],pq[0],pq[0] }
         // idx0 = (c0 - 1) * 4 + byte_offset
         // idx0 = c0 == 0 ? 255 : idx0
 
         // c1 = { pq[13],pq[13],pq[13],pq[13], pq[9],pq[9],pq[9],pq[9], pq[5],pq[5],pq[5],pq[5], pq[1],pq[1],pq[1],pq[1] }
         // idx1 = (c1 - 1) * 4 + byte_offset
         // idx1 = c1 == 0 ? 255 : idx1
 
         // c2 = { pq[14],pq[14],pq[14],pq[14], pq[10],pq[10],pq[10],pq[10], pq[6],pq[6],pq[6],pq[6], pq[2],pq[2],pq[2],pq[2] }
         // idx2 = (c2 - 1) * 4 + byte_offset
         // idx2 = c2 == 0 ? 255 : idx2
 
         uint8x16_t packedClusterIndex = (uint8x16_t &)s_fourCluster[i];
 
         const uint8x16_t byte_offsets = {0,1,2,3, 0,1,2,3, 0,1,2,3, 0,1,2,3};
         uint8x16_t indices = {0,0,0,0, 4,4,4,4, 8,8,8,8, 12,12,12,12};
 
         uint8x16_t c0 = vqtbl1q_u8(packedClusterIndex, indices);
         uint8x16_t idx0 = vaddq_u8(vshlq_n_u8(vsubq_u8(c0, vdupq_n_u8(1)), 2), byte_offsets); // (c0 - 1) * 4 + byte_offsets
         idx0 = vorrq_u8(idx0, vceqzq_u8(c0));  // idx0 = c0 == 0 ? 255 : idx0
 
         indices = vaddq_u8(indices, vdupq_n_u8(1));
         uint8x16_t c1 = vqtbl1q_u8(packedClusterIndex, indices);
         uint8x16_t idx1 = vaddq_u8(vshlq_n_u8(vsubq_u8(c1, vdupq_n_u8(1)), 2), byte_offsets); // (c1 - 1) * 4 + byte_offsets
         idx1 = vorrq_u8(idx1, vceqzq_u8(c1));  // idx1 = c1 == 0 ? 255 : idx1
 
         indices = vaddq_u8(indices, vdupq_n_u8(1));
         uint8x16_t c2 = vqtbl1q_u8(packedClusterIndex, indices);
         uint8x16_t idx2 = vaddq_u8(vshlq_n_u8(vsubq_u8(c2, vdupq_n_u8(1)), 2), byte_offsets); // (c2 - 1) * 4 + byte_offsets
         idx2 = vorrq_u8(idx2, vceqzq_u8(c2));  // idx2 = c2 == 0 ? 255 : bc0
         */
 
         // Loading the precomputed byte indices is faster than computing them on the fly.
         uint8x16_t idx0 = (uint8x16_t &)s_byte_vtl_index0_4[4*i];
         uint8x16_t idx1 = (uint8x16_t &)s_byte_vtl_index1_4[4*i];
         uint8x16_t idx2 = (uint8x16_t &)s_byte_vtl_index2_4[4*i];
 
         /*if (count <= 4) {
             uint8x16_t rsat1 = (uint8x16_t&)sat.r;
             uint8x16_t gsat1 = (uint8x16_t&)sat.g;
             uint8x16_t bsat1 = (uint8x16_t&)sat.b;
             uint8x16_t wsat1 = (uint8x16_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl1q_u8(rsat1, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl1q_u8(gsat1, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl1q_u8(bsat1, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl1q_u8(wsat1, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl1q_u8(rsat1, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl1q_u8(gsat1, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl1q_u8(bsat1, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl1q_u8(wsat1, idx1));
 
             x2.x = vreinterpretq_f32_u8(vqtbl1q_u8(rsat1, idx2));
             x2.y = vreinterpretq_f32_u8(vqtbl1q_u8(gsat1, idx2));
             x2.z = vreinterpretq_f32_u8(vqtbl1q_u8(bsat1, idx2));
             w2   = vreinterpretq_f32_u8(vqtbl1q_u8(wsat1, idx2));
         }
         else*/ 
         if (count <= 8) {
             uint8x16x2_t rsat2 = (uint8x16x2_t&)sat.r;
             uint8x16x2_t gsat2 = (uint8x16x2_t&)sat.g;
             uint8x16x2_t bsat2 = (uint8x16x2_t&)sat.b;
             uint8x16x2_t wsat2 = (uint8x16x2_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl2q_u8(rsat2, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl2q_u8(gsat2, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl2q_u8(bsat2, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl2q_u8(wsat2, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl2q_u8(rsat2, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl2q_u8(gsat2, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl2q_u8(bsat2, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl2q_u8(wsat2, idx1));
 
             x2.x = vreinterpretq_f32_u8(vqtbl2q_u8(rsat2, idx2));
             x2.y = vreinterpretq_f32_u8(vqtbl2q_u8(gsat2, idx2));
             x2.z = vreinterpretq_f32_u8(vqtbl2q_u8(bsat2, idx2));
             w2   = vreinterpretq_f32_u8(vqtbl2q_u8(wsat2, idx2));
         }
         else if (count <= 12) {
             uint8x16x3_t rsat3 = (uint8x16x3_t&)sat.r;
             uint8x16x3_t gsat3 = (uint8x16x3_t&)sat.g;
             uint8x16x3_t bsat3 = (uint8x16x3_t&)sat.b;
             uint8x16x3_t wsat3 = (uint8x16x3_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl3q_u8(rsat3, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl3q_u8(gsat3, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl3q_u8(bsat3, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl3q_u8(wsat3, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl3q_u8(rsat3, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl3q_u8(gsat3, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl3q_u8(bsat3, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl3q_u8(wsat3, idx1));
 
             x2.x = vreinterpretq_f32_u8(vqtbl3q_u8(rsat3, idx2));
             x2.y = vreinterpretq_f32_u8(vqtbl3q_u8(gsat3, idx2));
             x2.z = vreinterpretq_f32_u8(vqtbl3q_u8(bsat3, idx2));
             w2   = vreinterpretq_f32_u8(vqtbl3q_u8(wsat3, idx2));
         }
         else {
             uint8x16x4_t rsat4 = (uint8x16x4_t&)sat.r;
             uint8x16x4_t gsat4 = (uint8x16x4_t&)sat.g;
             uint8x16x4_t bsat4 = (uint8x16x4_t&)sat.b;
             uint8x16x4_t wsat4 = (uint8x16x4_t&)sat.w;
 
             x0.x = vreinterpretq_f32_u8(vqtbl4q_u8(rsat4, idx0));
             x0.y = vreinterpretq_f32_u8(vqtbl4q_u8(gsat4, idx0));
             x0.z = vreinterpretq_f32_u8(vqtbl4q_u8(bsat4, idx0));
             w0   = vreinterpretq_f32_u8(vqtbl4q_u8(wsat4, idx0));
 
             x1.x = vreinterpretq_f32_u8(vqtbl4q_u8(rsat4, idx1));
             x1.y = vreinterpretq_f32_u8(vqtbl4q_u8(gsat4, idx1));
             x1.z = vreinterpretq_f32_u8(vqtbl4q_u8(bsat4, idx1));
             w1   = vreinterpretq_f32_u8(vqtbl4q_u8(wsat4, idx1));
 
             x2.x = vreinterpretq_f32_u8(vqtbl4q_u8(rsat4, idx2));
             x2.y = vreinterpretq_f32_u8(vqtbl4q_u8(gsat4, idx2));
             x2.z = vreinterpretq_f32_u8(vqtbl4q_u8(bsat4, idx2));
             w2   = vreinterpretq_f32_u8(vqtbl4q_u8(wsat4, idx2));
         }
 
 #else
 
         // Scalar path
         x0.x = vzero(); x0.y = vzero(); x0.z = vzero(); w0 = vzero();
         x1.x = vzero(); x1.y = vzero(); x1.z = vzero(); w1 = vzero();
         x2.x = vzero(); x2.y = vzero(); x2.z = vzero(); w2 = vzero();
 
         for (int l = 0; l < VEC_SIZE; l++) {
             int c0 = s_fourCluster[i + l].c0 - 1;
             if (c0 >= 0) {
                 lane(x0.x, l) = sat.r[c0];
                 lane(x0.y, l) = sat.g[c0];
                 lane(x0.z, l) = sat.b[c0];
                 lane(w0, l) = sat.w[c0];
             }
 
             int c1 = s_fourCluster[i + l].c1 - 1;
             if (c1 >= 0) {
                 lane(x1.x, l) = sat.r[c1];
                 lane(x1.y, l) = sat.g[c1];
                 lane(x1.z, l) = sat.b[c1];
                 lane(w1, l) = sat.w[c1];
             }
 
             int c2 = s_fourCluster[i + l].c2 - 1;
             if (c2 >= 0) {
                 lane(x2.x, l) = sat.r[c2];
                 lane(x2.y, l) = sat.g[c2];
                 lane(x2.z, l) = sat.b[c2];
                 lane(w2, l) = sat.w[c2];
             }
         }
 
 #endif
 
         VFloat w3 = vbroadcast(w_sum) - w2;
         x2 = x2 - x1;
         x1 = x1 - x0;
         w2 = w2 - w1;
         w1 = w1 - w0;
 
 #if 1
         VFloat alpha2_sum = vmadd(w2, w1_9, vmadd(w1, w4_9, w0));
         VFloat beta2_sum  = vmadd(w1, w1_9, vmadd(w2, w4_9, w3));
 
         VFloat alphabeta_sum = (w1 + w2) * vbroadcast(w2_9);
         VFloat factor = vrcp(vm2sub(alpha2_sum, beta2_sum, alphabeta_sum, alphabeta_sum));
 
         VVector3 alphax_sum = vmadd(x2, w1_3, vmadd(x1, w2_3, x0));
 #else
         VFloat alpha2_sum = vmadd(w2, (1.0f / 9.0f), vmadd(w1, (4.0f / 9.0f), w0));
         VFloat beta2_sum  = vmadd(w1, (1.0f / 9.0f), vmadd(w2, (4.0f / 9.0f), w3));
 
         VFloat alphabeta_sum = (w1 + w2) * vbroadcast(2.0f / 9.0f);
         VFloat factor = vrcp(vm2sub(alpha2_sum, beta2_sum, alphabeta_sum, alphabeta_sum));
 
         VVector3 alphax_sum = vmadd(x2, (1.0f / 3.0f), vmadd(x1, (2.0f / 3.0f), x0));
 #endif
         VVector3 betax_sum = vbroadcast(r_sum, g_sum, b_sum) - alphax_sum;
 
         VVector3 a = vm2sub(alphax_sum, beta2_sum, betax_sum, alphabeta_sum) * factor;
         VVector3 b = vm2sub(betax_sum, alpha2_sum, alphax_sum, alphabeta_sum) * factor;
 
         // snap to the grid
         a = vround_ept(a);
         b = vround_ept(b);
 
         // compute the error
         VVector3 e2 = vm2sub(a, vmsub(b, alphabeta_sum, alphax_sum), b, betax_sum) * vbroadcast(2.0f);
         VVector3 e1 = vmadd(a * a, alpha2_sum, vmadd(b * b, beta2_sum, e2));
 
         // apply the metric to the error term
         VFloat error = vdot(e1, vbroadcast(metric_sqr));
 
         // keep the solution if it wins
         auto mask = (error < vbesterror);
 
         // We could mask the unused lanes here, but instead set the invalid SAT entries to FLT_MAX.
         //mask = (mask & (vbroadcast(total_order_count) >= tid8(i))); // This doesn't seem to help. Is it OK to consider elements out of bounds?
 
         vbesterror = vselect(mask, vbesterror, error);
         vbeststart = vselect(mask, vbeststart, a);
         vbestend = vselect(mask, vbestend, b);
     }
 
     int bestindex = reduce_min_index(vbesterror);
 
     start->x = lane(vbeststart.x, bestindex);
     start->y = lane(vbeststart.y, bestindex);
     start->z = lane(vbeststart.z, bestindex);
     end->x = lane(vbestend.x, bestindex);
     end->y = lane(vbestend.y, bestindex);
     end->z = lane(vbestend.z, bestindex);
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Palette evaluation.
 
 Decoder s_decoder = Decoder_D3D10;
 
 // D3D10
 inline void evaluate_palette4_d3d10(Color32 palette[4]) {
     palette[2].r = (2 * palette[0].r + palette[1].r) / 3;
     palette[2].g = (2 * palette[0].g + palette[1].g) / 3;
     palette[2].b = (2 * palette[0].b + palette[1].b) / 3;
     palette[2].a = 0xFF;
 
     palette[3].r = (2 * palette[1].r + palette[0].r) / 3;
     palette[3].g = (2 * palette[1].g + palette[0].g) / 3;
     palette[3].b = (2 * palette[1].b + palette[0].b) / 3;
     palette[3].a = 0xFF;
 }
 inline void evaluate_palette3_d3d10(Color32 palette[4]) {
     palette[2].r = (palette[0].r + palette[1].r) / 2;
     palette[2].g = (palette[0].g + palette[1].g) / 2;
     palette[2].b = (palette[0].b + palette[1].b) / 2;
     palette[2].a = 0xFF;
     palette[3].u = 0;
 }
 static void evaluate_palette_d3d10(Color16 c0, Color16 c1, Color32 palette[4]) {
     palette[0] = bitexpand_color16_to_color32(c0);
     palette[1] = bitexpand_color16_to_color32(c1);
     if (c0.u > c1.u) {
         evaluate_palette4_d3d10(palette);
     }
     else {
         evaluate_palette3_d3d10(palette);
     }
 }
 
 // NV
 inline void evaluate_palette4_nv(Color16 c0, Color16 c1, Color32 palette[4]) {
     int gdiff = palette[1].g - palette[0].g;
     palette[2].r = uint8(((2 * c0.r + c1.r) * 22) / 8);
     palette[2].g = uint8((256 * palette[0].g + gdiff / 4 + 128 + gdiff * 80) / 256);
     palette[2].b = uint8(((2 * c0.b + c1.b) * 22) / 8);
     palette[2].a = 0xFF;
 
     palette[3].r = uint8(((2 * c1.r + c0.r) * 22) / 8);
     palette[3].g = uint8((256 * palette[1].g - gdiff / 4 + 128 - gdiff * 80) / 256);
     palette[3].b = uint8(((2 * c1.b + c0.b) * 22) / 8);
     palette[3].a = 0xFF;
 }
 inline void evaluate_palette3_nv(Color16 c0, Color16 c1, Color32 palette[4]) {
     int gdiff = palette[1].g - palette[0].g;
     palette[2].r = uint8(((c0.r + c1.r) * 33) / 8);
     palette[2].g = uint8((256 * palette[0].g + gdiff / 4 + 128 + gdiff * 128) / 256);
     palette[2].b = uint8(((c0.b + c1.b) * 33) / 8);
     palette[2].a = 0xFF;
     palette[3].u = 0;
 }
 static void evaluate_palette_nv(Color16 c0, Color16 c1, Color32 palette[4]) {
     palette[0] = bitexpand_color16_to_color32(c0);
     palette[1] = bitexpand_color16_to_color32(c1);
 
     if (c0.u > c1.u) {
         evaluate_palette4_nv(c0, c1, palette);
     }
     else {
         evaluate_palette3_nv(c0, c1, palette);
     }
 }
 
 // AMD
 inline void evaluate_palette4_amd(Color32 palette[4]) {
     palette[2].r = uint8((43 * palette[0].r + 21 * palette[1].r + 32) >> 6);
     palette[2].g = uint8((43 * palette[0].g + 21 * palette[1].g + 32) >> 6);
     palette[2].b = uint8((43 * palette[0].b + 21 * palette[1].b + 32) >> 6);
     palette[2].a = 0xFF;
 
     palette[3].r = uint8((43 * palette[1].r + 21 * palette[0].r + 32) >> 6);
     palette[3].g = uint8((43 * palette[1].g + 21 * palette[0].g + 32) >> 6);
     palette[3].b = uint8((43 * palette[1].b + 21 * palette[0].b + 32) >> 6);
     palette[3].a = 0xFF;
 }
 inline void evaluate_palette3_amd(Color32 palette[4]) {
     palette[2].r = uint8((palette[0].r + palette[1].r + 1) / 2);
     palette[2].g = uint8((palette[0].g + palette[1].g + 1) / 2);
     palette[2].b = uint8((palette[0].b + palette[1].b + 1) / 2);
     palette[2].a = 0xFF;
     palette[3].u = 0;
 }
 static void evaluate_palette_amd(Color16 c0, Color16 c1, Color32 palette[4]) {
     palette[0] = bitexpand_color16_to_color32(c0);
     palette[1] = bitexpand_color16_to_color32(c1);
 
     if (c0.u > c1.u) {
         evaluate_palette4_amd(palette);
     }
     else {
         evaluate_palette3_amd(palette);
     }
 }
 
 // Intel
 inline void evaluate_palette4_intel(Color32 palette[4]) {
     palette[2].r = uint8((171 * palette[0].r + 85 * palette[1].r + 128) >> 8);
     palette[2].g = uint8((171 * palette[0].g + 85 * palette[1].g + 128) >> 8);
     palette[2].b = uint8((171 * palette[0].b + 85 * palette[1].b + 128) >> 8);
     palette[2].a = 0xFF;
 
     palette[3].r = uint8((171 * palette[1].r + 85 * palette[0].r + 128) >> 8);
     palette[3].g = uint8((171 * palette[1].g + 85 * palette[0].g + 128) >> 8);
     palette[3].b = uint8((171 * palette[1].b + 85 * palette[0].b + 128) >> 8);
     palette[3].a = 0xFF;
 }
 inline void evaluate_palette3_intel(Color32 palette[4]) {
     palette[2].r = uint8((128 * palette[0].r + 128 * palette[1].r + 128) >> 8);
     palette[2].g = uint8((128 * palette[0].g + 128 * palette[1].g + 128) >> 8);
     palette[2].b = uint8((128 * palette[0].b + 128 * palette[1].b + 128) >> 8);
     palette[2].a = 0xFF;
     palette[3].u = 0;
 }
 static void evaluate_palette_intel(Color16 c0, Color16 c1, Color32 palette[4]) {
     palette[0] = bitexpand_color16_to_color32(c0);
     palette[1] = bitexpand_color16_to_color32(c1);
 
     if (c0.u > c1.u) {
         evaluate_palette4_intel(palette);
     }
     else {
         evaluate_palette3_intel(palette);
     }
 }
 
 inline void evaluate_palette(Color16 c0, Color16 c1, Color32 palette[4], Decoder decoder = s_decoder) {
     if (decoder == Decoder_D3D10)         evaluate_palette_d3d10(c0, c1, palette);
     else if (decoder == Decoder_NVIDIA)   evaluate_palette_nv(c0, c1, palette);
     else if (decoder == Decoder_AMD)      evaluate_palette_amd(c0, c1, palette);
     else if (decoder == Decoder_Intel)    evaluate_palette_intel(c0, c1, palette);
 }
 
 inline void evaluate_palette4(Color16 c0, Color16 c1, Color32 palette[4], Decoder decoder = s_decoder) {
     palette[0] = bitexpand_color16_to_color32(c0);
     palette[1] = bitexpand_color16_to_color32(c1);
 
     if (decoder == Decoder_D3D10)         evaluate_palette4_d3d10(palette);
     else if (decoder == Decoder_NVIDIA)   evaluate_palette4_nv(c0, c1, palette);
     else if (decoder == Decoder_AMD)      evaluate_palette4_amd(palette);
     else if (decoder == Decoder_Intel)    evaluate_palette4_intel(palette);
 }
 
 static void evaluate_palette(Color16 c0, Color16 c1, Vector3 palette[4]) {
     Color32 palette32[4];
     evaluate_palette(c0, c1, palette32);
 
     for (int i = 0; i < 4; i++) {
         palette[i] = color_to_vector3(palette32[i]);
     }
 }
 
 static void evaluate_palette(uint8 alpha0, uint8  alpha1, uint8 alpha_palette[8], Decoder decoder) {
     // @@ 
     alpha_palette[0] = alpha0;
     alpha_palette[1] = alpha1;
     if (alpha0 > alpha1) {
         // @@ 
     }
     else {
         // @@
     }
 }
 
 
 static void decode_bc1(const BlockBC1 * block, unsigned char rgba_block[16 * 4], Decoder decoder) {
     Color32 palette[4];
     evaluate_palette(block->col0, block->col1, palette, decoder);
 
     for (int i = 0; i < 16; i++) {
         int index = (block->indices >> (2 * i)) & 3;
         Color32 c = palette[index];
         rgba_block[4 * i + 0] = c.r;
         rgba_block[4 * i + 1] = c.g;
         rgba_block[4 * i + 2] = c.b;
         rgba_block[4 * i + 3] = c.a;
     }
 }
  
 static void decode_bc3(const BlockBC3 * block, unsigned char rgba_block[16 * 4], Decoder decoder) {
     Color32 rgb_palette[4];
     evaluate_palette4(block->rgb.col0, block->rgb.col1, rgb_palette, decoder);
 
     uint8 alpha_palette[8];
     evaluate_palette(block->alpha.alpha0, block->alpha.alpha1, alpha_palette, decoder);
 
     for (int i = 0; i < 16; i++) {
         int rgb_index = (block->rgb.indices >> (2 * i)) & 3;
         Color32 c = rgb_palette[rgb_index];
         rgba_block[4 * i + 0] = c.r;
         rgba_block[4 * i + 1] = c.g;
         rgba_block[4 * i + 2] = c.b;
 
         int alpha_index = 0; // @@
         uint8 alpha = alpha_palette[alpha_index];
         rgba_block[4 * i + 3] = alpha;
     }
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Error evaluation.
 
 // Different ways of estimating the error.
 
 static float evaluate_mse(const Vector3 & p, const Vector3 & c, const Vector3 & w) {
     Vector3 d = (p - c) * w * 255;
     return dot(d, d);
 }
 
 static float evaluate_mse(const Color32 & p, const Vector3 & c, const Vector3 & w) {
     Vector3 d = (color_to_vector3(p) - c) * w * 255;
     return dot(d, d);
 }
 
 static int evaluate_mse(const Color32 & p, const Color32 & c) {
     return (square(int(p.r)-c.r) + square(int(p.g)-c.g) + square(int(p.b)-c.b));
 }
 
 static int evaluate_mse(const Color32 palette[4], const Color32 & c) {
     int e0 = evaluate_mse(palette[0], c);
     int e1 = evaluate_mse(palette[1], c);
     int e2 = evaluate_mse(palette[2], c);
     int e3 = evaluate_mse(palette[3], c);
     return min(min(e0, e1), min(e2, e3));
 }
 
 // Returns MSE error in [0-255] range.
 static int evaluate_mse(const BlockBC1 * output, Color32 color, int index) {
     Color32 palette[4];
     evaluate_palette(output->col0, output->col1, palette);
 
     return evaluate_mse(palette[index], color);
 }
 
 // Returns weighted MSE error in [0-255] range.
 static float evaluate_palette_error(Color32 palette[4], const Color32 * colors, const float * weights, int count) {
 
     float total = 0.0f;
     for (int i = 0; i < count; i++) {
         total += weights[i] * float(evaluate_mse(palette, colors[i]));
     }
 
     return total;
 }
 
 static float evaluate_palette_error(Color32 palette[4], const Color32 * colors, int count) {
 
     float total = 0.0f;
     for (int i = 0; i < count; i++) {
         total += float(evaluate_mse(palette, colors[i]));
     }
 
     return total;
 }
 
 static float evaluate_mse(const Vector4 input_colors[16], const float input_weights[16], const Vector3 & color_weights, const BlockBC1 * output) {
     Color32 palette[4];
     evaluate_palette(output->col0, output->col1, palette);
 
     // evaluate error for each index.
     float error = 0.0f;
     for (int i = 0; i < 16; i++) {
         int index = (output->indices >> (2 * i)) & 3;
         error += input_weights[i] * evaluate_mse(palette[index], input_colors[i].xyz, color_weights);
     }
     return error;
 }
 
 static float evaluate_mse(const Vector4 input_colors[16], const float input_weights[16], const Vector3& color_weights, Vector3 palette[4], uint32 indices) {
 
     // evaluate error for each index.
     float error = 0.0f;
     for (int i = 0; i < 16; i++) {
         int index = (indices >> (2 * i)) & 3;
         error += input_weights[i] * evaluate_mse(palette[index], input_colors[i].xyz, color_weights);
     }
     return error;
 }
 
 float evaluate_bc1_error(const uint8 rgba_block[16*4], const BlockBC1 * block, Decoder decoder) {
     Color32 palette[4];
     evaluate_palette(block->col0, block->col1, palette, decoder);
 
     // evaluate error for each index.
     float error = 0.0f;
     for (int i = 0; i < 16; i++) {
         int index = (block->indices >> (2 * i)) & 3;
         Color32 c;
         c.r = rgba_block[4 * i + 0];
         c.g = rgba_block[4 * i + 1];
         c.b = rgba_block[4 * i + 2];
         c.a = 255;
         error += float(evaluate_mse(palette[index], c));
     }
     return error;
 }
 
 float evaluate_bc3_error(const uint8 rgba_block[16*4], const BlockBC3 * block, bool alpha_blend, Decoder decoder) {
     Color32 rgb_palette[4];
     evaluate_palette4(block->rgb.col0, block->rgb.col1, rgb_palette, decoder);
 
     uint8 alpha_palette[8];
     evaluate_palette(block->alpha.alpha0, block->alpha.alpha1, alpha_palette, decoder);
 
     // evaluate error for each index.
     float error = 0.0f;
     for (int i = 0; i < 16; i++) {
         Color32 c;
         c.r = rgba_block[4 * i + 0];
         c.g = rgba_block[4 * i + 1];
         c.b = rgba_block[4 * i + 2];
         c.a = rgba_block[4 * i + 3];
 
         int rgb_index = (block->rgb.indices >> (2 * i)) & 3;
         int alpha_index = 0; // @@ 
 
         error += (alpha_blend ? (c.a * 1.0f/255.0f) : 1.0f) * float(evaluate_mse(rgb_palette[rgb_index], c));
         error += float(square(int(alpha_palette[alpha_index]) - c.a));
     }
     return error;
 }
 
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Index selection
 
 // @@ Can we interleave the two uint16 at once?
 inline uint32 interleave_uint16_with_zeros(uint32 input)  {
     uint32 word = input;
     word = (word ^ (word << 8 )) & 0x00ff00ff;
     word = (word ^ (word << 4 )) & 0x0f0f0f0f;
     word = (word ^ (word << 2 )) & 0x33333333;
     word = (word ^ (word << 1 )) & 0x55555555;
     return word;
 }
 
 // Interleave the bits. https://lemire.me/blog/2018/01/08/how-fast-can-you-bit-interleave-32-bit-integers/
 ICBC_FORCEINLINE uint32 interleave(uint32 a, uint32 b) {
 #if ICBC_BMI2
     return _pdep_u32(a, 0x55555555) | _pdep_u32(b, 0xaaaaaaaa);
 #else
     return interleave_uint16_with_zeros(a) | (interleave_uint16_with_zeros(b) << 1);
 #endif
 }
 
 static uint compute_indices4(const Vector4 input_colors[16], const Vector3 & color_weights, const Vector3 palette[4]) {
     uint indices0 = 0;
     uint indices1 = 0;
 
     VVector3 vw = vbroadcast(color_weights);
     VVector3 vp0 = vbroadcast(palette[0]) * vw;
     VVector3 vp1 = vbroadcast(palette[1]) * vw;
     VVector3 vp2 = vbroadcast(palette[2]) * vw;
     VVector3 vp3 = vbroadcast(palette[3]) * vw;
 
     for (int i = 0; i < 16; i += VEC_SIZE) {
         VVector3 vc = vload(&input_colors[i]) * vw;
 
         VFloat d0 = vlen2(vc - vp0);
         VFloat d1 = vlen2(vc - vp1);
         VFloat d2 = vlen2(vc - vp2);
         VFloat d3 = vlen2(vc - vp3);
 
         VMask b1 = d1 > d2;
         VMask b2 = d0 > d2;
         VMask x0 = b1 & b2;
 
         VMask b0 = d0 > d3;
         VMask b3 = d1 > d3;
         x0 = x0 | (b0 & b3);
 
         VMask b4 = d2 > d3;
         VMask x1 = b0 & b4;
 
         indices0 |= mask(x0) << i;
         indices1 |= mask(x1) << i;
     }
 
     return interleave(indices1, indices0);
 }
 
 static uint compute_indices3(const Vector4 input_colors[16], const Vector3 & color_weights, bool allow_transparent_black, const Vector3 palette[4]) {
     uint indices0 = 0;
     uint indices1 = 0;
 
     VVector3 vw = vbroadcast(color_weights);
     VVector3 vp0 = vbroadcast(palette[0]) * vw;
     VVector3 vp1 = vbroadcast(palette[1]) * vw;
     VVector3 vp2 = vbroadcast(palette[2]) * vw;
 
     if (allow_transparent_black) {
         for (int i = 0; i < 16; i += VEC_SIZE) {
             VVector3 vc = vload(&input_colors[i]) * vw;
 
             VFloat d0 = vlen2(vp0 - vc);
             VFloat d1 = vlen2(vp1 - vc);
             VFloat d2 = vlen2(vp2 - vc);
             VFloat d3 = vdot(vc, vc);
 
             VMask i1 = (d1 < d2);
             VMask i2 = (d2 <= d0) & (d2 <= d1);
             VMask i3 = (d3 <= d0) & (d3 <= d1) & (d3 <= d2);
 
             indices0 |= mask(i2 | i3) << i;
             indices1 |= mask(i1 | i3) << i;
         }
     }
     else {
         for (int i = 0; i < 16; i += VEC_SIZE) {
             VVector3 vc = vload(&input_colors[i]) * vw;
 
             VFloat d0 = vlen2(vc - vp0);
             VFloat d1 = vlen2(vc - vp1);
             VFloat d2 = vlen2(vc - vp2);
 
             VMask i1 = (d1 < d2);
             VMask i2 = (d2 <= d0) & (d2 <= d1);
 
             indices0 |= mask(i2) << i;
             indices1 |= mask(i1) << i;
         }
     }
 
     return interleave(indices1, indices0);
 }
 
 
 static uint compute_indices(const Vector4 input_colors[16], const Vector3 & color_weights, const Vector3 palette[4]) {
 #if 0
     Vector3 p0 = palette[0] * color_weights;
     Vector3 p1 = palette[1] * color_weights;
     Vector3 p2 = palette[2] * color_weights;
     Vector3 p3 = palette[3] * color_weights;
 
     uint indices = 0;
     for (int i = 0; i < 16; i++) {
         Vector3 ci = input_colors[i].xyz * color_weights;
         float d0 = lengthSquared(p0 - ci);
         float d1 = lengthSquared(p1 - ci);
         float d2 = lengthSquared(p2 - ci);
         float d3 = lengthSquared(p3 - ci);
 
         uint index;
         if (d0 < d1 && d0 < d2 && d0 < d3) index = 0;
         else if (d1 < d2 && d1 < d3) index = 1;
         else if (d2 < d3) index = 2;
         else index = 3;
 
         indices |= index << (2 * i);
     }
 
     return indices;
 #else
     uint indices0 = 0;
     uint indices1 = 0;
 
     VVector3 vw = vbroadcast(color_weights);
     VVector3 vp0 = vbroadcast(palette[0]) * vw;
     VVector3 vp1 = vbroadcast(palette[1]) * vw;
     VVector3 vp2 = vbroadcast(palette[2]) * vw;
     VVector3 vp3 = vbroadcast(palette[3]) * vw;
 
     for (int i = 0; i < 16; i += VEC_SIZE) {
         VVector3 vc = vload(&input_colors[i]) * vw;
 
         VFloat d0 = vlen2(vc - vp0);
         VFloat d1 = vlen2(vc - vp1);
         VFloat d2 = vlen2(vc - vp2);
         VFloat d3 = vlen2(vc - vp3);
 
         //VMask i0 = (d0 < d1) & (d0 < d2) & (d0 < d3);
         VMask i1 = (d1 <= d0) & (d1 < d2) & (d1 < d3);
         VMask i2 = (d2 <= d0) & (d2 <= d1) & (d2 < d3);
         VMask i3 = (d3 <= d0) & (d3 <= d1) & (d3 <= d2);
         //VFloat vindex = vselect(i0, vselect(i1, vselect(i2, vbroadcast(3), vbroadcast(2)), vbroadcast(1)), vbroadcast(0));
 
         indices0 |= mask(i2 | i3) << i;
         indices1 |= mask(i1 | i3) << i;
     }
 
     uint indices = interleave(indices1, indices0);
     return indices;
 #endif
 }
 
 
 static float output_block3(const Vector4 input_colors[16], const float input_weights[16], const Vector3 & color_weights, bool allow_transparent_black, const Vector3 & v0, const Vector3 & v1, BlockBC1 * block)
 {
     Color16 color0 = vector3_to_color16(v0);
     Color16 color1 = vector3_to_color16(v1);
 
     if (color0.u > color1.u) {
         swap(color0, color1);
     }
 
     Vector3 palette[4];
     evaluate_palette(color0, color1, palette);
 
     block->col0 = color0;
     block->col1 = color1;
     block->indices = compute_indices3(input_colors, color_weights, allow_transparent_black, palette);
 
     return evaluate_mse(input_colors, input_weights, color_weights, palette, block->indices);
 }
 
 static float output_block4(const Vector4 input_colors[16], const float input_weights[16], const Vector3 & color_weights, const Vector3 & v0, const Vector3 & v1, BlockBC1 * block)
 {
     Color16 color0 = vector3_to_color16(v0);
     Color16 color1 = vector3_to_color16(v1);
 
     if (color0.u < color1.u) {
         swap(color0, color1);
     }
 
     Vector3 palette[4];
     evaluate_palette(color0, color1, palette);
 
     block->col0 = color0;
     block->col1 = color1;
     block->indices = compute_indices4(input_colors, color_weights, palette);
 
     return evaluate_mse(input_colors, input_weights, color_weights, palette, block->indices);
 }
 
 
 // Least squares fitting of color end points for the given indices. @@ Take weights into account.
 static bool optimize_end_points4(uint indices, const Vector4 * colors, /*const float * weights,*/ int count, Vector3 * a, Vector3 * b)
 {
     float alpha2_sum = 0.0f;
     float beta2_sum = 0.0f;
     float alphabeta_sum = 0.0f;
     Vector3 alphax_sum = { 0,0,0 };
     Vector3 betax_sum = { 0,0,0 };
 
     for (int i = 0; i < count; i++)
     {
         const uint bits = indices >> (2 * i);
 
         float beta = float(bits & 1);
         if (bits & 2) beta = (1 + beta) / 3.0f;
         float alpha = 1.0f - beta;
 
         alpha2_sum += alpha * alpha;
         beta2_sum += beta * beta;
         alphabeta_sum += alpha * beta;
         alphax_sum += alpha * colors[i].xyz;
         betax_sum += beta * colors[i].xyz;
     }
 
     float denom = alpha2_sum * beta2_sum - alphabeta_sum * alphabeta_sum;
     if (equal(denom, 0.0f)) return false;
 
     float factor = 1.0f / denom;
 
     *a = saturate((alphax_sum * beta2_sum - betax_sum * alphabeta_sum) * factor);
     *b = saturate((betax_sum * alpha2_sum - alphax_sum * alphabeta_sum) * factor);
 
     return true;
 }
 
 // find minimum and maximum colors based on bounding box in color space
 inline static void fit_colors_bbox(const Vector3 * colors, int count, Vector3 * __restrict c0, Vector3 * __restrict c1)
 {
     *c0 = { 0,0,0 };
     *c1 = { 1,1,1 };
 
     for (int i = 0; i < count; i++) {
         *c0 = max(*c0, colors[i]);
         *c1 = min(*c1, colors[i]);
     }
 }
 
 inline static void select_diagonal(const Vector3 * colors, int count, Vector3 * __restrict c0, Vector3 * __restrict c1)
 {
     Vector3 center = (*c0 + *c1) * 0.5f;
 
     /*Vector3 center = colors[0];
     for (int i = 1; i < count; i++) {
         center = center * float(i-1) / i + colors[i] / i;
     }*/
     /*Vector3 center = colors[0];
     for (int i = 1; i < count; i++) {
         center += colors[i];
     }
     center /= count;*/
 
     float cov_xz = 0.0f;
     float cov_yz = 0.0f;
     for (int i = 0; i < count; i++) {
         Vector3 t = colors[i] - center;
         cov_xz += t.x * t.z;
         cov_yz += t.y * t.z;
     }
 
     float x0 = c0->x;
     float y0 = c0->y;
     float x1 = c1->x;
     float y1 = c1->y;
 
     if (cov_xz < 0) {
         swap(x0, x1);
     }
     if (cov_yz < 0) {
         swap(y0, y1);
     }
 
     *c0 = { x0, y0, c0->z };
     *c1 = { x1, y1, c1->z };
 }
 
 inline static void inset_bbox(Vector3 * __restrict c0, Vector3 * __restrict c1)
 {
     float bias = (8.0f / 255.0f) / 16.0f;
     Vector3 inset = saturate((*c0 - *c1) / 16.0f - scalar_to_vector3(bias));
     *c0 = saturate(*c0 - inset);
     *c1 = saturate(*c1 + inset);
 }
 
 
 
 // Single color lookup tables from:
 // https://github.com/nothings/stb/blob/master/stb_dxt.h
 static uint8 s_match5[256][2];
 static uint8 s_match6[256][2];
 
 static void PrepareOptTable5(uint8 * table, Decoder decoder)
 {
     uint8 expand[32];
     for (int i = 0; i < 32; i++) expand[i] = uint8((i << 3) | (i >> 2));
 
     for (int i = 0; i < 256; i++) {
         int bestErr = 256 * 100;
 
         for (int mn = 0; mn < 32; mn++) {
             for (int mx = 0; mx < 32; mx++) {
                 int mine = expand[mn];
                 int maxe = expand[mx];
 
                 int err;
 
                 int amd_r = (43 * maxe + 21 * mine + 32) >> 6;
                 int amd_err = abs(amd_r - i);
 
                 int nv_r = ((2 * mx + mn) * 22) / 8;
                 int nv_err = abs(nv_r - i);
 
                 int intel_r = (171 * maxe + 85 * mine + 128) >> 8;
                 int intel_err = abs(intel_r - i);
 
                 if (decoder == Decoder_D3D10) {
                     // DX10 spec says that interpolation must be within 3% of "correct" result,
                     // add this as error term. (normally we'd expect a random distribution of
                     // +-1.5% error, but nowhere in the spec does it say that the error has to be
                     // unbiased - better safe than sorry).
                     int r = (maxe * 2 + mine) / 3;
                     err = abs(r - i) * 100 + abs(mx - mn) * 3;
 
                     // Another approach is to consider the worst of AMD and NVIDIA errors.
                     err = max(max(amd_err, nv_err), intel_err);
                 }
                 else if (decoder == Decoder_NVIDIA) {
                     err = nv_err;
                 }
                 else if (decoder == Decoder_Intel) {
                     err = intel_err;
                 }
                 else /*if (decoder == Decoder_AMD)*/ {
                     err = amd_err;
                 }
 
                 if (err < bestErr) {
                     bestErr = err;
                     table[i * 2 + 0] = uint8(mx);
                     table[i * 2 + 1] = uint8(mn);
                 }
             }
         }
     }
 }
 
 static void PrepareOptTable6(uint8 * table, Decoder decoder)
 {
     uint8 expand[64];
     for (int i = 0; i < 64; i++) expand[i] = uint8((i << 2) | (i >> 4));
 
     for (int i = 0; i < 256; i++) {
         int bestErr = 256 * 100;
 
         for (int mn = 0; mn < 64; mn++) {
             for (int mx = 0; mx < 64; mx++) {
                 int mine = expand[mn];
                 int maxe = expand[mx];
 
                 int err;
 
                 int amd_g = (43 * maxe + 21 * mine + 32) >> 6;
                 int amd_err = abs(amd_g - i);
 
                 int nv_g = (256 * mine + (maxe - mine) / 4 + 128 + (maxe - mine) * 80) / 256;
                 int nv_err = abs(nv_g - i);
 
                 int intel_g = (171 * maxe + 85 * mine + 128) >> 8;
                 int intel_err = abs(intel_g - i);
 
                 if (decoder == Decoder_D3D10) {
                     // DX10 spec says that interpolation must be within 3% of "correct" result,
                     // add this as error term. (normally we'd expect a random distribution of
                     // +-1.5% error, but nowhere in the spec does it say that the error has to be
                     // unbiased - better safe than sorry).
                     int g = (maxe * 2 + mine) / 3;
                     err = abs(g - i) * 100 + abs(mx - mn) * 3;
 
                     // Another approach is to consider the worst of AMD and NVIDIA errors.
                     err = max(max(amd_err, nv_err), intel_err);
                 }
                 else if (decoder == Decoder_NVIDIA) {
                     err = nv_err;
                 }
                 else if (decoder == Decoder_Intel) {
                     err = intel_err;
                 }
                 else /*if (decoder == Decoder_AMD)*/ {
                     err = amd_err;
                 }
 
                 if (err < bestErr) {
                     bestErr = err;
                     table[i * 2 + 0] = uint8(mx);
                     table[i * 2 + 1] = uint8(mn);
                 }
             }
         }
     }
 }
 
 
 static void init_single_color_tables(Decoder decoder)
 {
     // Prepare single color lookup tables.
     PrepareOptTable5(&s_match5[0][0], decoder);
     PrepareOptTable6(&s_match6[0][0], decoder);
 }
 
 static void init_decoder_weights(Decoder decoder)
 {
     w1_3 = 1.0f / 3; 
     w2_3 = 2.0f / 3;
 
     if (decoder == Decoder_NVIDIA) {
         // @@ Choose weights that approximate better nvidia's decoder.
     }
     else if (decoder == Decoder_AMD) {
         w1_3 = 21.0f / 64;
         w2_3 = 43.0f / 64;
     }
     else if (decoder == Decoder_Intel) {
         w1_3 = 85.0f / 256;
         w2_3 = 171.0f / 256;
     }
 
     w1_9 = w1_3 * w1_3;
     w4_9 = w2_3 * w2_3;
     w2_9 = w1_3 * w2_3;
 }
 
 // Single color compressor, based on:
 // https://mollyrocket.com/forums/viewtopic.php?t=392
 static void compress_dxt1_single_color_optimal(Color32 c, BlockBC1 * output)
 {
     output->col0.r = s_match5[c.r][0];
     output->col0.g = s_match6[c.g][0];
     output->col0.b = s_match5[c.b][0];
     output->col1.r = s_match5[c.r][1];
     output->col1.g = s_match6[c.g][1];
     output->col1.b = s_match5[c.b][1];
     output->indices = 0xaaaaaaaa;
 
     if (output->col0.u < output->col1.u)
     {
         swap(output->col0.u, output->col1.u);
         output->indices ^= 0x55555555;
     }
 }
 
 
 static float compress_dxt1_cluster_fit(const Vector4 input_colors[16], const float input_weights[16], const Vector3 * colors, const float * weights, int count, const Vector3 & color_weights, bool three_color_mode, bool try_transparent_black, bool allow_transparent_black, BlockBC1 * output)
 {
     Vector3 metric_sqr = color_weights * color_weights;
 
     SummedAreaTable sat;
     int sat_count = compute_sat(colors, weights, count, &sat);
 
     Vector3 start, end;
     cluster_fit_four(sat, sat_count, metric_sqr, &start, &end);
 
     float best_error = output_block4(input_colors, input_weights, color_weights, start, end, output);
 
     if (three_color_mode) {
         if (try_transparent_black) {
             Vector3 tmp_colors[16];
             float tmp_weights[16];
             int tmp_count = skip_blacks(colors, weights, count, tmp_colors, tmp_weights);
             if (!tmp_count) return best_error;
 
             sat_count = compute_sat(tmp_colors, tmp_weights, tmp_count, &sat);
         }
 
         cluster_fit_three(sat, sat_count, metric_sqr, &start, &end);
 
         BlockBC1 three_color_block;
         float three_color_error = output_block3(input_colors, input_weights, color_weights, allow_transparent_black, start, end, &three_color_block);
 
         if (three_color_error < best_error) {
             best_error = three_color_error;
             *output = three_color_block;
         }
     }
 
     return best_error;
 }
 
 
 static float refine_endpoints(const Vector4 input_colors[16], const float input_weights[16], const Vector3 & color_weights, bool three_color_mode, float input_error, BlockBC1 * output) {
     // TODO:
     // - Optimize palette evaluation when updating only one channel.
     // - try all diagonals.
 
     // Things that don't help:
     // - Alternate endpoint updates.
     // - Randomize order.
     // - If one direction does not improve, test opposite direction next.
 
     static const int8 deltas[16][3] = {
         {1,0,0},
         {0,1,0},
         {0,0,1},
 
         {-1,0,0},
         {0,-1,0},
         {0,0,-1},
 
         {1,1,0},
         {1,0,1},
         {0,1,1},
 
         {-1,-1,0},
         {-1,0,-1},
         {0,-1,-1},
 
         {-1,1,0},
         //{-1,0,1},
 
         {1,-1,0},
         {0,-1,1},
 
         //{1,0,-1},
         {0,1,-1},
     };
 
     float best_error = input_error;
 
     int lastImprovement = 0;
     for (int i = 0; i < 256; i++) {
         BlockBC1 refined = *output;
         int8 delta[3] = { deltas[i % 16][0], deltas[i % 16][1], deltas[i % 16][2] };
 
         if ((i / 16) & 1) {
             refined.col0.r += delta[0];
             refined.col0.g += delta[1];
             refined.col0.b += delta[2];
         }
         else {
             refined.col1.r += delta[0];
             refined.col1.g += delta[1];
             refined.col1.b += delta[2];
         }
 
         if (!three_color_mode) {
             if (refined.col0.u == refined.col1.u) refined.col1.g += 1;
             if (refined.col0.u < refined.col1.u) swap(refined.col0.u, refined.col1.u);
         }
 
         Vector3 palette[4];
         evaluate_palette(output->col0, output->col1, palette);
 
         refined.indices = compute_indices(input_colors, color_weights, palette);
 
         float refined_error = evaluate_mse(input_colors, input_weights, color_weights, &refined);
         if (refined_error < best_error) {
             best_error = refined_error;
             *output = refined;
             lastImprovement = i;
         }
 
         // Early out if the last 32 steps didn't improve error.
         if (i - lastImprovement > 32) break;
     }
 
     return best_error;
 }
 
 struct Options {
     float threshold = 0.0f;
     bool box_fit = false;
     bool least_squares_fit = false;
     bool cluster_fit = false;
     bool cluster_fit_3 = false;
     bool cluster_fit_3_black_only = false;
     bool endpoint_refinement = false;
 };
 
 static Options setup_options(Quality level, bool enable_three_color_mode, bool enable_transparent_black) {
     Options opt;
 
     switch (level) {
         case Quality_Level0:            // Box fit only.
             opt.box_fit = true;
             opt.least_squares_fit = false;
             break;
 
         case Quality_Level1:            // Box fit + least squares fit.
             opt.box_fit = true;
             opt.least_squares_fit = true;
             break;
 
         case Quality_Level2:            // Cluster fit 4, threshold = 24.
             opt.box_fit = true;
             opt.least_squares_fit = true;
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 24;
             break;
 
         case Quality_Level3:            // Cluster fit 4, threshold = 32.
             opt.box_fit = true;
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 32;
             break;
 
         case Quality_Level4:            // Cluster fit 3+4, threshold = 48.
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 48;
             break;
 
         case Quality_Level5:            // Cluster fit 3+4, threshold = 64.
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 64;
             break;
 
         case Quality_Level6:            // Cluster fit 3+4, threshold = 96.
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 96;
             break;
 
         case Quality_Level7:            // Cluster fit 3+4, threshold = 128.
             opt.cluster_fit = true;
             opt.cluster_fit_3_black_only = enable_three_color_mode && enable_transparent_black;
             opt.threshold = 1.0f / 128;
             break;
 
         case Quality_Level8:            // Cluster fit 3+4, threshold = 256.
             opt.cluster_fit = true;
             opt.cluster_fit_3 = enable_three_color_mode;
             opt.threshold = 1.0f / 256;
             break;
 
         case Quality_Level9:           // Cluster fit 3+4, threshold = 256 + Refinement.
             opt.cluster_fit = true;
             opt.cluster_fit_3 = enable_three_color_mode;
             opt.threshold = 1.0f / 256;
             opt.endpoint_refinement = true;
             break;
     }
 
     return opt;
 }
 
 
 static float compress_bc1(Quality level, const Vector4 input_colors[16], const float input_weights[16], const Vector3 & color_weights, bool three_color_mode, bool three_color_black, BlockBC1 * output)
 {
     Options opt = setup_options(level, three_color_mode, three_color_black);
 
     Vector3 colors[16];
     float weights[16];
     bool any_black = false;
     int count;
     if (opt.cluster_fit) {
         count = reduce_colors(input_colors, input_weights, 16, opt.threshold, colors, weights, &any_black);
     }
     else {
         for (int i = 0; i < 16; i++) {
             colors[i] = input_colors[i].xyz;
         }
         count = 16;
     }
 
     if (count == 0) {
         // Output trivial block.
         output->col0.u = 0;
         output->col1.u = 0;
         output->indices = 0;
         return 0;
     }
 
     // Cluster fit cannot handle single color blocks, so encode them optimally.
     if (count == 1) {
         compress_dxt1_single_color_optimal(vector3_to_color32(colors[0]), output);
         return evaluate_mse(input_colors, input_weights, color_weights, output);
     }
 
     float error = FLT_MAX;
 
     // Quick end point selection.
     if (opt.box_fit) {
         Vector3 c0, c1;
         fit_colors_bbox(colors, count, &c0, &c1);
         inset_bbox(&c0, &c1);
         select_diagonal(colors, count, &c0, &c1);
         error = output_block4(input_colors, input_weights, color_weights, c0, c1, output);
 
         // Refine color for the selected indices.
         if (opt.least_squares_fit && optimize_end_points4(output->indices, input_colors, 16, &c0, &c1)) {
             BlockBC1 optimized_block;
             float optimized_error = output_block4(input_colors, input_weights, color_weights, c0, c1, &optimized_block);
 
             if (optimized_error < error) {
                 error = optimized_error;
                 *output = optimized_block;
             }
         }
     }
 
     if (opt.cluster_fit) {
         // @@ Use current endpoints as input for initial PCA approximation?
 
         bool use_three_color_black = any_black && three_color_black;
         bool use_three_color_mode = opt.cluster_fit_3 || (use_three_color_black && opt.cluster_fit_3_black_only);
 
         // Try cluster fit.
         BlockBC1 cluster_fit_output;
         float cluster_fit_error = compress_dxt1_cluster_fit(input_colors, input_weights, colors, weights, count, color_weights, use_three_color_mode, use_three_color_black, three_color_black, &cluster_fit_output);
         if (cluster_fit_error < error) {
             *output = cluster_fit_output;
             error = cluster_fit_error;
         }
     }
 
     if (opt.endpoint_refinement) {
         error = refine_endpoints(input_colors, input_weights, color_weights, three_color_mode, error, output);
     }
 
     return error;
 }
 
 
 // Public API
 
 void init(Decoder decoder) {
     s_decoder = decoder;
     init_single_color_tables(decoder);
     init_cluster_tables();
     init_decoder_weights(decoder);
 }
 
 void decode_bc1(const void * block, unsigned char rgba_block[16 * 4], Decoder decoder/*=Decoder_D3D10*/) {
     decode_bc1((const BlockBC1 *)block, rgba_block, decoder);
 }
 
 void decode_bc3(const void * block, unsigned char rgba_block[16 * 4], Decoder decoder/*=Decoder_D3D10*/) {
     decode_bc3((const BlockBC3 *)block, rgba_block, decoder);
 }
 
 // void decode_bc4(const void * block, bool snorm, float rgba_block[16 * 4], Decoder decoder/*=Decoder_D3D10*/) {
 //     decode_bc4((const BlockBC4 *)block, snorm, rgba_block, decoder);
 // }
 
 // void decode_bc5(const void * block, bool snorm, float rgba_block[16 * 4], Decoder decoder/*=Decoder_D3D10*/) {
 //     decode_bc5((const BlockBC5 *)block, snorm, rgba_block, decoder);
 // }
 
 
 float evaluate_bc1_error(const unsigned char rgba_block[16 * 4], const void * dxt_block, Decoder decoder/*=Decoder_D3D10*/) {
     return evaluate_bc1_error(rgba_block, (const BlockBC1 *)dxt_block, decoder);
 }
 
 float evaluate_bc3_error(const unsigned char rgba_block[16 * 4], const void * dxt_block, bool alpha_blend, Decoder decoder/*=Decoder_D3D10*/) {
     return evaluate_bc3_error(rgba_block, (const BlockBC3 *)dxt_block, alpha_blend, decoder);
 }
 
 
 float compress_bc1(Quality level, const float * input_colors, const float * input_weights, const float rgb[3], bool three_color_mode, bool three_color_black, void * output) {
     return compress_bc1(level, (Vector4*)input_colors, input_weights, { rgb[0], rgb[1], rgb[2] }, three_color_mode, three_color_black, (BlockBC1*)output);
 }
 
 float compress_bc3(Quality level, const float * input_colors, const float * input_weights, const float rgb[3], bool six_alpha_mode, void * output) {
     //return compress_bc3(level, (Vector4*)input_colors, input_weights, { rgb[0], rgb[1], rgb[2] }, six_alpha_mode, (BlockBC3*)output);
     return 0.0f; // @@ 
 }
 
 } // icbc
 
 // // Do not polute preprocessor definitions.
 // #undef ICBC_SIMD
 // #undef ICBC_ASSERT
 
 // #undef ICBC_SCALAR
 // #undef ICBC_SSE2
 // #undef ICBC_SSE41
 // #undef ICBC_AVX1
 // #undef ICBC_AVX2
 // #undef ICBC_AVX512
 // #undef ICBC_NEON
 // #undef ICBC_VMX
 
 // #undef ICBC_USE_FMA
 // #undef ICBC_USE_AVX2_PERMUTE2
 // #undef ICBC_USE_AVX512_PERMUTE
 // #undef ICBC_USE_NEON_VTL
 
 // #undef ICBC_PERFECT_ROUND
 
 #endif // ICBC_IMPLEMENTATION
 
 // Version History:
 // v1.00 - Initial release.
 // v1.01 - Added SPMD code path with AVX support.
 // v1.02 - Removed SIMD code path.
 // v1.03 - Quality levels. AVX512, Neon, Altivec, vectorized reduction and index selection.
 // v1.04 - Automatic compile-time SIMD selection. Specify hw decoder at runtime. More optimizations.
 // v1.05 - Bug fixes. Small optimizations.
 
 // Copyright (c) 2020 Ignacio Castano <castano@gmail.com>
 //
 // Permission is hereby granted, free of charge, to any person obtaining
 // a copy of this software and associated documentation files (the
 // "Software"), to	deal in the Software without restriction, including
 // without limitation the rights to use, copy, modify, merge, publish,
 // distribute, sublicense, and/or sell copies of the Software, and to
 // permit persons to whom the Software is furnished to do so, subject to
 // the following conditions:
 //
 // The above copyright notice and this permission notice shall be included
 // in all copies or substantial portions of the Software.
 //
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 // TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.