/* 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

void
compute_vec2s(in vec3 q, in vec3 dqdx, in vec3 dqdy,
	      out vec2 uv, out vec2 duvdx, out vec2 duvdy) 
	      
{
  float recip, common;

  recip = 1.0 / q.z;
  common = recip * recip

  uv.xy = q.xy * recip;
  duvdx = dqdx * recip - q.xy * (dqdx.z * common);
  duvdy = dqdy * recip - q.xy * (dqdy.z * common);
}

/* 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 uv, duvdx, duvdy;
  float m;

  compute_selector(p, dpdx, dpdy, q, dqdx, dqdy);
  compute_vec2s(q, dqdx, dqdy, uv, duvdx, duvdy);

  m = max(dot(duvdx, duvdx), dot(duvdy, duvdy));

  /*
     rho = 0.5 * L * max( ||duvdx||, ||duvdy||)
         = 0.5 * L * sqrt( max(dot(duvdx, duvdx), dot(duvdy, duvdy)) )
	 = 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(m * L * L);
  
}