#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// Class template for computing the barycentric coordinates of a point
// relative to a given triangle. Typically, this class will be
// instantiated for either <float> or <double> type, depending
// on the desired precision.
//
// This code precomputes and stores data that speed up computing
// the barycentric coordinates. This is particularly useful when
// repeatedly evaluating different points on the same triangle.
//
// (C) 2002 Ivan Neulander

template <class T> class Barycentric {
  int _proj; // which dimension we discard

  T _v0[3];  // original position of point that we
             // moved to the origin

  T _v1[2];  // the other two points with one dimension discarded
  T _v2[2]; 

  T _bv1[2]; // precomputed precursor to barycentric
  T _bv2[2];

public:
  Barycentric() { Init(); }

  Barycentric(const T *v0, const T *v1, const T *v2) 
  { Init(v0, v1, v2); }

  static void Diff(const T *a, const T *b, T *res) {
    res[0] = a[0] - b[0];
    res[1] = a[1] - b[1];
    res[2] = a[2] - b[2];
  }

  void Init() {
    _proj = -1;
    _v0[0] = _v0[1] = _v0[2] = 0;
    _v1[0] =_v1[1] = 0;
    _v2[0] =_v2[1] = 0;
    _bv1[0] =_bv1[1] = 0;
    _bv2[0] =_bv2[1] = 0;
  }
  
  void Init(const T *iv0, const T *iv1, const T *iv2);
  // Call this to specify the 3 verts of the triangle for which
  // barycentric coordinates are to be computed.


  bool IsValid() const { return _proj >= 0; }
  // True iff a non-degenerate triangle has been specified
  // using constructor or Init().
  
  bool Eval(const T *pos, T *bary) {
    // Attempts to compute barycentric coords of point 'pos' relative
    // to the triangle specified using constructor or Init().  This
    // method works even if the point lies outside triangle; however
    // the point must always be coplanar with the triangle. Returns
    // true on success, false on failure (i.e. triangle is
    // degenerate or uninitialized).

    T p0, p1;
    switch(_proj) {
    case 0:
      p0 = pos[1] - _v0[1];
      p1 = pos[2] - _v0[2];
      break;
    case 1:
      p0 = pos[2] - _v0[2];
      p1 = pos[0] - _v0[0];
      break;
    case 2:
      p0 = pos[0] - _v0[0];
      p1 = pos[1] - _v0[1];
      break;
    default:
      return false;
    }

    bary[1] = _bv2[1]*p0 - _bv2[0]*p1;
    bary[2] = _bv1[0]*p1 - _bv1[1]*p0;
    bary[0] = 1 - bary[1] - bary[2];

    return true;
  }

  bool Test(const T *pos, T *bary, T min = 0, T max = 1) {
    // This is an optimized version of Eval(). It only computes
    // barycentric coordinates if point 'pos' lies within the
    // triangle. If so, it returns true. Otherwise it returns
    // false. As with Eval(), it assumes that 'pos' is coplanar with
    // triangle.

    T p0, p1;

    switch(_proj) {
    case 0:
      p0 = pos[1] - _v0[1];
      p1 = pos[2] - _v0[2];
      break;
    case 1:
      p0 = pos[2] - _v0[2];
      p1 = pos[0] - _v0[0];
      break;
    case 2:
      p0 = pos[0] - _v0[0];
      p1 = pos[1] - _v0[1];
      break;
    default:
      return false;
    }

    bary[1] = _bv2[1]*p0 - _bv2[0]*p1;
    if (bary[1] < min || bary[1] > max) return false;

    bary[2] = _bv1[0]*p1 - _bv1[1]*p0;
    if (bary[2] < min || bary[2] > max) return false;

    bary[0] = 1 - bary[1] - bary[2];
    return bary[0] >= min;
  }
};

template <class T> void Barycentric<T>::
 Init(const T *iv0, const T *iv1, const T *iv2) { 
  // Invoke this method before evaluating or testing any barycentric
  // coords. It reads in the verts (iv0, iv1, iv2) of the triangle to
  // be used as a basis for barycentric coordinates, and precomputes
  // needed information. Returns 'true' if triangle is valid, false
  // otherwise.
  
  T v1[3], v2[3];
  
  Init();
  
  // copy iv0
  _v0[0] = iv0[0];
  _v0[1] = iv0[1];
  _v0[2] = iv0[2];
  
  // translate everything to place iv0 at origin
  Diff(iv1, iv0, v1);
  Diff(iv2, iv0, v2);
  
  // The system we're solving is overdetermined; We solve using one of
  // the following projections: (y,z), (z,x), (x,y): namely, the one
  // that yields the biggest determinant. So we just need to decide
  // which dimension to discard.
  
  T det0 = v1[1]*v2[2] - v2[1]*v1[2];
  T det1 = v1[2]*v2[0] - v2[2]*v1[0];
  T det2 = v1[0]*v2[1] - v2[0]*v1[1];
  
  T adet0 = fabs(det0);
  T adet1 = fabs(det1);
  T adet2 = fabs(det2);
  
  // Pick proj according to the biggest among adet0, adet1, adet2.
  // Note that if they're all zero, then proj will be assigned 2, so
  // that's the only case where we need to check before dividing by
  // adet2.
  int proj = adet0 > adet1 ? 
    (adet0 > adet2 ? 0 : 2) : (adet1 > adet2 ? 1 : 2);
  
  T detInv;
  switch(proj) {
  case 0:
    detInv = 1/det0;
    _proj = 0;
    
    _bv1[0] = detInv*v1[1];
    _bv1[1] = detInv*v1[2];
    _bv2[0] = detInv*v2[1];
    _bv2[1] = detInv*v2[2];
    
    _v1[0] = v1[1];
    _v1[1] = v1[2];
    _v2[0] = v2[1];
    _v2[1] = v2[2];
    
    break;
    
  case 1:
    detInv = 1/det1;
    _proj = 1;
    
    _bv1[0] = detInv*v1[2];
    _bv1[1] = detInv*v1[0];
    _bv2[0] = detInv*v2[2];
    _bv2[1] = detInv*v2[0];

    _v1[0] = v1[2];
    _v1[1] = v1[0];
    _v2[0] = v2[2];
    _v2[1] = v2[0];
    
    break;
    
  case 2:
    if (adet2 == 0) return;
    
    detInv = 1/det2;
    _proj = 2;
    
    _bv1[0] = detInv*v1[0];
    _bv1[1] = detInv*v1[1];
    _bv2[0] = detInv*v2[0];
    _bv2[1] = detInv*v2[1];
    
    _v1[0] = v1[0];
    _v1[1] = v1[1];
    _v2[0] = v2[0];
    _v2[1] = v2[1];
    
    break;
  }
}

/////////////////////////////////////////////////////////////////////

// This code exercises the above Barycentric class. The user must
// first enter three 3-vectors to define the vertices of the triangle.
// Each subsequent 3-vector input is interpreted as a point whose
// barycentric coordinates we wish to compute relative to that
// triangle.

bool syntaxError() {
  fprintf(stderr,"Invalid number!\n");
  return false;
}

bool isaNum(char s) {
  return (s >= '0' && s <= '9') || s=='-' || s=='.';
}

bool parseVector(const char *line, double *vec) {
  const char *s = line;
  int i;
  for(i = 0 ; i < 3 ; i++) {
    while (*s && !isaNum(*s)) s++;
    if (!*s) return syntaxError();
    vec[i] = atof(s);
    while (isaNum(*s)) s++;
  }
  while (*s && !isaNum(*s)) s++;
  if (*s) return syntaxError();

  return true;
}

int main() {
  const int LINE_LENGTH = 255;
  char line[LINE_LENGTH+1];
  double tri[3][3];

  printf("[Type Ctrl-D at any time to exit.]\n\n");

  for (int i = 0 ; i < 3 ; i++) {
    do {
      printf("Enter triangle v%d: ", i);
      if (!fgets(line, LINE_LENGTH, stdin)) return 0;
      line[LINE_LENGTH] = 0;
    } while (!parseVector(line, tri[i]));
  }
  
  Barycentric<double> bary(tri[0], tri[1], tri[2]);
  if (!bary.IsValid()) {
    printf("\nDegenerate triangle!\n");
    return 1;
  }

  while(1) {
    double pos[3];
    do {
      printf("Enter point: ");
      if (!fgets(line, LINE_LENGTH, stdin)) return 0;
      line[LINE_LENGTH] = 0;
    } while(!parseVector(line, pos));

    double b[3];
    if (bary.Eval(pos, b)) 
      printf("  bary ==> (%g %g %g)\n", b[0], b[1], b[2]);
    else
      printf("  Could not compute bary!\n");
  }
}