planeCoefficients

//Three doubles (x, y, z)
struct D3{
   double x;
   double y;
   double z;

   D3(double x_ = 0, double y_ = 0, double z_ = 0):
   x(x_), y(y_), z(z_){}

   D3 operator-(const D3 &rhs) const{
      return D3(x - rhs.x,
                y - rhs.y,
                z - rhs.z);
   }

   D3 operator-() const{
      return D3(-x, -y, -z);
   }

   D3 operator/=(double s) const{
      return D3(x/s, y/s, z/s);
   }
};

//Four doubles (a, b, c, d)
struct D4{
   double a;
   double b;
   double c;
   double d;

   D4(double a_ = 0, double b_ = 0, double c_ = 0, double d_ = 0):
   a(a_), b(b_), c(c_), d(d_){}
};

double determinantFromVectors(const D3 &x, const D3 &y, const D3 &z){
   return
      + x.x * y.y * z.z
      + y.x * z.y * x.z
      + z.x * x.y * y.z

      - z.x * y.y * x.z
      - y.x * x.y * z.z
      - x.x * z.y * y.z;
}

D3 crossProduct(const D3 &u, const D3 &v){
   return D3(u.y * v.z - u.z * v.y,
             u.z * v.x - u.x * v.z,
             u.x * v.y - u.y * v.z);
}

void normalize(D3 &v){
   double length = sqrt(v.x * v.x +
                        v.y * v.y +
                        v.z * v.z);
   v /= length;
}

D4 planeCoefficients(const D3 &p1, const D3 &p2, const D3 &p3){
   //Option 1: Matrix determinants
   {
      // v1 vX vY vZ
      // -----------
      //  1 x1 y1 z1
      //  1 x2 y2 z2
      //  1 x3 y3 z3
      D3
         v1(1, 1, 1),
         vX(p1.x, p2.x, p3.x),
         vY(p1.y, p2.y, p3.y),
         vZ(p1.z, p2.z, p3.z);

      return D4(determinantFromVectors(v1, vY, vZ),
                determinantFromVectors(vX, v1, vZ),
                determinantFromVectors(vX, v1, vZ),
               -determinantFromVectors(vX, vY, vZ));
   }

   //Option 2: Cross-product of vectors
   {
      //   |
      //   |
      //   |
      // p1|__________p2
      //    \    v2
      //     \v1
      //      \
      //       \p3
      D3
         v1 = p3 - p1,
         v2 = p2 - p1,
         cross = crossProduct(v1, v2);
      normalize(cross);

      return D4(cross.x,
                cross.y,
                cross.z,

                cross.x * p3.x +
                cross.y * p3.y +
                cross.z * p3.z);
   }
}