Add RGB - LogLuv conversion

* Also added missing f32x4x4 dot operators.
* Combined tone-mapping and color-space.

Author: cfnptr

include/math/color-space.hpp

@@ -23,6 +23,16 @@
 namespace math
 {
 
+/**
+ * @brief Default gamma correction value.
+ * 
+ * @details
+ * The default gamma correction value of 2.2 is chosen because it approximates the natural response curve of CRT 
+ * monitors, which were the most common display technology when gamma correction became standard. Modern displays 
+ * have inherited this standard, as it closely matches how human vision perceives brightness non-linearly.
+ */
+constexpr float defaultGamma = 2.2f;
+
 /**
  * @brief Converts linear RGBA color to the sRGB color space.
  * @param rgba target linear RGBA color
@@ -46,6 +56,65 @@ static f32x4 srgbToRgb(f32x4 sRGB) noexcept
 	return r;
 }
 
+/**
+ * @brief Applies gamma correction to the specified color.
+ * 
+ * @details
+ * Gamma correction is used to adjust the brightness of an image or display to match the 
+ * nonlinear response of display devices, such as monitors. It compensates for the fact that 
+ * displays do not linearly represent the light intensity of a color.
+ * 
+ * @param color target linear RGB color to gamma correct
+ * @param invGamma inverse gamma correction value (1.0/x)
+ */
+static f32x4 gammaCorrection(f32x4 color, float invGamma) noexcept
+{
+	return f32x4(pow(color, f32x4(invGamma)), color.getW());
+}
+/**
+ * @brief Applies gamma correction to the specified color.
+ * @details See the @ref gammaCorrection().
+ * @param color target linear RGB color to gamma correct
+ */
+static f32x4 gammaCorrection(f32x4 color) noexcept
+{
+	return f32x4(pow(color, f32x4(1.0f / defaultGamma)), color.getW());
+}
+
+/**
+ * @brief Applies gamma correction to the specified color. (Fast approximation)
+ * @details See the @ref gammaCorrection().
+ * 
+ * @param color target linear RGB color to gamma correct
+ * @param invGamma inverse gamma correction value (1.0/x)
+ */
+static f32x4 fastGammaCorrection(f32x4 color, float invGamma) noexcept
+{
+	return f32x4(fastPow(color, f32x4(invGamma)), color.getW());
+}
+/**
+ * @brief Applies gamma correction to the specified color. (Fast approximation)
+ * @details See the @ref gammaCorrection().
+ * @param color target linear RGB color to gamma correct
+ */
+static f32x4 fastGammaCorrection(f32x4 color) noexcept
+{
+	return f32x4(fastPow(color, f32x4(1.0f / defaultGamma)), color.getW());
+}
+
+/**
+ * @brief Calculates relative luminance of a color. (Rec. 709)
+ * @param x target linear color
+ */
+static float calcLum(f32x4 x) noexcept { return dot3(x, f32x4(0.2126f, 0.7152f, 0.0722f)); }
+/**
+ * @brief Calculates perceptual brightness (Luma) of a color. (Rec. 709)
+ * @param rgb target linear RGB color
+ */
+static float rgbToLuma(f32x4 rgb) noexcept { return calcLum(fastGammaCorrection(rgb)); }
+
+//**********************************************************************************************************************
+// Linear sRGB <-> CIE XYZ
 static const f32x4x4 rgbToXyzMat = f32x4x4
 (
 	0.41239079926595934f, 0.21263900587151027f, 0.01933081871559182f, 0.0f,
@@ -65,12 +134,12 @@ static const f32x4x4 xyzToRgbMat = f32x4x4
  * @brief Converts linear sRGB color to the CIE XYZ color space.
  * @param rgb target linear sRGB color
  */
-static f32x4 rgbToXyz(f32x4 rgb) noexcept { return multiply3x3(rgbToXyzMat, rgb); }
+static f32x4 rgbToXyz(f32x4 rgb) noexcept { return dot3x3(rgbToXyzMat, rgb); }
 /**
  * @brief Converts CIE XYZ color to the linear sRGB color space.
  * @param xyz target CIE XYZ color
  */
-static f32x4 xyzToRgb(f32x4 xyz) noexcept { return multiply3x3(xyzToRgbMat, xyz); }
+static f32x4 xyzToRgb(f32x4 xyz) noexcept { return dot3x3(xyzToRgbMat, xyz); }
 
 /**
  * @brief Converts CIE XYZ color to the CIE xyY color space.
@@ -102,4 +171,47 @@ static f32x4 rgbToXyy(f32x4 rgb) noexcept { return xyzToXyy(rgbToXyz(rgb)); }
  */
 static f32x4 xyyToRgb(f32x4 xyy) noexcept { return xyzToRgb(xyyToXyz(xyy)); }
 
+//**********************************************************************************************************************
+// Linear sRGB <-> LogLuv
+static const f32x4x4 rgbToLogLuvMat = f32x4x4
+(
+	0.2209f, 0.1138f, 0.0102f, 0.0f,
+	0.3390f, 0.6780f, 0.1130f, 0.0f,
+	0.4184f, 0.7319f, 0.2969f, 0.0f,
+	0.0f   , 0.0f   , 0.0f   , 0.0f
+);
+static const f32x4x4 logLuvToRgbMat = f32x4x4
+(
+	 6.0014f, -1.3320f,  0.3008f, 0.0f,
+	-2.7008f,  3.1029f, -1.0882f, 0.0f,
+	-1.7996f, -5.7721f,  5.6268f, 0.0f,
+	0.0f    ,  0.0f   ,  0.0f   , 0.0f
+);
+
+/**
+ * @brief Encodes linear RGB color (HDR) to the LogLuv format.
+ * @param rgb target linear RGB color
+ */
+static uint32 rgbToLogLuv(f32x4 rgb) noexcept
+{
+	auto luv = max(dot3x3(rgbToLogLuvMat, rgb), f32x4(1e-6f));
+	auto uv = (uint2)fma(saturate((float2)luv / luv.getZ()), float2(255.0f), float2(0.5f));
+	auto logLuv = (uint32)std::fma(saturate(std::fma(std::log2(
+		luv.getY()), 1.0f / 64.0f, 0.5f)), 65535.0f, 0.5f);
+	logLuv |= (uv.x << 24u) | (uv.y << 16u);
+	return dot3(rgb, rgb) > 0.0f ? logLuv : 0;
+}
+/**
+ * @brief Decodes linear RGB color (HDR) from the LogLuv format.
+ * @param logLuv target LogLuv encoded color
+ */
+static f32x4 logLuvToRgb(uint32 logLuv)
+{
+	if (logLuv == 0) return f32x4::zero;
+	f32x4 luv; auto uv = float2(uint2(logLuv >> 24u, logLuv >> 16u) & 255u) * (1.0f / 255.0f);
+	luv.floats.y = std::exp2(std::fma(logLuv & 65535u, (1.0f / 65535.0f) * 64.0f, -32.0f));
+	luv.floats.z = luv.floats.y / uv.y; luv.floats.x = luv.floats.z * uv.x;
+	return max(dot3x3(logLuvToRgbMat, luv), f32x4::zero);
+}
+
 } // namespace math

include/math/matrix/transform.hpp

@@ -60,7 +60,7 @@ static f32x4x4 translate(f32x4 t) noexcept
  */
 static f32x4x4 translate(const f32x4x4& m, f32x4 t) noexcept
 {
-	return f32x4x4(m.c0, m.c1, m.c2, f32x4(m.c3 + multiply3x3(m, t), m.c3.getW()));
+	return f32x4x4(m.c0, m.c1, m.c2, f32x4(m.c3 + dot3x3(m, t), m.c3.getW()));
 }
 /**
  * @brief Applies translation transformation to an object in 3D space. [r = translate(t) * m]
@@ -320,7 +320,7 @@ static quat lookAtQuat(f32x4 direction, f32x4 up = f32x4::top) noexcept
 static f32x4x4 inverseTransRot(const f32x4x4& m)
 {
 	auto t = transpose4x4(m);
-	setTranslation(t, -multiply3x3(t, getTranslation(m)));
+	setTranslation(t, -dot3x3(t, getTranslation(m)));
 	return t;
 }
 

include/math/simd/matrix/float.hpp

@@ -187,7 +187,7 @@ struct [[nodiscard]] alignas(MATH_SIMD_VECTOR_ALIGNMENT) f32x4x4
 	}
 
 	/**
-	 * @brief Calculates dot product between two SIMD matrices.
+	 * @brief Calculates 4x4 dot product between two SIMD matrices.
 	 * @param[in] m target SIMD matrix to dot by
 	 */
 	f32x4x4 operator*(const f32x4x4& m) const noexcept
@@ -214,12 +214,12 @@ struct [[nodiscard]] alignas(MATH_SIMD_VECTOR_ALIGNMENT) f32x4x4
 			result[i].data = r;
 		}
 		#else
-		result = f32x4x4((*((const float4x4*)this) * (*((const float4x4*)&m))));
+		result = f32x4x4(*((const float4x4*)this) * (*((const float4x4*)&m)));
 		#endif
 		return result;
 	}
 	/**
-	 * @brief Calculates dot product between SIMD matrix and vector.
+	 * @brief Calculates 4x4 dot product between SIMD matrix and vector.
 	 * @param v target SIMD vector to dot by
 	 */
 	f32x4 operator*(f32x4 v) const noexcept
@@ -237,7 +237,7 @@ struct [[nodiscard]] alignas(MATH_SIMD_VECTOR_ALIGNMENT) f32x4x4
 		r = vmlaq_f32(r, c3.data, vdupq_laneq_f32(v.data, 3));
 		return r;
 		#else
-		return f32x4((*((const float4x4*)this) * (*((const float4*)&v))));
+		return f32x4(*((const float4x4*)this) * (*((const float4*)&v)));
 		#endif
 	}
 
@@ -306,12 +306,38 @@ static bool isBinaryLess(const f32x4x4& a, const f32x4x4& b) noexcept
 	return memcmp(&a, &b, sizeof(f32x4x4)) < 0;
 }
 
+/**
+ * @brief Calculates 4x4 dot product between vector and SIMD matrix. (v * m)
+ *
+ * @param v target SIMD vector to use
+ * @param[in] m target SIMD matrix to dot by
+ */
+static f32x4 operator*(f32x4 v, const f32x4x4& m) noexcept
+{
+	#if defined(MATH_SIMD_SUPPORT_SSE) || defined(MATH_SIMD_SUPPORT_AVX2)
+	auto r = _mm_mul_ps(m.c0.data, v.data);
+	r = MATH_SIMD_FMA(m.c1.data, v.data, r);
+	r = MATH_SIMD_FMA(m.c2.data, v.data, r);
+	r = MATH_SIMD_FMA(m.c3.data, v.data, r);
+	return r;
+	#elif defined(MATH_SIMD_SUPPORT_NEON)
+	auto r = vmulq_f32(m.c0.data, v.data);
+	r = vmlaq_f32(r, m.c1.data, v.data);
+	r = vmlaq_f32(r, m.c2.data, v.data);
+	r = vmlaq_f32(r, m.c3.data, v.data);
+	return r;
+	#else
+	return f32x4(*((const float4*)&v) * (*((const float4x4*)this)));
+	#endif
+}
+
 /**
  * @brief Calculates 3x3 dot product between two SIMD matrices.
+ *
  * @param[in] a first SIMD matrix to use
  * @param[in] b second SIMD matrix to use
  */
-static f32x4x4 multiply3x3(const f32x4x4& a, const f32x4x4& b) noexcept
+static f32x4x4 dot3x3(const f32x4x4& a, const f32x4x4& b) noexcept
 {
 	f32x4x4 result; 
 	#if defined(MATH_SIMD_SUPPORT_SSE) || defined(MATH_SIMD_SUPPORT_AVX2)
@@ -338,28 +364,48 @@ static f32x4x4 multiply3x3(const f32x4x4& a, const f32x4x4& b) noexcept
 	return result;
 }
 /**
- * @brief Calculates 3x3 dot product between SIMD matrix and vector.
+ * @brief Calculates 3x3 dot product between SIMD matrix and vector. (m * v)
+ *
  * @param[in] m target SIMD matrix to use
  * @param v target SIMD vector to dot by
  */
-static f32x4 multiply3x3(const f32x4x4& m, f32x4 v) noexcept
+static f32x4 dot3x3(const f32x4x4& m, f32x4 v) noexcept
 {
-	f32x4 result;
 	#if defined(MATH_SIMD_SUPPORT_SSE) || defined(MATH_SIMD_SUPPORT_AVX2)
 	auto r = _mm_mul_ps(m.c0.data, _mm_shuffle_ps(v.data, v.data, _MM_SHUFFLE(0, 0, 0, 0)));
 	r = MATH_SIMD_FMA(m.c1.data, _mm_shuffle_ps(v.data, v.data, _MM_SHUFFLE(1, 1, 1, 1)), r);
 	r = MATH_SIMD_FMA(m.c2.data, _mm_shuffle_ps(v.data, v.data, _MM_SHUFFLE(2, 2, 2, 2)), r);
-	result = r;
+	return f32x4(r).swizzle<SwX, SwY, SwZ>();
 	#elif defined(MATH_SIMD_SUPPORT_NEON)
 	auto r = vmulq_f32(m.c0.data, vdupq_laneq_f32(v.data, 0));
 	r = vmlaq_f32(r, m.c1.data, vdupq_laneq_f32(v.data, 1));
 	r = vmlaq_f32(r, m.c2.data, vdupq_laneq_f32(v.data, 2));
-	result = r;
+	return f32x4(r).swizzle<SwX, SwY, SwZ>();
+	#else
+	return f32x4((float3)v * float3x3((float3)m.c0, (float3)m.c1, (float3)m.c2));
+	#endif
+}
+/**
+ * @brief Calculates 3x3 dot product between SIMD matrix and vector. (v * m)
+ *
+ * @param v target SIMD vector to use
+ * @param[in] m target SIMD matrix to dot by
+ */
+static f32x4 dot3x3(f32x4 v, const f32x4x4& m) noexcept
+{
+	#if defined(MATH_SIMD_SUPPORT_SSE) || defined(MATH_SIMD_SUPPORT_AVX2)
+	auto r = _mm_mul_ps(m.c0.data, v.data);
+	r = MATH_SIMD_FMA(m.c1.data, v.data, r);
+	r = MATH_SIMD_FMA(m.c2.data, v.data, r);
+	return f32x4(r).swizzle<SwX, SwY, SwZ>();
+	#elif defined(MATH_SIMD_SUPPORT_NEON)
+	auto r = vmulq_f32(m.c0.data, v.data);
+	r = vmlaq_f32(r, m.c1.data, v.data);
+	r = vmlaq_f32(r, m.c2.data, v.data);
+	return f32x4(r).swizzle<SwX, SwY, SwZ>();
 	#else
-	auto t = float3x3((float3)m.c0, (float3)m.c1, (float3)m.c2) * (float3)v;
-	result = f32x4(t.x, t.y, t.z, t.z);
+	return f32x4((float3)v * float3x3((float3)m.c0, (float3)m.c1, (float3)m.c2));
 	#endif
-	return result.swizzle<SwX, SwY, SwZ>();
 }
 
 /**