ConvertToPlanarRGB YUV->RGB bits=32 keep integer arith. accuracy before conversion

Author: pinterf

avs_core/convert/convert_planar.cpp

@@ -890,7 +890,15 @@ static void convert_yuv444p16_to_rgb16_c(BYTE* dstp, const BYTE* srcY, const BYT
 template<typename pixel_t, bool lessthan16bit, bool need_float_conversion>
 static void convert_yuv_to_planarrgb_uintN_c(BYTE* dstp[3], int dstPitch[3], const BYTE* srcp[3], const int srcPitch[3], int width, int height, const ConversionMatrix& matrix, const int bits_per_pixel)
 {
- auto round_mask_plus_rgb_offset_i = 4096 + (matrix.offset_rgb << 13);
+  // Float conversion extra rules
+  // - no rounding
+  // - no integer scaling back 13 bits
+  // - no clamping to bit depth limits
+  // - the 13-bit scaling factor is integrated into the bits_per_pixel shift
+
+  constexpr int ROUNDER = need_float_conversion ? 0 : (1 << 12); // 4096
+
+  auto round_mask_plus_rgb_offset_i = ROUNDER + (matrix.offset_rgb << 13);
 
   // int64_t needed for exact 16 bit pixels.
   // unlike SIMD versions which optimize it by moving into signed int16 then back to uint16, taking care of overflows.
@@ -907,7 +915,8 @@ static void convert_yuv_to_planarrgb_uintN_c(BYTE* dstp[3], int dstPitch[3], con
     max_pixel_value = (1 << bits_per_pixel) - 1;
   }
 
-  const float inv_max = 1.0f / max_pixel_value;
+  constexpr int int_arithmetic_shift = 1 << 13;
+  const float scale_f = 1.0f / static_cast<float>(int_arithmetic_shift * max_pixel_value);
 
   for (int y = 0; y < height; y++) {
     const pixel_t* srcY = reinterpret_cast<const pixel_t*>(srcp[0]);
@@ -920,17 +929,21 @@ static void convert_yuv_to_planarrgb_uintN_c(BYTE* dstp[3], int dstPitch[3], con
       const int U = static_cast<int>(srcU[x]) - half_pixel_value;
       const int V = static_cast<int>(srcV[x]) - half_pixel_value;
 
-      int b = static_cast<int>(((safe_int_t)matrix.y_b * Y + (safe_int_t)matrix.u_b * U + (safe_int_t)matrix.v_b * V + round_mask_plus_rgb_offset_i) >> 13);
-      int g = static_cast<int>(((safe_int_t)matrix.y_g * Y + (safe_int_t)matrix.u_g * U + (safe_int_t)matrix.v_g * V + round_mask_plus_rgb_offset_i) >> 13);
-      int r = static_cast<int>(((safe_int_t)matrix.y_r * Y + (safe_int_t)matrix.u_r * U + (safe_int_t)matrix.v_r * V + round_mask_plus_rgb_offset_i) >> 13);
+      int b = static_cast<int>(((safe_int_t)matrix.y_b * Y + (safe_int_t)matrix.u_b * U + (safe_int_t)matrix.v_b * V + round_mask_plus_rgb_offset_i));
+      int g = static_cast<int>(((safe_int_t)matrix.y_g * Y + (safe_int_t)matrix.u_g * U + (safe_int_t)matrix.v_g * V + round_mask_plus_rgb_offset_i));
+      int r = static_cast<int>(((safe_int_t)matrix.y_r * Y + (safe_int_t)matrix.u_r * U + (safe_int_t)matrix.v_r * V + round_mask_plus_rgb_offset_i));
 
       if constexpr (need_float_conversion) {
-        // no clamp when converting to float
-        reinterpret_cast<float*>(d0)[x] = static_cast<float>(g) * inv_max;
-        reinterpret_cast<float*>(d1)[x] = static_cast<float>(b) * inv_max;
-        reinterpret_cast<float*>(d2)[x] = static_cast<float>(r) * inv_max;
+        // no clamp when converting to float, scale back integrated into scale_f
+        reinterpret_cast<float*>(d0)[x] = static_cast<float>(g) * scale_f;
+        reinterpret_cast<float*>(d1)[x] = static_cast<float>(b) * scale_f;
+        reinterpret_cast<float*>(d2)[x] = static_cast<float>(r) * scale_f;
       }
       else {
+        // scale back
+        g = g >> 13;
+        b = b >> 13;
+        r = r >> 13;
         // Clamp
         g = clamp(g, 0, max_pixel_value);
         b = clamp(b, 0, max_pixel_value);

avs_core/convert/intel/convert_planar_avx2.cpp

@@ -375,8 +375,9 @@ void convert_yuv_to_planarrgb_uintN_avx2(BYTE* (&dstp)[3], int(&dstPitch)[3], co
     4. Output rgb offset is also added to the precalculated patch.
   */
   const __m256i half = _mm256_set1_epi16((short)(1 << (bits_per_pixel - 1)));  // 128
-  const __m256i limit = _mm256_set1_epi16((short)((1 << bits_per_pixel) - 1)); // 255
-  const __m256i offset = _mm256_set1_epi16((short)(m.offset_y));
+  const int max_pixel_value = (1 << bits_per_pixel) - 1;
+  const __m256i limit = _mm256_set1_epi16((short)max_pixel_value); // 255
+  const __m256i offset = _mm256_set1_epi16((short)m.offset_y);
 
   // to be able to use it as signed 16 bit in madd; 4096+(16<<13) would not fit into i16
   // multiplier is 4096 instead of 1:
@@ -387,19 +388,27 @@ void convert_yuv_to_planarrgb_uintN_avx2(BYTE* (&dstp)[3], int(&dstPitch)[3], co
   constexpr int ROUND_SCALE = 4096; // 1 << 12, for 13 bit integer arithmetic: "0.5"
   const __m256i m256i_round_scale = _mm256_set1_epi16(ROUND_SCALE);
 
+  // Float conversion extra rules
+  // - no rounding
+  // - no integer scaling back 13 bits
+  // - no clamping to bit depth limits
+  // - the 13-bit scaling factor is integrated into the bits_per_pixel shift
+
+  constexpr int ROUNDER = need_float_conversion ? 0 : (1 << 12); // 4096
+
   int round_mask_plus_rgb_offset_scaled_i;
   __m256i v_patch_G, v_patch_B, v_patch_R;
   __m256i sign_flip_mask = _mm256_set1_epi16((short)0x8000); // for 16 bit pivot
 
   if constexpr (lessthan16bit) {
     // 8-14 bit
-    round_mask_plus_rgb_offset_scaled_i = (4096 + (m.offset_rgb << 13)) / ROUND_SCALE;
+    round_mask_plus_rgb_offset_scaled_i = (ROUNDER + (m.offset_rgb << 13)) / ROUND_SCALE;
     v_patch_G = v_patch_B = v_patch_R = _mm256_setzero_si256(); // No patch needed
   }
   else {
     // exact 16 bit
     // keep madd simple: only handle the rounding (ROUND_SCALE * 1)
-    round_mask_plus_rgb_offset_scaled_i = 1; // effectively 1
+    round_mask_plus_rgb_offset_scaled_i = ROUNDER / ROUND_SCALE; // effectively 1 or 0 (need_float_conversion)
 
     // move BOTH the pivot and the output RGB offset to the 32-bit patch
     // Since we have to do the patching anyway, we can combine both adjustments here
@@ -413,7 +422,8 @@ void convert_yuv_to_planarrgb_uintN_avx2(BYTE* (&dstp)[3], int(&dstPitch)[3], co
   }
 
   __m256i zero = _mm256_setzero_si256();
-  const __m256 scale_f = _mm256_set1_ps(1.0f / static_cast<float>((1u << bits_per_pixel) - 1u));
+  constexpr int int_arithmetic_shift = 1 << 13;
+  const __m256 scale_f = _mm256_set1_ps(1.0f / static_cast<float>(int_arithmetic_shift * max_pixel_value));
 
   const __m256i m_uy_G = _mm256_set1_epi32(int((static_cast<uint16_t>(m.y_g) << 16) | static_cast<uint16_t>(m.u_g))); // y and u
   const __m256i m_vR_G = _mm256_set1_epi32(int((static_cast<uint16_t>(round_mask_plus_rgb_offset_scaled_i) << 16) | static_cast<uint16_t>(m.v_g))); // rounding 13 bit >> 1 and v
@@ -459,7 +469,7 @@ void convert_yuv_to_planarrgb_uintN_avx2(BYTE* (&dstp)[3], int(&dstPitch)[3], co
       // Avisynth FRAME_ALIGN == 64, so when need_float_conversion, we cannot process 32 pixels at once at the end of the line
       bool safe_last_16float_64bytes;
       if constexpr (need_float_conversion)
-        safe_last_16float_64bytes = (x + 64) <= rowsize;
+        safe_last_16float_64bytes = (x + 16 * sizeof(pixel_t)) <= rowsize;
       else
         safe_last_16float_64bytes = false;
 
@@ -508,7 +518,9 @@ void convert_yuv_to_planarrgb_uintN_avx2(BYTE* (&dstp)[3], int(&dstPitch)[3], co
           __m256i sum = _mm256_add_epi32(_mm256_madd_epi16(m_uy, uy), _mm256_madd_epi16(m_vr, vr));
           // 16-bit adjustment (signed patch, offset, output rgb offset)
           if constexpr (!lessthan16bit) sum = _mm256_add_epi32(sum, v_patch);
-          return _mm256_srai_epi32(sum, 13); // 13 bit fixed point shift
+          if constexpr (!need_float_conversion)
+            sum = _mm256_srai_epi32(sum, 13); // 13 bit fixed point shift
+          return sum;
           };
 
         __m256i res1_lo = madd_scale(uy1_lo, vr1_lo); // Pixels 0-3, 8-11

avs_core/convert/intel/convert_planar_sse.cpp

@@ -1072,8 +1072,9 @@ void convert_yuv_to_planarrgb_uintN_sse2(BYTE *(&dstp)[3], int (&dstPitch)[3], c
   */
 
   __m128i half = _mm_set1_epi16((short)(1 << (bits_per_pixel - 1)));  // 128
-  __m128i limit = _mm_set1_epi16((short)((1 << bits_per_pixel) - 1)); // 255
-  __m128i offset = _mm_set1_epi16((short)(m.offset_y));
+  const int max_pixel_value = (1 << bits_per_pixel) - 1;
+  __m128i limit = _mm_set1_epi16((short)max_pixel_value); // 255
+  __m128i offset = _mm_set1_epi16((short)m.offset_y);
 
   // to be able to use it as signed 16 bit in madd; 4096+(16<<13) would not fit into i16
   // multiplier is 4096 instead of 1:
@@ -1084,19 +1085,27 @@ void convert_yuv_to_planarrgb_uintN_sse2(BYTE *(&dstp)[3], int (&dstPitch)[3], c
   constexpr int ROUND_SCALE = 4096; // 1 << 12, for 13 bit integer arithmetic: "0.5"
   const __m128i m128i_round_scale = _mm_set1_epi16(ROUND_SCALE);
 
+  // Float conversion extra rules
+  // - no rounding
+  // - no integer scaling back 13 bits
+  // - no clamping to bit depth limits
+  // - the 13-bit scaling factor is integrated into the bits_per_pixel shift
+
+  constexpr int ROUNDER = need_float_conversion ? 0 : (1 << 12); // 4096
+
   int round_mask_plus_rgb_offset_scaled_i;
   __m128i v_patch_G, v_patch_B, v_patch_R;
   __m128i sign_flip_mask = _mm_set1_epi16((short)0x8000); // for 16 bit pivot
 
   if constexpr (lessthan16bit) {
     // 8-14 bit
-    round_mask_plus_rgb_offset_scaled_i = (4096 + (m.offset_rgb << 13)) / ROUND_SCALE;
+    round_mask_plus_rgb_offset_scaled_i = (ROUNDER + (m.offset_rgb << 13)) / ROUND_SCALE;
     v_patch_G = v_patch_B = v_patch_R = _mm_setzero_si128(); // No patch needed
   }
   else {
     // exact 16 bit
     // keep madd simple: only handle the rounding (ROUND_SCALE * 1)
-    round_mask_plus_rgb_offset_scaled_i = 1; // effectively 1
+    round_mask_plus_rgb_offset_scaled_i = ROUNDER / ROUND_SCALE; // effectively 1 or 0 (need_float_conversion)
 
     // move BOTH the pivot and the output RGB offset to the 32-bit patch
     // Since we have to do the patching anyway, we can combine both adjustments here
@@ -1110,7 +1119,8 @@ void convert_yuv_to_planarrgb_uintN_sse2(BYTE *(&dstp)[3], int (&dstPitch)[3], c
   }
 
   __m128i zero = _mm_setzero_si128();
-  const __m128 scale_f_sse2 = _mm_set1_ps(1.0f / static_cast<float>((1u << bits_per_pixel) - 1u));
+  constexpr int int_arithmetic_shift = 1 << 13;
+  const __m128 scale_f_sse2 = _mm_set1_ps(1.0f / static_cast<float>(int_arithmetic_shift * max_pixel_value));
 
   const __m128i m_uy_G = _mm_set1_epi32((static_cast<uint16_t>(m.y_g) << 16) | static_cast<uint16_t>(m.u_g));
   const __m128i m_vr_G = _mm_set1_epi32((static_cast<uint16_t>(round_mask_plus_rgb_offset_scaled_i) << 16) | static_cast<uint16_t>(m.v_g));
@@ -1170,11 +1180,13 @@ void convert_yuv_to_planarrgb_uintN_sse2(BYTE *(&dstp)[3], int (&dstPitch)[3], c
         auto process_plane_sse2 = [&](BYTE* plane_ptr, __m128i m_uy, __m128i m_vr, __m128i v_patch) {
 
           auto madd_shift = [&](__m128i uy, __m128i vr) {
-            __m128i res = _mm_add_epi32(_mm_madd_epi16(m_uy, uy), _mm_madd_epi16(m_vr, vr));
+            __m128i sum = _mm_add_epi32(_mm_madd_epi16(m_uy, uy), _mm_madd_epi16(m_vr, vr));
             // 16-bit adjustment (signed patch, offset, output rgb offset)
-            if constexpr (!lessthan16bit) res = _mm_add_epi32(res, v_patch);
-            return _mm_srai_epi32(res, 13); // 13 bit fixed point shift
-            };
+            if constexpr (!lessthan16bit) sum = _mm_add_epi32(sum, v_patch);
+            if constexpr(!need_float_conversion)
+              sum = _mm_srai_epi32(sum, 13); // 13 bit fixed point shift
+            return sum;
+          };
 
           __m128i res_lo = madd_shift(uy_lo, vr_lo); // Pixels 0, 1, 2, 3
           __m128i res_hi = madd_shift(uy_hi, vr_hi); // Pixels 4, 5, 6, 7