4
\$\begingroup\$

The cofactor of a 4×4 matrix can be used to convert a "regular geometry" matrix into the matrix that transforms the normals. It's an alternative to the common inverse-transpose pattern. In this post I use row-major matrixes.

Here is reference scalar code from https://github.com/graphitemaster/normals_revisited

float minor(const float m[16], int r0, int r1, int r2, int c0, int c1, int c2) {
  return m[4*r0+c0] * (m[4*r1+c1] * m[4*r2+c2] - m[4*r2+c1] * m[4*r1+c2]) -
         m[4*r0+c1] * (m[4*r1+c0] * m[4*r2+c2] - m[4*r2+c0] * m[4*r1+c2]) +
         m[4*r0+c2] * (m[4*r1+c0] * m[4*r2+c1] - m[4*r2+c0] * m[4*r1+c1]);
}

void cofactor(const float src[16], float dst[16]) {
  dst[ 0] =  minor(src, 1, 2, 3, 1, 2, 3);
  dst[ 1] = -minor(src, 1, 2, 3, 0, 2, 3);
  dst[ 2] =  minor(src, 1, 2, 3, 0, 1, 3);
  dst[ 3] = -minor(src, 1, 2, 3, 0, 1, 2);
  dst[ 4] = -minor(src, 0, 2, 3, 1, 2, 3);
  dst[ 5] =  minor(src, 0, 2, 3, 0, 2, 3);
  dst[ 6] = -minor(src, 0, 2, 3, 0, 1, 3);
  dst[ 7] =  minor(src, 0, 2, 3, 0, 1, 2);
  dst[ 8] =  minor(src, 0, 1, 3, 1, 2, 3);
  dst[ 9] = -minor(src, 0, 1, 3, 0, 2, 3);
  dst[10] =  minor(src, 0, 1, 3, 0, 1, 3);
  dst[11] = -minor(src, 0, 1, 3, 0, 1, 2);
  dst[12] = -minor(src, 0, 1, 2, 1, 2, 3);
  dst[13] =  minor(src, 0, 1, 2, 0, 2, 3);
  dst[14] = -minor(src, 0, 1, 2, 0, 1, 3);
  dst[15] =  minor(src, 0, 1, 2, 0, 1, 2);
}

With SSE it could be implemented as:

void cofactorSSE(const float src[16], float dst[16]) {
    __m128 r0 = _mm_load_ps(&src[0]);
    __m128 r1 = _mm_load_ps(&src[4]);
    __m128 r2 = _mm_load_ps(&src[8]);
    __m128 r3 = _mm_load_ps(&src[12]);

    __m128 r0_0001 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(0, 0, 0, 1));
    __m128 r1_0001 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(0, 0, 0, 1));
    __m128 r2_0001 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(0, 0, 0, 1));
    __m128 r3_0001 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(0, 0, 0, 1));
    __m128 r0_1122 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(1, 1, 2, 2));
    __m128 r1_1122 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(1, 1, 2, 2));
    __m128 r2_1122 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(1, 1, 2, 2));
    __m128 r3_1122 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(1, 1, 2, 2));
    __m128 r0_2333 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(2, 3, 3, 3));
    __m128 r1_2333 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(2, 3, 3, 3));
    __m128 r2_2333 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(2, 3, 3, 3));
    __m128 r3_2333 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(2, 3, 3, 3));
    __m128 odd = _mm_set_ps(0.0, -0.0, 0.0, -0.0);
    __m128 even = _mm_set_ps(-0.0, 0.0, -0.0, 0.0);

    __m128 res0 = _mm_mul_ps(_mm_mul_ps(r1_0001, r2_1122), r3_2333);
    __m128 res1 = _mm_mul_ps(_mm_mul_ps(r0_0001, r2_1122), r3_2333);
    __m128 res2 = _mm_mul_ps(_mm_mul_ps(r0_0001, r1_1122), r3_2333);
    __m128 res3 = _mm_mul_ps(_mm_mul_ps(r0_0001, r1_1122), r2_2333);
    res0 = _mm_add_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_1122, r2_2333), r3_0001));
    res1 = _mm_add_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_1122, r2_2333), r3_0001));
    res2 = _mm_add_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_2333), r3_0001));
    res3 = _mm_add_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_2333), r2_0001));
    res0 = _mm_add_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_2333, r2_0001), r3_1122));
    res1 = _mm_add_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_2333, r2_0001), r3_1122));
    res2 = _mm_add_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_0001), r3_1122));
    res3 = _mm_add_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_0001), r2_1122));
    res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_2333, r2_1122), r3_0001));
    res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_2333, r2_1122), r3_0001));
    res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_1122), r3_0001));
    res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_1122), r2_0001));
    res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_1122, r2_0001), r3_2333));
    res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_1122, r2_0001), r3_2333));
    res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_0001), r3_2333));
    res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_0001), r2_2333));
    res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_0001, r2_2333), r3_1122));
    res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_0001, r2_2333), r3_1122));
    res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_0001, r1_2333), r3_1122));
    res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_0001, r1_2333), r2_1122));

    _mm_store_ps(&dst[0], _mm_xor_ps(res0, even));
    _mm_store_ps(&dst[4], _mm_xor_ps(res1, odd));
    _mm_store_ps(&dst[8], _mm_xor_ps(res2, even));
    _mm_store_ps(&dst[12], _mm_xor_ps(res3, odd));
}

But it's not that fast compared to scalar code, only about 2.6 times faster (compiled with MSVC 2010, tested on Haswell). Of course, there are many shuffles. Also a couple of spills, which don't help, but there are not as many as it naively looks like there might be (MSVC changes the order and interleaves the shuffles and their uses somewhat).

Can this be done more efficiently? For example, a trick to save some shuffles? SSE through SSE3 can definitely be used, AVX(2) answers are also interesting (but I will have to dual-wield two versions of the function). AVX512 as a bonus but I doubt it will see much use.

\$\endgroup\$
  • \$\begingroup\$ You could use restrict (C11) in your parameters: void cofactorSSE(const float [static restrict 16], float [static restrict 16]). Maybe it would help optimizations. I would be surprised if MSVC supported it, though. \$\endgroup\$ – Cacahuete Frito Aug 6 at 16:40
  • \$\begingroup\$ @CacahueteFrito isn't that a C thing technically? sorry for the C tag by the way someone else added it \$\endgroup\$ – harold Aug 6 at 16:43
  • \$\begingroup\$ Ah sorry then. I saw the C tag :). GCC has it also in C++ as an extension (as __restrict__, but you can easily do the macro). Maybe other compilers also provide it as an extension. \$\endgroup\$ – Cacahuete Frito Aug 6 at 16:43
  • \$\begingroup\$ Also, if you care about performance, why not program the function in C (or Fortran or assembly) and just export it to C++ with an extern "C"? \$\endgroup\$ – Cacahuete Frito Aug 6 at 16:48
  • 1
    \$\begingroup\$ @CacahueteFrito I've tried __restrict now but it didn't do anything. I don't have an assembly version that's actually better, if you have one I can consider that \$\endgroup\$ – harold Aug 6 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.