//---------------------------------------------------------------------------


#pragma hdrstop

#include "Trian.h"
#include "MC.h"
#include <math.h>
#include "MT.h"
//---------------------------------------------------------------------------

void Initialize(TProcess *p);
void Refine(TProcess *p);
//---------------------------------------------------------------------------

void Triangulate(TProcess *p){

  p->Stats.All.Start = clock();
  Initialize(p);
  CreateCubes(p);
  if(p->quit) return;
  CreateTriangles(p);
  if(p->quit) return;
  Refine(p);
  
  p->Stats.All.End = clock();
//
}

void Initialize(TProcess *p){

  MCALLBACK(p->Callback, p, S_INIT, A_IN_START)

//Initialize config
//  p->Config.STEP = 0.032;
//  p->Config.SUBDIVIDE = 1;
//  p->Config.MAX_DEPTH = 4;
//  p->Config.RC_RATIO = pow(2, p->Config.MAX_DEPTH);
//  p->Config.BIG_STEP = p->Config.STEP * (double) p->Config.RC_RATIO;

//  p->Config.ALPHA_SUBDIVISON = 0.7;

//  p->Config.ALPHA_DIVIDE = 0.95;
//  p->Config.ALPHA_STEP = 0.95;

//  p->Config.MAX_STEP3 = 0.178;
//  p->Config.MAX_STEP3 = 0.4;
//  p->Config.MAX_STEP3_QUAD = p->Config.MAX_STEP3 * p->Config.MAX_STEP3;
//  p->Config.MIN_STEP3 = 0.07;
//  p->Config.MIN_STEP3 = 0.3;
//  p->Config.MIN_STEP3_QUAD = p->Config.MIN_STEP3 * p->Config.MIN_STEP3;
//  p->Config.MAX_DISTANCE = 0.0;
//  p->Config.MAX_TN = 1.7;

//  p->Config.MIN_VALUE = 0.00000001;

//Initialize hashtables

  p->FoundVertices = (TCubeVertices **) mycalloc(HASHSIZE, sizeof(TCubeVertices *));
  p->FoundCubes = (TCubesHash **) mycalloc(HASHSIZE, sizeof(TCubesHash *));
  p->FoundEdges = (TCubeEdges **) mycalloc(HASHSIZE, sizeof(TCubeEdge *));
  p->TriangleEdges = (TTriangleEdges **) mycalloc(HASHSIZE, sizeof(TTriangleEdges *));

//Clear Lists


  p->ActiveCubes = NULL;


//p.Vertices = NULL;

//p.Patches = NULL;
//p.Triangles = NULL;
//p.TrPatches = NULL;

//Clear Statistics
  p->Stats.Patches.Count = 0;
  p->Stats.Triangles.Count = 0;
  p->Stats.Vertices.Count = 0;

  MCALLBACK(p->Callback, p, S_INIT, A_IN_END)

}

void Refine(TProcess *p){
  T3DPoints *AP = p->Vertices;
  T3DPoint P1;
  T3DPoint P2;
  double Value;

  for(T3DNormals *AN = p->Normals; AN != NULL; AN = AN->Next){
    P1 = AP->Point;
    while(p->Function(P1.X, P1.Y, P1.Z) > 0.0){
      P1 = P1 + (0.01) * AN->Normal;
    }
    P2 = P1 + (-0.01) * AN->Normal;
    while(p->Function(P2.X, P2.Y, P2.Z) < 0.0){
      P2 = P2 + (-0.01) * AN->Normal;
    }


    while(fabs(Value = p->Function(AP->Point.X, AP->Point.Y, AP->Point.Z)) > 0.001){
      if(Value > 0.0){
        P2 = AP->Point;
      }
      else {
        P1 = AP->Point;
      }
      AP->Point = (P1 + P2) / 2.0;
    }

    AN->Normal.X = (p->Function(AP->Point.X + 0.00001, AP->Point.Y, AP->Point.Z)
              - p->Function(AP->Point.X - 0.00001, AP->Point.Y, AP->Point.Z)) / 0.00002;
    AN->Normal.Y = (p->Function(AP->Point.X , AP->Point.Y + 0.00001, AP->Point.Z)
              - p->Function(AP->Point.X , AP->Point.Y - 0.00001, AP->Point.Z)) / 0.00002;
    AN->Normal.Z = (p->Function(AP->Point.X , AP->Point.Y, AP->Point.Z + 0.00001)
              - p->Function(AP->Point.X , AP->Point.Y, AP->Point.Z - 0.00001)) / 0.00002;

    AN->Normal = AN->Normal / sqrt(AN->Normal * AN->Normal);

    AN->Normal = (-1.0) * AN->Normal;

    AP = AP->Next;
  }
}


//---------------------------------------------------------------------------


#pragma package(smart_init)
