/* -----------------------------------------------------------------------------
	GSFramework
	Copyright 2001-2013 Emmanuel Julien. All Rights Reserved.
----------------------------------------------------------------------------- */


		#include <cmath>
		#include <cstring>
		#include "core/iso_surface.h"
		#include "core/geometry.h"
		#include "rand/rand.h"

		using namespace GS;
		using namespace GS::Core;


const int Isosurface::EdgeArray[256]={0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };

const int	Isosurface::TriTable[256][16]=
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};

/*static*/ uint Isosurface::m_MaxNbVertex = 3000;
/*static*/ Isosurface::TVertex* Isosurface::m_TempVertexBuffer = NULL; // Vertex buffer containing the geometry
/*static*/ int Isosurface::m_MaxIdx = 3000;
/*static*/ int* Isosurface::m_TempIdxBuffer = NULL; // index buffer containing the geometry

//***************************************************

void Isosurface::CalcNormalX(float &_x, float &_y, float &_z, Vector4 &_Normal)
{
	float l_Xplus1 = _x+1;
	float l_Yplus1 = _y+1;
	float l_Zplus1 = _z+1;
	float l_Xplus2 = _x+2;
	float l_Xmoins1 = _x-1;
	float l_Ymoins1 = _y-1;
	float l_Zmoins1 = _z-1;
	// so interpolation is on X, the composant x has a special treatment
	_Normal.x = interpolateVal(m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus2,_y,_z)]- m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)]- m_Grid[GridIndex(l_Xmoins1,_y,_z)]);

	_Normal.y = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(l_Xplus1,l_Yplus1,_z)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(_x,l_Ymoins1,_z)], m_Grid[GridIndex(l_Xplus1,l_Ymoins1,_z)]);

	_Normal.z = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(l_Xplus1,_y,l_Zplus1)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zmoins1)], m_Grid[GridIndex(l_Xplus1,_y,l_Zmoins1)]);
}

//***************************************************

void Isosurface::CalcNormalY(float &_x, float &_y, float &_z, Vector4 &_Normal)
{
	float l_Xplus1 = _x+1;
	float l_Yplus1 = _y+1;
	float l_Zplus1 = _z+1;
	float l_Yplus2 = _y+2;
	float l_Xmoins1 = _x-1;
	float l_Ymoins1 = _y-1;
	float l_Zmoins1 = _z-1;
	// so interpolation is on Y, the composant Y has a special treatment
	_Normal.y = interpolateVal(m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus2,_z)]- m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)]- m_Grid[GridIndex(_x,l_Ymoins1,_z)]);

	_Normal.x = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(l_Xplus1,l_Yplus1,_z)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(l_Xmoins1,_y,_z)], m_Grid[GridIndex(l_Xmoins1,l_Yplus1,_z)]);

	_Normal.z = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(_x,l_Yplus1,l_Zplus1)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(_x,_y,l_Zmoins1)], m_Grid[GridIndex(_x,l_Yplus1,l_Zmoins1)]);
}

//***************************************************

void Isosurface::CalcNormalZ(float &_x, float &_y, float &_z, Vector4 &_Normal)
{
	float l_Xplus1 = _x+1;
	float l_Yplus1 = _y+1;
	float l_Zplus1 = _z+1;
	float l_Zplus2 = _z+2;
	float l_Xmoins1 = _x-1;
	float l_Ymoins1 = _y-1;
	float l_Zmoins1 = _z-1;
	// so interpolation is on Z, the composant Z has a special treatment
	_Normal.z = interpolateVal(m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus2)]- m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)]- m_Grid[GridIndex(_x,_y,l_Zmoins1)]);

	_Normal.y = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(_x,l_Yplus1,_z)], m_Grid[GridIndex(_x,l_Yplus1,l_Zplus1)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(_x,l_Ymoins1,_z)], m_Grid[GridIndex(_x,l_Ymoins1,l_Zplus1)]);

	_Normal.x = interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(l_Xplus1,_y,_z)], m_Grid[GridIndex(l_Xplus1,_y,l_Zplus1)]) -
		interpolateVal(m_Grid[GridIndex(_x,_y,_z)], m_Grid[GridIndex(_x,_y,l_Zplus1)], m_Grid[GridIndex(l_Xmoins1,_y,_z)], m_Grid[GridIndex(l_Xmoins1,_y,l_Zplus1)]);
}

//***************************************************

// function permit to linear interpolate between vector
void Isosurface::interpolateVect(Vector4 &_Vect1, Vector4 &_Vect2, float &_Val1, float &_Val2, Vector4&_Vect)
{
	// _Vect1 and _Vect2 are extremities (with _val1 and _val2)
	// _Vect is the position to find

// 	if(fabsf(m_IsoValue - _Val1) < 0.00001)
// 	{
// 		// to don't have arround error
// 		_Vect = _Vect1;
// 	}
// 
// 	if(fabsf(m_IsoValue - _Val2) < 0.00001)
// 	{
// 		// to don't have arround error
// 		_Vect = _Vect2;
// 	}
// 
// 	if(fabsf(_Val1 - _Val2) < 0.00001)
// 	{
// 		// to don't have arround error
// 		_Vect = _Vect1;
// 	}

	float l_Coef = (m_IsoValue - _Val1) / (_Val2 - _Val1);

	_Vect.x = _Vect1.x + l_Coef * (_Vect2.x - _Vect1.x);
	_Vect.y = _Vect1.y + l_Coef * (_Vect2.y - _Vect1.y);
	_Vect.z = _Vect1.z + l_Coef * (_Vect2.z - _Vect1.z);

}

//***************************************************

// function permit to linear interpolate between val
float Isosurface::interpolateVal(float &_Val1, float &_Val2, float _Val_cible1, float _Val_cible2)
{
	if(fabsf(m_IsoValue - _Val1) < 0.00001)
	{
		return _Val_cible1;
	}

	if(fabsf(m_IsoValue - _Val2) < 0.00001)
	{
		return _Val_cible2;
	}

	if(fabsf(_Val1 - _Val2) < 0.00001)
	{
		return _Val_cible1;
	}

	if(m_IsoValue - _Val1 != 0)
	{
		float l_Coef = (m_IsoValue - _Val1) / (_Val2 - _Val1);

		return  _Val_cible1 + l_Coef * (_Val_cible2 - _Val_cible1);
	}
	else
		return _Val_cible1;
}

//***************************************************

// Set the value of the vertex
void Isosurface::EvalPos(float &_x, float &_y, float &_z, Vector4 &_CasePosition)
{
	_CasePosition.x = (_x/* - m_NbGridCase.x*0.5f*/) * m_CaseSizeDivNbCase.x /*+m_Pos.x*/;
	_CasePosition.y = (_y/* - m_NbGridCase.y*0.5f*/) * m_CaseSizeDivNbCase.y /*+m_Pos.y*/;
	_CasePosition.z = (_z/* - m_NbGridCase.z*0.5f*/) * m_CaseSizeDivNbCase.z /*+m_Pos.z*/;
}

//***************************************************

// return the good float with the 3d parameter in the grid
int Isosurface::GridIndex(float &_x, float &_y, float &_z)const
{
	return (int)(_x+_y*m_NbGridCase.x+_z*m_NbGridCase.x*m_NbGridCase.y);
}

//***************************************************

/// calcul the polygon for one case, and return the number of triangle for the polygon
int		Isosurface::CalculPolygon(float &_x, float &_y, float &_z, STriangle* _TriangleList)
{
	float l_GridValues[8];		// values for each vertex of the case
	Vector4  l_GridPositions[8];	// case vertex position

	// vertex List for the vertex in the final polygon
	Vector4 l_ListVertex[12];

	float l_Xplus1 = _x+1;
	float l_Yplus1 = _y+1;
	float l_Zplus1 = _z+1;

	// calcul values of the vertex of the case
	l_GridValues[0] = m_Grid[GridIndex(_x,_y,_z)];
	l_GridValues[1] = m_Grid[GridIndex(l_Xplus1,_y,_z)];
	l_GridValues[2] = m_Grid[GridIndex(l_Xplus1,_y,l_Zplus1)];
	l_GridValues[3] = m_Grid[GridIndex(_x,_y,l_Zplus1)];
	l_GridValues[4] = m_Grid[GridIndex(_x,l_Yplus1,_z)];
	l_GridValues[5] = m_Grid[GridIndex(l_Xplus1,l_Yplus1,_z)];
	l_GridValues[6] = m_Grid[GridIndex(l_Xplus1,l_Yplus1,l_Zplus1)];
	l_GridValues[7] = m_Grid[GridIndex(_x,l_Yplus1,l_Zplus1)];

	// Calcul vertex Position
	EvalPos(_x, _y, _z, l_GridPositions[0]);
	EvalPos(l_Xplus1, _y, _z, l_GridPositions[1]);
	EvalPos(l_Xplus1, _y, l_Zplus1, l_GridPositions[2]);
	EvalPos(_x, _y, l_Zplus1, l_GridPositions[3]);
	EvalPos(_x, l_Yplus1, _z, l_GridPositions[4]);
	EvalPos(l_Xplus1, l_Yplus1, _z, l_GridPositions[5]);
	EvalPos(l_Xplus1, l_Yplus1, l_Zplus1, l_GridPositions[6]);
	EvalPos(_x, l_Yplus1, l_Zplus1, l_GridPositions[7]);

	// find the index in the edge Array to know wich side the surface intersect
	int l_Index = 0;
	if(l_GridValues[0] < m_IsoValue)
		l_Index |= 1;		// put the bit 0 at 1

	if(l_GridValues[1] < m_IsoValue)
		l_Index |= 2;		
	if(l_GridValues[2] < m_IsoValue)
		l_Index |= 4;		
	if(l_GridValues[3] < m_IsoValue)
		l_Index |= 8;		
	if(l_GridValues[4] < m_IsoValue)
		l_Index |= 16;		
	if(l_GridValues[5] < m_IsoValue)
		l_Index |= 32;		
	if(l_GridValues[6] < m_IsoValue)
		l_Index |= 64;		
	if(l_GridValues[7] < m_IsoValue)
		l_Index |= 128;		

	// Calcul vertex position where the surface intersect the case
	if(EdgeArray[l_Index] == 0)		// the case is out the surface
		return 0;
	if(EdgeArray[l_Index] & 1)
		interpolateVect(l_GridPositions[0], l_GridPositions[1], l_GridValues[0], l_GridValues[1], l_ListVertex[0]);
	if(EdgeArray[l_Index] & 2)
		interpolateVect(l_GridPositions[1], l_GridPositions[2], l_GridValues[1], l_GridValues[2], l_ListVertex[1]);
	if(EdgeArray[l_Index] & 4)
		interpolateVect(l_GridPositions[2], l_GridPositions[3], l_GridValues[2], l_GridValues[3], l_ListVertex[2]);
	if(EdgeArray[l_Index] & 8)
		interpolateVect(l_GridPositions[3], l_GridPositions[0], l_GridValues[3], l_GridValues[0], l_ListVertex[3]);
	if(EdgeArray[l_Index] & 16)
		interpolateVect(l_GridPositions[4], l_GridPositions[5], l_GridValues[4], l_GridValues[5], l_ListVertex[4]);
	if(EdgeArray[l_Index] & 32)
		interpolateVect(l_GridPositions[5], l_GridPositions[6], l_GridValues[5], l_GridValues[6], l_ListVertex[5]);
	if(EdgeArray[l_Index] & 64)
		interpolateVect(l_GridPositions[6], l_GridPositions[7], l_GridValues[6], l_GridValues[7], l_ListVertex[6]);
	if(EdgeArray[l_Index] & 128)
		interpolateVect(l_GridPositions[7], l_GridPositions[4], l_GridValues[7], l_GridValues[4], l_ListVertex[7]);
	if(EdgeArray[l_Index] & 256)
		interpolateVect(l_GridPositions[0], l_GridPositions[4], l_GridValues[0], l_GridValues[4], l_ListVertex[8]);
	if(EdgeArray[l_Index] & 512)
		interpolateVect(l_GridPositions[1], l_GridPositions[5], l_GridValues[1], l_GridValues[5], l_ListVertex[9]);
	if(EdgeArray[l_Index] & 1024)
		interpolateVect(l_GridPositions[2], l_GridPositions[6], l_GridValues[2], l_GridValues[6], l_ListVertex[10]);
	if(EdgeArray[l_Index] & 2048)
		interpolateVect(l_GridPositions[3], l_GridPositions[7], l_GridValues[3], l_GridValues[7], l_ListVertex[11]);

	// Calcul the triangles
	int l_NbTriangles = 0;
	for(int i=0; TriTable[l_Index][i]!=-1; i+=3)
	{
		_TriangleList[l_NbTriangles].m_Vertex[0] = l_ListVertex[TriTable[l_Index][i]];
		_TriangleList[l_NbTriangles].m_Num[0] = TriTable[l_Index][i];

		_TriangleList[l_NbTriangles].m_Vertex[1] = l_ListVertex[TriTable[l_Index][i+1]];
		_TriangleList[l_NbTriangles].m_Num[1] = TriTable[l_Index][i+1];

		_TriangleList[l_NbTriangles].m_Vertex[2] = l_ListVertex[TriTable[l_Index][i+2]];
		_TriangleList[l_NbTriangles].m_Num[2] = TriTable[l_Index][i+2];

		++l_NbTriangles;
	}

	return l_NbTriangles;
}

//***************************************************

// draw the triangle in the cellule x,y, z
void		Isosurface::RenderCell(uint &vtx_count, float &_x, float &_y, float &_z)
{
//	nVector l_U, l_V;	// 2 vector permit to define 1 face
	//	float l_Color[3];	//  color of the current vertex

	Vector4 *l_Pos;		// vertex pos

	// Polygon triangle
	STriangle l_TriangleList[12];

	int l_NbFace= CalculPolygon(_x, _y, _z, l_TriangleList);	// nombre de faces in the current polygone

	for( int Face=0; Face< l_NbFace; ++Face)
	{
		m_Mutex.Lock();

		// compute the normal vector
		for(int i=0; i< 3; ++i)
		{
			// calcul of the normale for each vertex for each new triangle
		//	nVector l_Normal; // current Normale

			l_Pos = &l_TriangleList[Face].m_Vertex[i];

		/*	if(_x>0 && _y>0 && _z>0)
			{
				float l_Xplus1 = _x+1;
				float l_Yplus1 = _y+1;
				float l_Zplus1 = _z+1;

				switch(l_TriangleList[Face].m_Num[i])
				{
				case 0:
					CalcNormalX(_x, _y, _z, l_Normal);
					break;
				case 1:
					CalcNormalZ(l_Xplus1, _y, _z, l_Normal);
					break;
				case 2:
					CalcNormalX(_x, _y, l_Zplus1, l_Normal);
					break;
				case 3:
					CalcNormalZ(_x, _y, _z, l_Normal);
					break;
				case 4:
					CalcNormalX(_x, l_Yplus1, _z, l_Normal);
					break;
				case 5:
					CalcNormalZ(l_Xplus1, l_Yplus1, _z, l_Normal);
					break;
				case 6:
					CalcNormalX(_x, l_Yplus1, l_Zplus1, l_Normal);
					break;
				case 7:
					CalcNormalZ(_x, l_Yplus1, _z, l_Normal);
					break;
				case 8:
					CalcNormalY(_x, _y, _z, l_Normal);
					break;
				case 9:
					CalcNormalY(l_Xplus1, _y, _z, l_Normal);
					break;
				case 10:
					CalcNormalY(l_Xplus1, _y, l_Zplus1, l_Normal);
					break;
				case 11:
					CalcNormalY(_x, _y, l_Zplus1, l_Normal);
					break;
				}
			}
*/
			// if don't use the cube mapping
			//		if(texture == NULL)
			{
				// compute vertex color
				/*	l_Color[0] = (l_Pos->x/ m_GridSize.x +1) *0.5f;
				l_Color[1] = (l_Pos->y/ m_GridSize.y +1) *0.5f;
				l_Color[2] = (l_Pos->z/ m_GridSize.z +1) *0.5f;

				glMaterialfv(GL_FRONT, GL_SPECULAR, l_Color);
				glMaterialf(GL_FRONT, GL_SHININESS, 50.0f);
				glMaterialfv(GL_FRONT, GL_AMBIENT, l_Color);
				glMaterialfv(GL_FRONT, GL_DIFFUSE, l_Color);*/
			}

			// re- Normalise normal
		//	l_Normal = l_Normal.Normalize();

			//			glNormal3f(-l_Normal.x, -l_Normal.y, -l_Normal.z);
			//			glVertex3f(l_Pos->x, l_Pos->y, l_Pos->z);

			if(vtx_count >= m_MaxNbVertex)
			{ // not enough vertice increase the temp array
				TVertex* l_NewTempVertexBuffer = new TVertex[m_MaxNbVertex+3000];
				memcpy(l_NewTempVertexBuffer, m_TempVertexBuffer, sizeof(TVertex)*m_MaxNbVertex);
				delete []m_TempVertexBuffer;
				m_MaxNbVertex += 3000;
				m_TempVertexBuffer = l_NewTempVertexBuffer;
			}

			if(m_CountIdxBuffer >= m_MaxIdx)
			{ // not enough vertice increase the temp array
				int* l_NewTempVertexBuffer = new int[m_MaxIdx+3000];
				memcpy(l_NewTempVertexBuffer, m_TempIdxBuffer, sizeof(int)*m_MaxIdx);
				delete []m_TempIdxBuffer;
				m_MaxIdx += 3000;
				m_TempIdxBuffer = l_NewTempVertexBuffer;
			}			
/*
			m_TempVertexBuffer[vtx.GetCount()].nx = -l_Normal.x;
			m_TempVertexBuffer[vtx.GetCount()].ny = -l_Normal.y;
			m_TempVertexBuffer[vtx.GetCount()].nz = -l_Normal.z;
*/

			// ugly stuff to find id
			int l_IdFind = -1;
			for(uint id = 0; id<vtx_count; ++id)
				if(		m_TempVertexBuffer[id].x == l_Pos->x
					&&	m_TempVertexBuffer[id].y == l_Pos->y
					&&	m_TempVertexBuffer[id].z == l_Pos->z)
				{
					l_IdFind = id;
					break;
				}

			if(l_IdFind == -1)	
			{
				m_TempVertexBuffer[vtx_count].x = l_Pos->x;
				m_TempVertexBuffer[vtx_count].y = l_Pos->y;
				m_TempVertexBuffer[vtx_count].z = l_Pos->z;

				l_IdFind = vtx_count;
				++vtx_count;
			}
			
			m_TempIdxBuffer[m_CountIdxBuffer++] = l_IdFind;

		}	
		m_Mutex.Unlock();	
	}
}

//*********************************************

void		ClipVector(Vector4 &_VectorIndexClip, Vector4 &_MaxIndex)
{
	if(_VectorIndexClip.x < 0)
		_VectorIndexClip.x = 0;
	if(_VectorIndexClip.x >= _MaxIndex.x)
		_VectorIndexClip.x = _MaxIndex.x-1;
	if(_VectorIndexClip.y < 0)
		_VectorIndexClip.y = 0;
	if(_VectorIndexClip.y >= _MaxIndex.y)
		_VectorIndexClip.y = _MaxIndex.y-1;
	if(_VectorIndexClip.z < 0)
		_VectorIndexClip.z = 0;
	if(_VectorIndexClip.z >= _MaxIndex.z)
		_VectorIndexClip.z = _MaxIndex.z-1;
}

//*********************************************
void		Isosurface::ComputeMetaball(int nb_metaball, Vector4* pos_metaball, float* value_metaball)
{
	Vector4 l_Vect;	// vector between the ball center and the case of the grid
	Vector4 l_Pos, l_PosMetaBall;		// case center position
	float l_Dist;	// distance between ball center and the case
	float l_TempRayon, l_TempRayon2, l_TempDistDivRayon;  // rayon of the metaball and distance divide by the rayon
	int l_Index, l_IndexMax = (int)(m_NbGridCase.x*m_NbGridCase.y*m_NbGridCase.z);

	float x, y, z;
	int i;		// count paramater
	Vector4 l_PosMetInGridDeb, l_PosMetInGridFin;

	Vector4 l_TempSize(m_GridSize.x/(m_NbGridCase.x-2), m_GridSize.y/(m_NbGridCase.y-2), m_GridSize.z/(m_NbGridCase.z-2));
	Vector4 l_TempSizeInverse((m_NbGridCase.x-2)/m_GridSize.x, (m_NbGridCase.y-2)/m_GridSize.y, (m_NbGridCase.z-2)/m_GridSize.z);

	// clean all the grid
	int l_NbCase = (int)(m_NbGridCase.x*m_NbGridCase.y*m_NbGridCase.z);
	memset(m_Grid, 0, sizeof(float)*l_NbCase);

	//	memset(m_GridMetaBall->m_TabIndexValid, 0, sizeof(nVector)*m_NbGridCase.x*m_NbGridCase.y*m_NbGridCase.z);
	int l_NbIndexValid = 0;
	
	float l_MetaballIntensity = 1.0f;

	for(i=0; i< nb_metaball; ++i)
	{
		l_PosMetaBall = pos_metaball[i];
		l_TempRayon = value_metaball[i];
		l_TempRayon2 = l_TempRayon*l_TempRayon;

 		if((l_PosMetaBall.x - l_TempRayon < m_Pos.x + m_GridSize.x && l_PosMetaBall.x + l_TempRayon > m_Pos.x)
 			&& (l_PosMetaBall.y - l_TempRayon < m_Pos.y + m_GridSize.y && l_PosMetaBall.y + l_TempRayon > m_Pos.y)
 			&& (l_PosMetaBall.z - l_TempRayon < m_Pos.z + m_GridSize.z && l_PosMetaBall.z + l_TempRayon > m_Pos.z))
		{
			l_PosMetInGridDeb.x = (((l_PosMetaBall.x - l_TempRayon)- m_Pos.x)*l_TempSizeInverse.x);
			l_PosMetInGridDeb.y = (((l_PosMetaBall.y - l_TempRayon)- m_Pos.y)*l_TempSizeInverse.y);
			l_PosMetInGridDeb.z = (((l_PosMetaBall.z - l_TempRayon)- m_Pos.z)*l_TempSizeInverse.z);

			l_PosMetInGridFin.x = (((l_PosMetaBall.x + l_TempRayon)- m_Pos.x)*l_TempSizeInverse.x);
			l_PosMetInGridFin.y = (((l_PosMetaBall.y + l_TempRayon)- m_Pos.y)*l_TempSizeInverse.y);
			l_PosMetInGridFin.z = (((l_PosMetaBall.z + l_TempRayon)- m_Pos.z)*l_TempSizeInverse.z);

			ClipVector(l_PosMetInGridDeb, m_NbGridCase);
			ClipVector(l_PosMetInGridFin, m_NbGridCase);

			// find only the voxel touch by the metaball

			for(z=l_PosMetInGridDeb.z; z< l_PosMetInGridFin.z; ++z)
				for(y=l_PosMetInGridDeb.y; y< l_PosMetInGridFin.y; ++y)
					for(x=l_PosMetInGridDeb.x; x< l_PosMetInGridFin.x; ++x)
					{
						// calcul distance between the case and the center of the metal ball
						l_Vect.x = (x*l_TempSize.x) + m_Pos.x - l_PosMetaBall.x;
						l_Vect.y = (y*l_TempSize.y) + m_Pos.y - l_PosMetaBall.y;
						l_Vect.z = (z*l_TempSize.z) + m_Pos.z - l_PosMetaBall.z;

						l_Dist = l_Vect.Len2();
						// take place for the metaball in the grid
						if(l_Dist <= l_TempRayon2)
							//if(l_Dist <= l_TempRayon2  && (y*l_TempSize.y) + m_Pos.y - l_PosMetaBall.y <= l_TempRayon2)
						{
							l_Index = (int)(x+ y*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y);

							if(l_Index >=0 && l_Index <l_IndexMax)
							{
								l_TempDistDivRayon = 1 - l_Dist / l_TempRayon2;
								m_Grid[l_Index] += l_MetaballIntensity * l_TempDistDivRayon* l_TempDistDivRayon;						
							}
						}
					}
		}
	}

	// to go faster after, just stock the used voxel
	// add the index valid , to don't look the other empty voxels
	for(z=0; z< m_NbGridCase.z-1; ++z)
		for(y=0; y< m_NbGridCase.y-1; ++y)
			for(x=0; x< m_NbGridCase.x-1; ++x)
			{
				if(m_Grid[(int)(x+ y*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)((x+1)+ y*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)((x+1)+ (y+1)*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)((x+1)+ (y+1)*m_NbGridCase.x + (z+1)*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)((x+1)+ y*m_NbGridCase.x + (z+1)*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)(x+ (y+1)*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)(x+ (y+1)*m_NbGridCase.x + (z+1)*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f
					|| m_Grid[(int)(x+ y*m_NbGridCase.x + (z+1)*m_NbGridCase.x*m_NbGridCase.y)] > 0.0f)
				{
					m_TabIndexValid[l_NbIndexValid].x = x;
					m_TabIndexValid[l_NbIndexValid].y = y;
					m_TabIndexValid[l_NbIndexValid].z = z;

					++l_NbIndexValid;
				}		
			}

	m_NbIndexValid = l_NbIndexValid;
}

/// Triangularize the iso-surface
void        Isosurface::Triangularize(Geometry* geometry, int nb_metaball, Vector4* pos_metaball, float* value_metaball)
{
// 	point[0].x += 0.1f;
// 	if(point[0].x > 25)
// 		point[0].x = 0;
// 	point[1].z -= 0.1f;
// 	point[2].y -= 0.1f;
// 	if(point[1].z < -5)
// 		point[1].z = 12;
// 	if(point[2].y < -15)
// 		point[2].y = 15;

	ComputeMetaball(nb_metaball, pos_metaball, value_metaball);

	m_CountIdxBuffer = 0;

    uint vtx_count = 0;

	#pragma omp parallel for
    for(int i=0; i< m_NbIndexValid; ++i)
    {
        RenderCell(vtx_count, m_TabIndexValid[i].x, m_TabIndexValid[i].y, m_TabIndexValid[i].z);
    }

    if(vtx_count >0 && vtx_count <m_MaxNbVertex)
    {
        // initialize Cube vertex
        // nombre de faces et de sommets

        geometry->AllocateVertex(vtx_count);
        geometry->AllocatePolygon(m_CountIdxBuffer/3);
	//	geometry->vtx_normal = new nVector[vtx.GetCount()];

        // fill the table
		#pragma omp parallel for
        for(int i=0; i< static_cast<int> (vtx_count); ++i)
        {
				geometry->vtx[i] = Vector4(m_TempVertexBuffer[i].x, m_TempVertexBuffer[i].y, m_TempVertexBuffer[i].z);
		//		geometry->vtx_normal[i*3 + 2-n] = nVector(m_TempVertexBuffer[i*3 + n].nx, m_TempVertexBuffer[i*3 + n].ny, m_TempVertexBuffer[i*3 + n].nz);
		}

		#pragma omp parallel for
        for(int i=0; i< static_cast<int> (geometry->pol.GetCount()); ++i)
        {
            geometry->pol[i].vtx_count = 3;
            geometry->pol[i].material = 0;
        }

        geometry->AllocatePolygonBinding();

		#pragma omp parallel for
        for(int i=0; i< static_cast<int> (geometry->pol.GetCount()); ++i)
        {
            geometry->pol[i].vtx_count = 3;
			for (int n = 0; n < 3; ++n)
				geometry->pol[i].binding[2-n] =  m_TempIdxBuffer[i * 3 + n];
            geometry->pol[i].material = 0;
        }

		geometry->ComputeVertexNormal(true);
	}
}

//********************************************************************

bool		Isosurface::Init(const Vector4 &_pos, const Vector4 &size, const Vector4 &step)
{
	m_Pos = _pos;

	m_NbGridCase = step;
	m_GridSize = size;

	delete []m_Grid;

	m_Grid = new float[int(ceil(m_NbGridCase.x)*ceil(m_NbGridCase.y)*ceil(m_NbGridCase.z))];

	// clean all the grid
	int l_NbCase = (int)(m_NbGridCase.x*m_NbGridCase.y*m_NbGridCase.z);
	memset(m_Grid, 0, sizeof(float)*l_NbCase);

	delete []m_TabIndexValid;

	m_TabIndexValid = new Vector4[int(ceil(m_NbGridCase.x)*ceil(m_NbGridCase.y)*ceil(m_NbGridCase.z))];

	m_NbIndexValid = 0;

	for(float z=0; z< m_NbGridCase.z-1; ++z)
		for(float y=0; y< m_NbGridCase.y-1; ++y)
			for(float x=0; x< m_NbGridCase.x-1; ++x)
			{
				//if(m_Grid[x+ y*m_NbGridCase.x + z*m_NbGridCase.x*m_NbGridCase.y])
				{
					m_TabIndexValid[m_NbIndexValid].x = x;
					m_TabIndexValid[m_NbIndexValid].y = y;
					m_TabIndexValid[m_NbIndexValid].z = z;

					++m_NbIndexValid;
				}		
			}

	m_CaseSizeDivNbCase = Vector4(m_GridSize.x/m_NbGridCase.x, m_GridSize.y/m_NbGridCase.y, m_GridSize.z/m_NbGridCase.z);

	return true;
}

/// Get the iso-surface field size
void		Isosurface::GetFieldSize(int &x, int &y, int &z)
{
	x = (int)(m_GridSize.x);
	y = (int)(m_GridSize.y);
	z = (int)(m_GridSize.z);
}

/// Get the iso-surface field
void		Isosurface::GetField(float *_Grid)
{
	_Grid = m_Grid;
}

//-------------------------------------------------------
Isosurface::Isosurface()
{
	m_Pos = Vector4(0,0,0);
	m_Grid = NULL;
	m_TabIndexValid = NULL;

	m_IsoValue =0.99f;

	m_CountIdxBuffer = 0;
	m_TempIdxBuffer = new int[m_MaxIdx];

	m_TempVertexBuffer = new TVertex[m_MaxNbVertex];
}

//-------------------------
Isosurface::~Isosurface()
{	
	delete []m_Grid;
	delete []m_TempVertexBuffer;
}