#include <d4minboundary.h>

#include <iostream>
#include <algorithm>
using namespace std;

#include <mathlib.h>
#include <print.h>


#include <d4tess.h>

/*

d4mingreedy::d4mingreedy(d4tess & _tess)
  : d4minoperator(_tess) 
{
}

bool const d4mingreedy::eval(uintc a, uintc b)
{
  return tess.tet2to3_(a,b);
}
*/





d4minboundary::d4minboundary(d4tess & _tess)
  : d4minoperator(_tess)
{
}






bool const d4minboundary::eval(uintc a, uintc b)
{
MessageGlobal mg;

  mg() << "a=" << a << " b=" << b << " convex " << tess.isconvex(a,b) << endl;

  if (tess.isconvex(a,b)==false)
    return false;

  uint aj[5];

  d4tri const A(tess.vi[a]);
  aj[0] = A.niInverse(b);

  d4tri const B(tess.vi[b]);
  aj[4] = B.niInverse(a);

  A.getclockwiseface(aj[1],aj[2],aj[3],aj[0]);

  uint pj[5];
  for (uint i=0; i<4; ++i)
    pj[i] = A.pi[aj[i]];
  pj[4] = B.pi[aj[4]];

for (uint i=0; i<5; ++i)
  mg() << "aj[" << i << "]=" << aj[i] << " ";
mg() << endl;
for (uint i=0; i<5; ++i)
  mg() << "pj[" << i << "]=" << pj[i] << " ";
mg() << endl;

  // Find which combination of tetrahedrons has 
  // the minimal boundary.

  // There are two possible tetrahedron permuations.
  // The current on has one inner partition through p0,p1,p2.
  //
  // The other permutation has three blades off line p0 to p4, to points p1,p2,p3.
  //
  // Choose the permutation with the minimum boundary area.

  double area1;
  trianglearea
  (
    area1,
    tess.pt[pj[1]],
    tess.pt[pj[2]],
    tess.pt[pj[3]]
  );

mg() << SHOW(area1) << endl;

  double area2=0.0;

  double area;

  for (uint i=1; i<4; ++i)
  {
    trianglearea
    (
      area,
      tess.pt[pj[4]],
      tess.pt[pj[i]],
      tess.pt[pj[0]]
    );

    area2 += area;
  }

mg() << SHOW(area2) << endl;

  if (area1<area2)
    return false;

mg() << "FUCK" << endl;


  // nj[0] is ignored.
  // the indexes correspond with the points : the neighbour is opposite the point
  // with the same index.  For neighbours >3 simply subtract 3 and that is the
  // index for the second tetrahedron.
  uint nj[7];
  
  for (uint i=1; i<=3; ++i)
    nj[i] = A.ni[aj[i]];

  for (uint i=4; i<=6; ++i)
    nj[i] = B.ni[ B.piInverse(pj[i-3]) ];

for (uint i=1; i<=6; ++i)
  mg() << "nj[" << i << "]=" << nj[i] << " ";
mg() << endl;

  uint k1(a), k2(b), k3(tess.vi.size());

  tess.vi[a] = d4tri( pj[4],pj[2],pj[0],pj[3], nj[1],k3,nj[4],k2 );
  tess.vi[b] = d4tri( pj[4],pj[0],pj[2],pj[1], nj[3],nj[6],k3,k1 );
  tess.vi.push_back( d4tri( pj[4],pj[0],pj[1],pj[3], nj[2],nj[5],k1,k2 ) );

  if (nj[1]!=0)
    tess.vi[nj[1]].niReset(a,k1);
  if (nj[2]!=0)
    tess.vi[nj[2]].niReset(a,k3);
  if (nj[3]!=0)
    tess.vi[nj[3]].niReset(a,k2);
  if (nj[4]!=0)
    tess.vi[nj[4]].niReset(b,k1);
  if (nj[5]!=0)
    tess.vi[nj[5]].niReset(b,k3);
  if (nj[6]!=0)
    tess.vi[nj[6]].niReset(b,k2);

  tess.debugcheck();

  return true;
}


bool const d4mingreedy2::eval(uintc a, uintc b)
{
  if (tess.isconvex(a,b)==false)
    return false;

  d4tri A(tess.vi[a]);
  d4tri B(tess.vi[b]);

  uint ai = A.niInverse(b);
  uint bi = B.niInverse(a);

  double const z = dist(A.pi[ai], B.pi[bi]);

  double l[3];
  uint u[3];
  A.getanticlockwiseface(u[0],u[1],u[2],ai);
  l[0] = dist(u[0],u[1]);
  l[1] = dist(u[1],u[2]);
  l[2] = dist(u[2],u[0]);



  // Order l[0],l[1],l[2] from smallest to largest.

  sort(l,l+3);

cout << SHOW(l[0]) << " ";
cout << SHOW(l[1]) << " ";
cout << SHOW(l[2]) << " ";
cout << SHOW(z) << endl;

  assert(l[2]!=0.0);
  
  // This is the ideal length of z when the tetrahedrons
  // have equal lengths and side lengths of 1.
  //double const e = 1.32287;

  // Ratio of l[2] to z;
  double const r = 1.5;

  if (z>l[2]*r)
  {
    cout << "large z" << endl;
    return twotetrahedronMInvert(a,ai,b,bi);
  }

  // Test for thin triangle boundary.
  if (approximatelyequal(l[2],l[1]))
  {
    cout << "thin triangle" << endl; 
    return twotetrahedronMInvert(a,ai,b,bi);
  }

  cout << "tet2to3" << endl;

  return tess.tet2to3(a,b);
}