Quadric3.cpp

#include "Quadric3.hh"

using namespace BIAS;
using namespace std;

namespace BIAS {

Quadric3::Quadric3(double a, double b, double c, double d, double area){
  a2 = a*a;
  ab = a*b;
  ac = a*c;
  ad = a*d;
  b2 = b*b;
  bc = b*c;
  bd = b*d;
  c2 = c*c;
  cd = c*d;
  d2 = d*d;

  r = area;
}

Vector3<double> Quadric3::GetVector(){
  return Vector3<double>(ad,bd,cd);
}

Matrix3x3<double> Quadric3::GetTensor(){
  return Matrix3x3<double>(a2,ab,ac,
                   ab,b2,bc,
                   ac,bc,c2);
}

Matrix4x4<double> Quadric3::GetHomogenous(){
  return Matrix4x4<double>(a2,ab,ac,ad,
                   ab,b2,bc,bd,
                   ac,bc,c2,cd,
                   ad,bd,cd,d2);
}

double Quadric3::Evaluate(double x, double y, double z){
  return x*x*a2 + 2*x*y*ab + 2*x*z*ac + 2*x*ad
                + y*y*b2   + 2*y*z*bc + 2*y*bd
                           + z*z*c2   + 2*z*cd
                                      + d2;
}

bool Quadric3::Optimize(Vector3<double>& v){
  Matrix3x3<double> Ainv;
  Ainv = GetTensor();
  Ainv.InvertIP();
  double det = Ainv.GetDeterminant();

  if(det < 1e-12)
    return false;

  Ainv.Mult(GetVector(),v);
  v.MultIP(-1);
  return true;
}

Quadric3& Quadric3::
operator=(const Quadric3& Q){
    r = Q.r;

    a2 = Q.a2;  ab = Q.ab;  ac = Q.ac;  ad = Q.ad;
    b2 = Q.b2;  bc = Q.bc;  bd = Q.bd;
    c2 = Q.c2;  cd = Q.cd;
    d2 = Q.d2;

    return *this;
}

Quadric3& Quadric3::
operator+=(const Quadric3& Q)
{
    // Accumulate area
    r += Q.r;

    // Accumulate coefficients
    a2 += Q.a2;  ab += Q.ab;  ac += Q.ac;  ad += Q.ad;
                 b2 += Q.b2;  bc += Q.bc;  bd += Q.bd;
                              c2 += Q.c2;  cd += Q.cd;
                                           d2 += Q.d2;

    return *this;
}

Quadric3& Quadric3::
operator*=(double s)
{
    // Accumulate coefficients
    a2 *= s;  ab *= s;  ac *= s;  ad *= s;
              b2 *= s;  bc *= s;  bd *= s;
                        c2 *= s;  cd *= s;
                                  d2 *= s;

    return *this;
}

}//end namespace BIAS