/* Output:
     q is p reorded so that the component with the largest magnitude
     is in the z-coordinate. In addition dqdx and dqdy are the reorderings
     of dpdx and dpdy to mirror the reordering in q.
 */
void
compute_selector_fast_and_brave(in vec3 p, in vec3 dpdx, in vec3 dpdy,
				out vec3 q, out vec3 dqdx, out vec3 dqdy)
{
  vec3 absP;
  float m;

  absP = abs(p);
  m = max(absP.x, max(absP.y, absP.z));
  if (m == absP.x)
    {
      q = p.yzx;
      dqdx = dpdx.yzx;
      dqdy = dpdy.yzx;
    } 
  else if (m == absP.y)
    {
      q = p.xzy;
      dqdx = dpdx.xzy;
      dqdy = dpdy.xzy;
    }
  else
    {
      q = p;
      dqdx = dpdx;
      dqdy = dpdy;
    }
}

void
compute_selector_slow_and_careful(in vec3 p, in vec3 dpdx, in vec3 dpdy,
				  out vec3 q, out vec3 dqdx, out vec3 dqdy)
{
  vec3 absP;
  float m;

  absP = abs(p);
  if (absP.x >= max(absP.y, absP.z))
    {
      q = p.yzx;
      dqdx = dpdx.yzx;
      dqdy = dpdy.yzx;
    } 
  else if (absP.y >= max(absP.x, absP.z))
    {
      q = p.xzy;
      dqdx = dpdx.xzy;
      dqdy = dpdy.xzy;
    }
  else
    {
      q = p;
      dqdx = dpdx;
      dqdy = dpdy;
    }
}


/* Choose the compute_selector method
 */
#ifdef FAST_AND_BRAVE
#define compute_selector compute_selector_fast_and_brave
#else
#define compute_selector compute_selector_slow_and_careful
#endif

/* Perform the quoteint rule to compute the derivative
   of q.xy / q.z
 */
void
quotient_rule(in vec3 q, in vec3 dqdx, in vec3 dqdy,
	      out vec2 dx, out vec2 dy) 
	      
{
  float recip;

  recip = 1.0 / q.z;

  dx = recip * ( dqdx.xy - q.xy * (dqdx.z * recip) );
  dy = recip * ( dqdy.xy - q.xy * (dqdy.z * recip) );

}

/* given position and derivatives of position
   compute the LOD for a cubemap texture look
   up where the dimension of the cubemap is L
 */
float
compute_textureCubeLOD(in vec3 p, in vec3 dpdx, in vec3 dpdy,
		       in float L)
{
  vec3 q, dqdx, dqdy;
  vec2 dx, dy;
  float m;

  compute_selector(p, dpdx, dpdy, q, dqdx, dqdy);
  quotient_rule(q, dqdx, dqdy, dx, dy);

  m = max(dot(dx, dx), dot(dy, dy));

  /*
     rho = 0.5 * L * max( ||dx||, ||dy||)
         = 0.5 * L * sqrt( max(dot(dx, dx), dot(dy, dy)) )
	 = 0.5 * sqrt(L * L * m)

     LOD = log2(rho)
         = log2(0.5 * sqrt(L * L * m) ); 
	 = -1.0 + 0.5 * log2(L * L  * m);
  
   */
  return -1.0 + 0.5 * log2(L * L * m);
  
}