/*
 *  Ray++ - Object-oriented ray tracing library
 *  Copyright (C) 1998-2001 Martin Reinecke and others.
 *  See the AUTHORS file for more information.
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Library General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Library General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public
 *  License along with this library; if not, write to the Free
 *  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  See the README file for more information.
 */

#include "shapes/quadric.h"

namespace RAYPP {

const float8 QUADRIC::Epsilon = Small_float8*Small_float8;

VECTOR QUADRIC::Get_Normal (const VECTOR &loc) const
  {
  VECTOR tmpvec (0,0,0);

  if (Square_Terms)
    tmpvec += VECTOR (2*xx*loc.x, 2*yy*loc.y, 2*zz*loc.z);
  if (Mixed_Terms)
    tmpvec += VECTOR (xy*loc.y + xz*loc.z,
                      xy*loc.x + yz*loc.z,
                      xz*loc.x + yz*loc.y);
  if (Linear_Terms)
    tmpvec += VECTOR (x,y,z); 

  if (tmpvec.SquaredLength() < Small_dist*Small_dist) tmpvec = VECTOR(0,1,0);
  tmpvec.Normalize();
  if (Inverted) tmpvec.Flip();
  return tmpvec;
  }

void QUADRIC::Get_Coeffs (const GEOM_RAY &Ray, float8 &square_coeff,
                          float8 &linear_coeff, float8 &const_coeff) const
  {
  square_coeff = linear_coeff = const_coeff = 0;

  if (Square_Terms)
    {
    square_coeff += xx*Ray.dir.x*Ray.dir.x; 
    square_coeff += yy*Ray.dir.y*Ray.dir.y; 
    square_coeff += zz*Ray.dir.z*Ray.dir.z; 

    linear_coeff += 2*xx*Ray.start.x*Ray.dir.x;
    linear_coeff += 2*yy*Ray.start.y*Ray.dir.y;
    linear_coeff += 2*zz*Ray.start.z*Ray.dir.z;

    const_coeff  += xx*Ray.start.x*Ray.start.x;
    const_coeff  += yy*Ray.start.y*Ray.start.y;
    const_coeff  += zz*Ray.start.z*Ray.start.z;
    }

  if (Mixed_Terms)
    {
    square_coeff +=  xy*Ray.dir.x*Ray.dir.y; 
    square_coeff +=  xz*Ray.dir.x*Ray.dir.z; 
    square_coeff +=  yz*Ray.dir.y*Ray.dir.z; 

    linear_coeff += xy * (Ray.start.x*Ray.dir.y + Ray.start.y*Ray.dir.x);
    linear_coeff += xz * (Ray.start.x*Ray.dir.z + Ray.start.z*Ray.dir.x);
    linear_coeff += yz * (Ray.start.y*Ray.dir.z + Ray.start.z*Ray.dir.y);

    const_coeff  +=  xy*Ray.start.x*Ray.start.y;
    const_coeff  +=  xz*Ray.start.x*Ray.start.z;
    const_coeff  +=  yz*Ray.start.y*Ray.start.z;
    }

  if (Linear_Terms)
    {
    linear_coeff +=
      x*Ray.dir.x + y*Ray.dir.y + z*Ray.dir.z;
    const_coeff  +=
      x*Ray.start.x+y*Ray.start.y+z*Ray.start.z;
    }

  const_coeff  += C;
  }

QUADRIC::QUADRIC ()
  {
  xx = yy = zz = 1;
  xy = xz = yz = x = y = z = 0;
  C = -1;
  }

QUADRIC::QUADRIC (float4 pxx, float4 pyy, float4 pzz,
                  float4 pxy, float4 pxz, float4 pyz,
                  float4 px,  float4 py,  float4 pz,
                  float4 pC)
  {
  xx = pxx; yy = pyy; zz = pzz;
  xy = pxy; xz = pxz; yz = pyz;
  x  = px;  y  = py;  z  = pz;
  C  = pC;
  }

//virtual
void QUADRIC::Init ()
  {
  if (initialized) return;
  Square_Terms = ((xx*xx + yy*yy + zz*zz) > Epsilon);
  Mixed_Terms = ((xy*xy + xz*xz + yz*yz) > Epsilon);
  Linear_Terms = ((x*x + y*y + z*z) > Epsilon);
  initialized = true;
  }

//virtual
void QUADRIC::Transform (const TRANSFORM &trans)
  {
  const TRANSMAT &t = trans.Inverse();
  float4 tmpmat[4][4], tmpmat2[4][4], tmpmat3[4][4];
  uint2 m, n, i;

  tmpmat[0][0] = xx;
  tmpmat[1][1] = yy;
  tmpmat[2][2] = zz;

  tmpmat[0][1] = xy;
  tmpmat[0][2] = xz;
  tmpmat[1][2] = yz;

  tmpmat[1][0] = tmpmat[2][0] = tmpmat[2][1] = 0;

  tmpmat[3][0] = x;
  tmpmat[3][1] = y;
  tmpmat[3][2] = z;

  tmpmat[0][3] = tmpmat[1][3] = tmpmat[2][3] = 0;

  tmpmat[3][3] = C;

  for (m=0; m<4; m++)
    for (n=0; n<3; n++)
      tmpmat2 [m][n] = t.entry[m][n];

  tmpmat2[0][3] = tmpmat2[1][3] = tmpmat2[2][3] = 0;
  tmpmat2[3][3] = 1;

  for (m=0; m<4; m++)
    for (n=0; n<4; n++)
      {
      tmpmat3[m][n] = 0;
      for (i=0; i<4; i++)
        tmpmat3[m][n] += tmpmat[m][i]*tmpmat2[n][i];
      }

  for (m=0; m<4; m++)
    for (n=0; n<4; n++)
      {
      tmpmat[m][n] = 0;
      for (i=0; i<4; i++)
        tmpmat[m][n] += tmpmat2[m][i]*tmpmat3[i][n];
      }

  xx = tmpmat[0][0];   
  yy = tmpmat[1][1];   
  zz = tmpmat[2][2];

  xy = tmpmat[0][1] + tmpmat[1][0];
  xz = tmpmat[0][2] + tmpmat[2][0];
  yz = tmpmat[1][2] + tmpmat[2][1];

  x = tmpmat[3][0] + tmpmat[0][3];
  y = tmpmat[3][1] + tmpmat[1][3];
  z = tmpmat[3][2] + tmpmat[2][3];

  C = tmpmat[3][3];
  }

//virtual
AXISBOX QUADRIC::BBox () const
  {
// this could be improved
  return AXISBOX (VECTOR (-Huge_float8, -Huge_float8, -Huge_float8),
                  VECTOR ( Huge_float8,  Huge_float8,  Huge_float8));
  }

//virtual
bool QUADRIC::Test (const GEOM_RAY &Ray, float8 &dist, bool &realhit) const
  {
  float8 square_coeff, linear_coeff, const_coeff;

  Get_Coeffs (Ray, square_coeff, linear_coeff, const_coeff);

  if (abs(square_coeff) < Epsilon)
    {
    if (abs(linear_coeff) < Epsilon) return false;
    dist = -const_coeff/linear_coeff;
    if ((dist < Ray.mindist) || (dist > Ray.maxdist)) return false;
    }
  else
    {
    float8 discriminant, tmp1;
    discriminant = linear_coeff*linear_coeff - 4*square_coeff*const_coeff;
    if (discriminant < Epsilon) return false;

    tmp1 = sqrt (discriminant);
    if (square_coeff < 0) tmp1=-tmp1;
    dist = (-linear_coeff - tmp1)/(2*square_coeff);

    if (dist > Ray.maxdist) return false;

    if (dist < Ray.mindist)
      {
      dist = (-linear_coeff + tmp1)/(2*square_coeff);
      if ((dist < Ray.mindist) || (dist > Ray.maxdist)) return false;
      }
    }

  realhit = true;
  return true;
  }

//virtual
bool QUADRIC::Inside (const VECTOR &loc) const
  {
  float8 dens;

  dens = C;
  if (Square_Terms)
    dens += xx*loc.x*loc.x + yy*loc.y*loc.y + zz*loc.z*loc.z;
  if (Mixed_Terms)
    dens += xy*loc.x*loc.y + xz*loc.x*loc.z + xz*loc.y*loc.z;
  if (Linear_Terms)
    dens += x*loc.x        + y*loc.y        + z*loc.z;

  if (dens<0) return (!Inverted);
  return Inverted; 
  }

//virtual
bool QUADRIC::Intersect
  (const GEOM_RAY &Ray, float8 &dist, VECTOR &Normal) const
  {
  bool dummy;
  if (!Test (Ray, dist, dummy)) return false;

  Normal = Get_Normal (Ray.eval(dist));
  return true;
  }

//virtual
void QUADRIC::All_Intersections
  (const GEOM_RAY &Ray, vector<INTER> &Inter) const
  {
  float8 dist;
  float8 square_coeff, linear_coeff, const_coeff;

  Get_Coeffs (Ray, square_coeff, linear_coeff, const_coeff);

  if (abs(square_coeff) < Epsilon)
    {
    if (abs(linear_coeff) < Epsilon) return;
    dist = -const_coeff/linear_coeff;
    if ((dist < Ray.mindist) || (dist > Ray.maxdist)) return;
    Inter.push_back (INTER (dist, Get_Normal (Ray.eval(dist))));
    }
  else
    {
    float8 discriminant;
    discriminant = linear_coeff*linear_coeff - 4*square_coeff*const_coeff;
    if (discriminant < Epsilon) return;

    dist = (-linear_coeff - sqrt(discriminant))/(2*square_coeff);
    if ((dist>Ray.mindist) && (dist<Ray.maxdist))
      {
      Inter.push_back (INTER (dist, Get_Normal (Ray.eval(dist))));
      }
    dist = (-linear_coeff + sqrt(discriminant))/(2*square_coeff);
    if ((dist>Ray.mindist) && (dist<Ray.maxdist))
      {
      Inter.push_back (INTER (dist, Get_Normal (Ray.eval(dist))));
      }
    }
  }

void QUADRIC::Set_Params (float4 pxx, float4 pyy, float4 pzz,
                          float4 pxy, float4 pxz, float4 pyz,
                          float4 px,  float4 py,  float4 pz,
                          float4 pC)
  {
  xx = pxx; yy = pyy; zz = pzz;
  xy = pxy; xz = pxz; yz = pyz;
  x  = px;  y  = py;  z  = pz;
  C  = pC;
  }

} // namespace RAYPP