3D matrix rotation in homogeneous coordinate space

Question

Assuming that the framework is in place to handle the difference between row and column major matrices, I am curious to know if in a header based library such implementation is semantically and syntactically correct and or appealing. I am not using PIMPL here because it is strictly a header library.

template <typename T, matrix_layout ML>                                                                                                                                                                                                                                                                                    
matrix_4X4<T, ML> rotate_by(const matrix_4X4<T, ML> &m, const T &angle,                                                                                                                                                                                                                                                    
  const vector_3d<T> &v, angle_mode mode = radians)                                                                                                                                                                                                                                                                        
{                                                                                                                                                                                                                                                                                                                          
  T a = angle;                                                                                                                                                                                                                                                                                                             
  if (mode == degrees) {                                                                                                                                                                                                                                                                                                   
    a = degrees_to_radians(angle);                                                                                                                                                                                                                                                                                         
  }                                                                                                                                                                                                                                                                                                                        
  const T c = cos<T, radians>(a);                                                                                                                                                                                                                                                                                          
  const T s = sin<T, radians>(a);                                                                                                                                                                                                                                                                                          
  const T k = (T(1) - c);                                                                                                                                                                                                                                                                                                  

  vector_3d<T> axis(v);                                                                                                                                                                                                                                                                                                    
  axis.normalize();                                                                                                                                                                                                                                                                                                        
  const T &x = axis[0];                                                                                                                                                                                                                                                                                                    
  const T &y = axis[1];                                                                                                                                                                                                                                                                                                    
  const T &z = axis[2];                                                                                                                                                                                                                                                                                                    

  // The matrix to rotate m by.                                                                                                                                                                                                                                                                                            
  matrix_4X4<T, ML> n(null);                                                                                                                                                                                                                                                                                               
  matrix_4X4<T, ML> b(m);                                                                                                                                                                                                                                                                                                  

  if (ML == column) {                                                                                                                                                                                                                                                                                                      
    /* [col][row] */                                                                                                                                                                                                                                                                                                       
    n[0][0] = (x * x * k) + (c);                                                                                                                                                                                                                                                                                           
    n[0][1] = (y * x * k) + (z * s);                                                                                                                                                                                                                                                                                       
    n[0][2] = (x * z * k) - (y * s);                                                                                                                                                                                                                                                                                       

    n[1][0] = (x * y * k) - (z * s);                                                                                                                                                                                                                                                                                       
    n[1][1] = (y * y * k) + (c);                                                                                                                                                                                                                                                                                           
    n[1][2] = (y * z * k) + (x * s);                                                                                                                                                                                                                                                                                       

    n[2][0] = (x * z * k) + (y * s);                                                                                                                                                                                                                                                                                       
    n[2][1] = (y * z * k) - (x * s);                                                                                                                                                                                                                                                                                       
    n[2][2] = (z * z * k) + (c);                                                                                                                                                                                                                                                                                           

    b[0] = m[0] * n[0][0] + m[1] * n[0][1] + m[2] * n[0][2];                                                                                                                                                                                                                                                               
    b[1] = m[0] * n[1][0] + m[1] * n[1][1] + m[2] * n[1][2];                                                                                                                                                                                                                                                               
    b[2] = m[0] * n[2][0] + m[1] * n[2][1] + m[2] * n[2][2];                                                                                                                                                                                                                                                               
    b[3] = m[3];                                                                                                                                                                                                                                                                                                           
  } else if (ML == row) {                                                                                                                                                                                                                                                                                                  
    /* [row][col] */                                                                                                                                                                                                                                                                                                       
    n[0][0] = (x * x * k) + (c);                                                                                                                                                                                                                                                                                           
    n[0][1] = (x * y * k) - (z * s);                                                                                                                                                                                                                                                                                       
    n[0][2] = (x * z * k) + (y * s);                                                                                                                                                                                                                                                                                       

    n[1][0] = (y * x * k) + (z * s);                                                                                                                                                                                                                                                                                       
    n[1][1] = (y * y * k) + (c);                                                                                                                                                                                                                                                                                           
    n[1][2] = (y * z * k) - (x * s);                                                                                                                                                                                                                                                                                       

    n[2][0] = (x * z * k) - (y * s);                                                                                                                                                                                                                                                                                       
    n[2][1] = (y * z * k) + (x * s);                                                                                                                                                                                                                                                                                       
    n[2][2] = (z * z * k) + (c);                                                                                                                                                                                                                                                                                           

    b[0] = m[0] * n[0][0] + m[1] * n[1][0] + m[2] * n[2][0];                                                                                                                                                                                                                                                               
    b[1] = m[0] * n[0][1] + m[1] * n[1][1] + m[2] * n[2][1];                                                                                                                                                                                                                                                               
    b[2] = m[0] * n[0][2] + m[1] * n[1][2] + m[2] * n[2][2];                                                                                                                                                                                                                                                               
    b[3] = m[3];                                                                                                                                                                                                                                                                                                           
  }                                                                                                                                                                                                                                                                                                                        

  return b;                                                                                                                                                                                                                                                                                                                
}                                                                                                                                                                                                                                                                                           

Answer 1

You have a matrix class. Presumably, the class is part of a library that supports basic matrix operations such as multiplication and addition.

It's odd, therefore, that a rotate_by() function would implement matrix operations from scratch. What you want to express is \$\mathbf{b} = \mathrm{N}^T \mathbf{m}\$, and it should be written that way. All of the details of how to perform that multiplication, including the row-major vs. column-major layout, should be taken care of elsewhere.

The definitions of the elements of \$\mathrm{N}\$ look suspicious to me, as the units of the addends don't agree.