using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Windows.Forms;
using Tao.OpenGl;
using Tao.Platform.Windows;
using MyGeometry;
using System.IO;

using System.Runtime.InteropServices;


namespace FunctionalSubstructure
{
	public unsafe class MeshView : SimpleOpenGlControl
	{

		[DllImport("opengl32", EntryPoint = "wglUseFontBitmaps", CallingConvention=CallingConvention.Winapi)]
		public static extern bool wglUseFontBitmaps(
		IntPtr hDC,
		[MarshalAs(UnmanagedType.U4)] UInt32 first,
		[MarshalAs(UnmanagedType.U4)] UInt32 count,
		[MarshalAs(UnmanagedType.U4)] UInt32 listBase
		);

		[DllImport("GDI32.DLL", EntryPoint = "SelectObject",
		CallingConvention = CallingConvention.Winapi)]
		public static extern IntPtr SelectObject(
		[In] IntPtr hDC,
		[In] IntPtr font
		);


		public interface IMeshDisplay
		{
			void Display();
			void SetData();
		}
		public enum EnumOperationMode { Viewing, Selection};
		public delegate void OnPickingHandler(object sender, int index);


		// private fields
		private EnumOperationMode currentMode = EnumOperationMode.Viewing;


		// font
		private bool initFont = false;
		public UInt32 fontBase = 0;
		private void BuildFont(PaintEventArgs pe)
		{
			IntPtr dc = pe.Graphics.GetHdc();
			IntPtr oldFontH = IntPtr.Zero;
			System.Drawing.Font font =
				new Font(
				"Arial",
				10F,
				System.Drawing.FontStyle.Bold,
				System.Drawing.GraphicsUnit.Point,
				((System.Byte)(0)));

			fontBase = (uint)Gl.glGenLists(128);

			IntPtr fontH = font.ToHfont();
			oldFontH = SelectObject(dc, fontH);

			bool ret = wglUseFontBitmaps(
				dc,
				0,
				128,
				fontBase);

			SelectObject(dc, oldFontH);						// Selects The Font We Want

			pe.Graphics.ReleaseHdc(dc);

			if (!ret) throw new Exception();
		}

		// display properties
		private Color pointColor = Color.Red;
		private Color lineColor = Color.Black;

		public int selectedFace = -1;
		public int gindex = 0;
		public bool showWireFrame = false;


		// for viewing
		private Trackball ball = new Trackball(100, 100);
		private double scaleRatio;
		private Vector2d currMousePosition = new Vector2d();
		private Vector2d mouseDownPosition = new Vector2d();
		private bool isMouseDown = false;
		private Matrix4d currTransformation = Matrix4d.IdentityMatrix();
		private bool useGLOrthoProjection = false;



		// public properties
		public EnumOperationMode CurrentMode
		{
			get { return currentMode; }
			set { currentMode = value; }
		}


		public MeshView()
		{
			InitializeComponent();

			this.AutoMakeCurrent = true;
			this.InitializeContexts();
			this.InitOpenGL();

			LoadSpatialGraphs();
			LoadComponentSlots();
			LoadGroundContactSlots();
			LoadSymmetryTriplets();
			ComputeGraphNodeFunctionalities();
		}
		private void InitializeComponent()
		{
			this.SuspendLayout();
			this.ResumeLayout(false);
		}

		protected override void OnSizeChanged(EventArgs e)
		{
			base.OnSizeChanged(e);

			Size s = this.Size;
			if (s.Width == 0 || s.Height == 0) return;

			this.scaleRatio = (s.Width > s.Height) ? s.Height : s.Width;
		//	this.InitMatrix();
			this.ball.SetBounds(s.Width * 1.0, s.Height * 1.0);
			this.Refresh();
		}
		protected override void OnMouseDown(System.Windows.Forms.MouseEventArgs e)
		{
			base.OnMouseDown(e);

			Program.CurrentActiveView = this;



			this.currMousePosition = new Vector2d(e.X, this.Height - e.Y);
			this.mouseDownPosition = currMousePosition;
			this.isMouseDown = true;

			switch (this.currentMode)
			{
				case EnumOperationMode.Viewing:
					switch (e.Button)
					{
						case MouseButtons.Left: ball.Click(currMousePosition, Trackball.MotionType.Rotation); break;
						case MouseButtons.Middle: ball.Click(currMousePosition / this.scaleRatio, Trackball.MotionType.Pan); break;
						case MouseButtons.Right: ball.Click(currMousePosition, Trackball.MotionType.Scale); break;
					}
					break;
			}
		}
		protected override void OnMouseMove(System.Windows.Forms.MouseEventArgs e)
		{
			base.OnMouseMove(e);

			this.currMousePosition = new Vector2d(e.X, this.Height - e.Y);

			switch (this.currentMode)
			{
				case EnumOperationMode.Viewing:
					if (isMouseDown)
					{
						switch (e.Button)
						{
							case MouseButtons.Left: ball.Drag(currMousePosition); break;
							case MouseButtons.Middle: ball.Drag(currMousePosition / this.scaleRatio); break;
							case MouseButtons.Right: ball.Drag(currMousePosition); break;
						}
						this.Refresh();
					}
					break;

			}
		}
		protected override void OnMouseUp(System.Windows.Forms.MouseEventArgs e)
		{
			base.OnMouseUp(e);

			this.currMousePosition = new Vector2d(e.X, this.Height - e.Y);
			this.isMouseDown = false;

			switch (this.currentMode)
			{
				case EnumOperationMode.Viewing:
					if (currMousePosition == mouseDownPosition) break;
					Matrix4d m = ball.GetMatrix();
					this.currTransformation = m * this.currTransformation;
					this.ball.End();

					this.Refresh();
					break;
			}
		}
		protected override void OnPaint(PaintEventArgs e)
		{

			if (!initFont)
			{
				BuildFont(e);
				initFont = true;
			}


			if (this.DesignMode == true)
			{
				base.OnPaint(e); 
				return;
			}

			this.MakeCurrent();
			this.InitMatrix();

			
			Matrix4d m = ball.GetMatrix() * currTransformation;
			Gl.glClearColor(1.0f, 1.0f, 1.0f, 0.0f);
			Gl.glClear(Gl.GL_COLOR_BUFFER_BIT | Gl.GL_DEPTH_BUFFER_BIT);
			Gl.glMatrixMode(Gl.GL_MODELVIEW);
			if (!this.useGLOrthoProjection)
			{
				Gl.glPushMatrix();
				Gl.glMultMatrixd(m.Transpose().ToArray());
			}
			else
			{
				Gl.glLoadMatrixd(m.Transpose().ToArray());
			}

			Gl.glEnable(Gl.GL_POLYGON_OFFSET_FILL);

			
			Gl.glTranslated(0, -0.5, 0);


			if (this.synthesizedCuboidMeshComponents.Count > 0 && this._curr_mesh_displayindex < this.synthesizedCuboidMeshComponents.Count)
			{
				int j = this._curr_mesh_displayindex;

				List<CuboidMeshComponent> cmcs = this.synthesizedCuboidMeshComponents[j];
				foreach (CuboidMeshComponent cmc in cmcs)
				{
					if (this.showWireFrame)
						this.DrawWireframe(cmc.mesh, Color.Green);
					this.DrawTransparent(cmc.mesh,Color.Tan);
				}
	
			}


			if (!this.useGLOrthoProjection)
				Gl.glPopMatrix();

			this.SwapBuffers();
		}


		private void InitOpenGL()
		{
		//	this.MakeCurrent();
			
			Color c = SystemColors.Control;
			Gl.glClearColor(c.R / 255f, c.G / 255f, c.B / 255f, 0.0f);     // BackGround Color         
			Gl.glShadeModel(Gl.GL_SMOOTH);         // Set Smooth Shading                 
			Gl.glClearDepth(1.0f);                 // Depth buffer setup             
			Gl.glEnable(Gl.GL_DEPTH_TEST);         // Enables Depth Testing             
			Gl.glDepthFunc(Gl.GL_LEQUAL);             // The Type Of Depth Test To Do     
			Gl.glHint(Gl.GL_PERSPECTIVE_CORRECTION_HINT, Gl.GL_NICEST);     /* Really Nice Perspective Calculations */
			Gl.glEnable(Gl.GL_CULL_FACE);
			Gl.glPolygonOffset(1f, 1f);
			Gl.glBlendFunc(Gl.GL_SRC_ALPHA, Gl.GL_ONE_MINUS_SRC_ALPHA);

			Gl.glEnable(Gl.GL_POINT_SMOOTH);
			Gl.glEnable(Gl.GL_LINE_SMOOTH);
			Gl.glEnable(Gl.GL_BLEND);

			float[] LightDiffuse = { 1f, 1f, 1f, 0f };
			float[] LightSpecular = { 0.7f, 0.7f, 0.7f, 1f };
			float[] SpecularRef = { 0.5f, 0.5f, 0.5f, 1f };
			Gl.glLightfv(Gl.GL_LIGHT0, Gl.GL_DIFFUSE, LightDiffuse);
			Gl.glLightfv(Gl.GL_LIGHT0, Gl.GL_SPECULAR, LightSpecular);
			Gl.glEnable(Gl.GL_LIGHT0);
			Gl.glEnable(Gl.GL_COLOR_MATERIAL);
			Gl.glColorMaterial(Gl.GL_FRONT, Gl.GL_DIFFUSE);
			Gl.glMaterialfv(Gl.GL_FRONT, Gl.GL_SPECULAR, SpecularRef);
			Gl.glMateriali(Gl.GL_FRONT, Gl.GL_SHININESS, 128);

			Gl.glClearColor(0, 0, 0, 0);
			Gl.glMatrixMode(Gl.GL_PROJECTION);
			Gl.glLoadIdentity();
			Gl.glOrtho(0.0, 1.0, 0.0, 1.0, -1.0, 1.0);

	//		Size s = this.Size;
	//		ball = new Trackball(s.Width * 1.0, s.Height * 1.0);
		}
		private void InitMatrix_Ortho()
		{
			Size s = this.Size;
			double w = s.Width;
			double h = s.Height;

		//	this.MakeCurrent();
			
			Gl.glViewport(0, 0, s.Width, s.Height);
			Gl.glMatrixMode(Gl.GL_PROJECTION);
			Gl.glLoadIdentity();

			if (w > h)
			{
				double ratio = (w / h) / 2.0;
				Gl.glOrtho(-ratio, ratio, -0.5, 0.5, -100, 100);
			}
			else
			{
				double ratio = (h / w) / 2.0;
				Gl.glOrtho(-0.5, 0.5, -ratio, ratio, -100, 100);
			}

			Gl.glMatrixMode(Gl.GL_MODELVIEW);
			Gl.glLoadIdentity();
		}
		private void InitMatrix()
		{
			if (this.useGLOrthoProjection)
			{
				InitMatrix_Ortho();
				return;
			}

			double w = Size.Width;
			double h = Size.Height;
			Gl.glViewport(0, 0, Size.Width, Size.Height);
			Gl.glMatrixMode(Gl.GL_PROJECTION);
			Gl.glLoadIdentity();
			
			Glu.gluPerspective(60.0, w / h, 0.01, 100);

			Gl.glMatrixMode(Gl.GL_MODELVIEW);
			Gl.glLoadIdentity();

			Glu.gluLookAt(0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
		}


		private void DrawWireframe(Mesh m, Color c)
		{
			Gl.glPolygonMode(Gl.GL_FRONT_AND_BACK, Gl.GL_LINE);
			Gl.glDisable(Gl.GL_CULL_FACE);
			Gl.glDisable(Gl.GL_LIGHTING);
			Gl.glEnable(Gl.GL_LINE_SMOOTH);
			Gl.glColor3ub(c.R, c.G, c.B);

			Gl.glLineWidth(1.0f);
			Gl.glEnableClientState(Gl.GL_VERTEX_ARRAY);
			fixed (double* vp = m.VertexPos)
			fixed (int* index = m.FaceIndex)
			{
				Gl.glVertexPointer(3, Gl.GL_DOUBLE, 0, new IntPtr(vp));
				Gl.glDrawElements(Gl.GL_TRIANGLES, m.FaceCount * 3, Gl.GL_UNSIGNED_INT, new IntPtr(index));
			}
			Gl.glDisableClientState(Gl.GL_VERTEX_ARRAY);
			Gl.glEnable(Gl.GL_CULL_FACE);
			Gl.glEnable(Gl.GL_LIGHTING);
		}
		private void DrawFlatShaded(Mesh m, Color c)
		{
			Gl.glShadeModel(Gl.GL_FLAT);
			Gl.glPolygonMode(Gl.GL_FRONT_AND_BACK, Gl.GL_FILL);
			Gl.glEnable(Gl.GL_LIGHTING);
			Gl.glEnable(Gl.GL_NORMALIZE);

			Gl.glColor3ub(c.R, c.G, c.B);
			fixed (double* vp = m.VertexPos)
			fixed (double* np = m.FaceNormal)
			fixed (int* fi = m.FaceIndex)
			{
				Gl.glBegin(Gl.GL_TRIANGLES);
				for (int i = 0, j = 0; i < m.FaceCount; i++, j += 3)
				{
					Gl.glNormal3dv(new IntPtr(np + j));
					Gl.glVertex3dv(new IntPtr(vp + fi[j] * 3));
					Gl.glVertex3dv(new IntPtr(vp + fi[j + 1] * 3));
					Gl.glVertex3dv(new IntPtr(vp + fi[j + 2] * 3));
				}
				Gl.glEnd();
			}
			Gl.glDisable(Gl.GL_LIGHTING);
		}
		private void DrawSmoothShaded(Mesh m)
		{
			Gl.glShadeModel(Gl.GL_SMOOTH);
			Gl.glPolygonMode(Gl.GL_FRONT_AND_BACK, Gl.GL_FILL);
			Gl.glEnable(Gl.GL_LIGHTING);
			Gl.glEnable(Gl.GL_NORMALIZE);

			//Color c = this.meshColor;
			Color c = Color.Pink;
			Gl.glColor3ub(c.R, c.G, c.B);
			Gl.glEnableClientState(Gl.GL_VERTEX_ARRAY);
			Gl.glEnableClientState(Gl.GL_NORMAL_ARRAY);
			fixed (double* vp = m.VertexPos)
			fixed (double* np = m.VertexNormal)
			fixed (int* index = m.FaceIndex)
			{
				Gl.glVertexPointer(3, Gl.GL_DOUBLE, 0, new IntPtr(vp));
				Gl.glNormalPointer(Gl.GL_DOUBLE, 0, new IntPtr(np));
				Gl.glDrawElements(Gl.GL_TRIANGLES, m.FaceCount * 3, Gl.GL_UNSIGNED_INT, new IntPtr(index));
			}
			Gl.glDisableClientState(Gl.GL_VERTEX_ARRAY);
			Gl.glDisableClientState(Gl.GL_NORMAL_ARRAY);
			Gl.glDisable(Gl.GL_LIGHTING);
		}
		private void DrawTransparent(Mesh m, Color c)
		{
			byte alphaChannel = 60;

            Gl.glPushClientAttrib(Gl.GL_CLIENT_VERTEX_ARRAY_BIT);
            Gl.glPushAttrib(Gl.GL_ALL_ATTRIB_BITS);

            Gl.glShadeModel(Gl.GL_SMOOTH);
            Gl.glEnable(Gl.GL_NORMALIZE);
            Gl.glDisable(Gl.GL_LIGHTING);
            Gl.glEnable(Gl.GL_BLEND);
            Gl.glBlendFunc(Gl.GL_SRC_ALPHA, Gl.GL_ONE_MINUS_SRC_ALPHA);
            Gl.glEnable(Gl.GL_CULL_FACE);
            Gl.glDisable(Gl.GL_DEPTH_TEST);
            Gl.glPolygonMode(Gl.GL_FRONT_AND_BACK, Gl.GL_FILL);

            Gl.glColor4ub(c.R, c.G, c.B, alphaChannel);

            Gl.glEnableClientState(Gl.GL_VERTEX_ARRAY);
            Gl.glEnableClientState(Gl.GL_NORMAL_ARRAY);
            fixed (double* vp = m.VertexPos)
            fixed (double* np = m.VertexNormal)
            fixed (int* index = m.FaceIndex)
            {
              Gl.glVertexPointer(3, Gl.GL_DOUBLE, 0, (IntPtr)vp);
              Gl.glNormalPointer(Gl.GL_DOUBLE, 0, (IntPtr)np);
              Gl.glDrawElements(Gl.GL_TRIANGLES, m.FaceCount * 3, Gl.GL_UNSIGNED_INT, (IntPtr)index);
            }
            Gl.glDisableClientState(Gl.GL_VERTEX_ARRAY);
            Gl.glDisableClientState(Gl.GL_NORMAL_ARRAY);

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

			Gl.glEnable(Gl.GL_LIGHTING);
			Gl.glDepthFunc(Gl.GL_EQUAL);

			float[] LightDiffuse = { 1.0f, 1.0f, 1.0f, 1.0f };
			float[] LightSpecular = { 1.0f, 1.0f, 1.0f, 1.0f };
			Gl.glLightfv(Gl.GL_LIGHT0, Gl.GL_DIFFUSE, LightDiffuse);
			Gl.glLightfv(Gl.GL_LIGHT0, Gl.GL_SPECULAR, LightSpecular);
			Gl.glEnable(Gl.GL_LIGHT0);

			float[] Ambient = { 0.0f, 0.0f, 0.0f, 1.0f };
			float[] Diffuse = { 0.0f, 0.0f, 0.0f, 1.0f };
			float[] SpecularRef = { 1.0f, 1.0f, 1.0f, 1.0f };
			
			Gl.glColorMaterial(Gl.GL_FRONT, Gl.GL_DIFFUSE);
			Gl.glMaterialfv(Gl.GL_FRONT, Gl.GL_DIFFUSE, Diffuse);
			Gl.glMaterialfv(Gl.GL_FRONT, Gl.GL_AMBIENT, Ambient);
			Gl.glMaterialfv(Gl.GL_FRONT, Gl.GL_SPECULAR, SpecularRef);
			Gl.glMateriali(Gl.GL_FRONT, Gl.GL_SHININESS, 128);

			Gl.glColor4ub(c.R, c.G, c.B, alphaChannel);

			Gl.glEnableClientState(Gl.GL_VERTEX_ARRAY);
			Gl.glEnableClientState(Gl.GL_NORMAL_ARRAY);
			fixed (double* vp = m.VertexPos)
			fixed (double* np = m.VertexNormal)
			fixed (int* index = m.FaceIndex)
			{
				Gl.glVertexPointer(3, Gl.GL_DOUBLE, 0, (IntPtr)vp);
				Gl.glNormalPointer(Gl.GL_DOUBLE, 0, (IntPtr)np);
				Gl.glDrawElements(Gl.GL_TRIANGLES, m.FaceCount * 3, Gl.GL_UNSIGNED_INT, (IntPtr)index);
			}
			Gl.glDisableClientState(Gl.GL_VERTEX_ARRAY);
			Gl.glDisableClientState(Gl.GL_NORMAL_ARRAY);


			Gl.glDisable(Gl.GL_NORMALIZE);
            Gl.glDisable(Gl.GL_LIGHTING);
            Gl.glDisable(Gl.GL_BLEND);
            Gl.glDisable(Gl.GL_CULL_FACE);
            Gl.glEnable(Gl.GL_DEPTH_TEST);

			Gl.glPopAttrib();
			Gl.glPopClientAttrib();

		}
		private void DrawCuboid3D(PolyProxy cuboid, Color c, float width)
		{
			if (cuboid == null) return;

			Gl.glEnable(Gl.GL_BLEND);
			Gl.glDisable(Gl.GL_CULL_FACE);
			Gl.glBlendFunc(Gl.GL_SRC_ALPHA, Gl.GL_ONE_MINUS_SRC_ALPHA);
			Color temp = Color.FromArgb(30, c);
			Gl.glColor4ub(temp.R, temp.G, temp.B, temp.A);
			Gl.glEnable(Gl.GL_CULL_FACE);


			Gl.glDisable(Gl.GL_LIGHTING);
			Gl.glColor3ub(c.R, c.G, c.B);
			Gl.glEnable(Gl.GL_LINE_SMOOTH);
			Gl.glLineWidth(width);
			Gl.glBegin(Gl.GL_LINES);
			foreach (Polygon poly in cuboid.polygons)
			{
				for (int i = 0; i < poly.points3.Count; ++i)
				{
					Point3 p = poly.points3[i];
					Point3 q = poly.points3[(i + 1) % poly.points3.Count];
					Gl.glVertex3d(p.pos.x, p.pos.y, p.pos.z);
					Gl.glVertex3d(q.pos.x, q.pos.y, q.pos.z);
				}
			}
			Gl.glEnd();
			Gl.glLineWidth(1.0f);
			Gl.glEnable(Gl.GL_LIGHTING);

		}
		

		public void Clear()
		{
			this.synthesizedCuboidMeshComponents.Clear();
			this.synthesisTargetsFiles.Clear();
			this._curr_mesh_displayindex = 0;
			this.Refresh();
		}

	
		public List<List<CuboidMeshComponent>> synthesizedCuboidMeshComponents = new List<List<CuboidMeshComponent>>();
		public int _curr_mesh_displayindex = 0;
		public int _nmeshes_to_display = 4;
		private List<string> synthesisTargetsFiles = new List<string>();
		public bool HideBox = false;


		public void NextModel()
		{
			if (this.synthesizedCuboidMeshComponents.Count > 0)
				_curr_mesh_displayindex = (_curr_mesh_displayindex + 1) % this.synthesizedCuboidMeshComponents.Count;
			this.Refresh();
		}
		public void PrevModel()
		{
			if (this.synthesizedCuboidMeshComponents.Count > 0)
				_curr_mesh_displayindex = (_curr_mesh_displayindex - 1 + this.synthesizedCuboidMeshComponents.Count) % this.synthesizedCuboidMeshComponents.Count;
			this.Refresh();
		}


		// graph based object synthesis
		public bool IsComponentsSlotsCompatible(ComponentGraphNode node1, ComponentGraphNode node2, Slot ref_slot)
		{
			List<Slot> slots1 = new List<Slot>();
			foreach (Slot s in node1.contactSlots)
			{
				if (s.symmetrySupport == ref_slot.symmetrySupport)
				{
					slots1.Add(s);
				}
			}
			List<Slot> slots2 = new List<Slot>();
			foreach (Slot s in node2.contactSlots)
			{
				if (s.symmetrySupport == ref_slot.symmetrySupport)
				{
					slots2.Add(s);
				}
			}
			if (slots1.Count != slots2.Count) return false;
			for (int i = 0; i < slots1.Count; ++i)
			{
				int k = slots1[i].contactPoints.Count;
				for (int j = 0; j < slots2.Count; ++j)
				{
					if (slots2[j].contactPoints.Count == k)
					{
						slots2.RemoveAt(j);
					}
				}
			}
			return slots2.Count == 0;
		}
		public void GetCompatibleSlots(ComponentGraphNode node1, ComponentGraphNode node2, out List<Slot> sourceSlots, out List<Slot> targetSlots)
		{
			sourceSlots = new List<Slot>();
			targetSlots = new List<Slot>();
			foreach (ComponentGraphEdge e in node1.edges)
			{
				if (e.type == ComponentEdgeType.Symmetry) continue;
				
				ComponentGraphNode other = e.start == node1 ? e.end : e.start;
				if (other.replace != null) // already replaced
				{
					Slot targetSlot = null, sourceSlot = null;

					// find the corresponding slot
					ComponentGraphNode prevnode = other, prevreplacednode = other.replace;

					// find the slot connect to #other#, the already replaced node
					Slot slot = null;
					foreach (Slot slt in node1.contactSlots)
					{
					    if (slt.start == prevnode || slt.end == prevnode)
					    {
					        slot = slt;
					        break;
					    }
					}
					if (slot == null) throw new Exception();
					bool targetCompatible = this.IsComponentsSlotsCompatible(prevnode, prevreplacednode, slot);
					bool sourceCompatible = this.IsComponentsSlotsCompatible(node1, node2, slot);
					if (!targetCompatible) // need create slots
					{
			//			Program.OutputText("creating target slots ", true);

						targetSlot = this.CreateCorrespondingSlot(prevnode, prevreplacednode, slot, false);
						targetSlot.symmetry = this.CreateCorrespondingSlot(prevnode, prevreplacednode, slot.symmetry, false);
					}
					else
					{
						//// find the target slot in #other#'s replaced node (which the current node will connected to)
						double mindis = double.MaxValue;
						foreach (Slot slt in prevreplacednode.contactSlots)
						{
						    if (slt.used) continue;

						    if (!IsSlotCompatible(slt, slot)) continue;
						    double d = (slt.oldcenter - slot.oldcenter).Length();
						    if (d < mindis)
						    {
						        targetSlot = slt;
						        mindis = d;
						    }
						}
					}

					if (targetSlot != null)
					{
						if (!sourceCompatible)
						{
			//				Program.OutputText("creating source slots ", true);
							Slot ss = targetCompatible ? slot : targetSlot;
							sourceSlot = this.CreateCorrespondingSlot(node1, node2, ss, false);
							sourceSlot.symmetry = this.CreateCorrespondingSlot(node1, node2, ss.symmetry, false);
						}
						else
						{
							Vector3d t12 = node1.pos - node2.pos;
							// find the slot of node2 where it should connect to the targetSlot
							double mindis = double.MaxValue;
							foreach (Slot slt in node2.contactSlots)
							{
								if (!IsSlotCompatible(slt, slot)) continue;
								double dis = (slt.oldcenter - targetSlot.oldcenter + t12).Length();
								if (dis < mindis)
								{
									sourceSlot = slt;
									mindis = dis;
								}
							}
						}
					}

					if (sourceSlot == null || targetSlot == null) throw new Exception();

					sourceSlots.Add(sourceSlot);
					targetSlots.Add(targetSlot);

				}
			}
		}
		public void ReplaceNode(ComponentGraphNode node1, ComponentGraphNode node2, out List<CuboidMeshComponent> cmc, out List<CuboidMeshComponent> symmcmc)
		{
			cmc = this.CreateCuboidComponent(node2);
			symmcmc = null;
			if (node2.symmetry != null)
			{
				symmcmc = this.CreateCuboidComponent(node2.symmetry);
			}

			// replace node1 with node2, and create the corresponding cuboid component
			// 1. find contact and target slots in adjacent replaced nodes
			List<Slot> sourceSlots = null;
			List<Slot> targetSlots = null;
			
			this.GetCompatibleSlots(node1, node2, out sourceSlots, out targetSlots);

			foreach (Slot slot in sourceSlots) slot.used = true;
			foreach (Slot slot in targetSlots) slot.used = true;
			
			int n = sourceSlots.Count;
			List<Vector3d> srcpts = new List<Vector3d>();
			List<Vector3d> tarpts = new List<Vector3d>();
		//	double s = Math.Max(Math.Max(node1.scale.x, node1.scale.y), node1.scale.z) / Math.Max(Math.Max(node2.scale.x, node2.scale.y), node2.scale.z);
			double s = (node1.scale.x + node1.scale.y + node1.scale.z) / (node2.scale.x+ node2.scale.y + node2.scale.z);
			if (n == 0)
			{
				srcpts.Add(node2.pos);
				tarpts.Add(node1.pos);
			}
			else
			{
				bool usingICP = true;
				for (int i = 0; i < n; ++i)
				{
					if (sourceSlots[i].contactPoints.Count == 1 || sourceSlots[i].contactPoints.Count != targetSlots[i].contactPoints.Count)
					{
						srcpts.Add(sourceSlots[i].oldcenter);
						tarpts.Add(targetSlots[i].center);
					}
					else if (usingICP)
					{
						List<Vector3d> temppts = new List<Vector3d>();
						Vector3d tt = sourceSlots[i].oldcenter - targetSlots[i].center;
						foreach (Vector3d pt in sourceSlots[i].oldContactPoints)
						{
							srcpts.Add(pt);
							double mindis = double.MaxValue; Vector3d qq = new Vector3d();
							foreach (Vector3d qt in targetSlots[i].contactPoints)
							{
								if (temppts.Contains(qt)) continue;
								double d = (qt - pt + tt).Length();
								if (d < mindis)
								{
									mindis = d;
									qq = qt;
								}
							}
							tarpts.Add(qq);
							temppts.Add(qq);
						}
					}
					else
					{
						srcpts.AddRange(sourceSlots[i].contactPoints);
						tarpts.AddRange(targetSlots[i].contactPoints);
					}
				}
			}


			Vector3d scale = new Vector3d(1.0, 1.0, 1.0);
			Vector3d trans = new Vector3d();
			Matrix4d Q = null, S = null, T = null;
			n = srcpts.Count;
			if (n < 2)
			{
				if (!node1.isgroundcontact)
				{
					scale = new Vector3d(s, s, s);
					trans = tarpts[0] - srcpts[0];
					S = Matrix4d.ScalingMatrix(scale.x, scale.y, scale.z);
					Q = Matrix4d.TranslationMatrix(tarpts[0]) * S * Matrix4d.TranslationMatrix(new Vector3d() - srcpts[0]);
				}
				else
				{
					srcpts.Add(new Vector3d(srcpts[0].x, 0, srcpts[0].z));
					tarpts.Add(new Vector3d(tarpts[0].x, 0, tarpts[0].z));
					
					this.FindBestTransform(srcpts, tarpts, out S, out T, out Q);
				}
			}
			else if (n == 2 && !node1.isgroundcontact) // find translation, scaling, and rotation
			{
				this.FindBestTransform(srcpts, tarpts, out S, out T, out Q);
			}
			else if (n == 3)
			{
			    Vector3d t1 = new Vector3d();
			    Vector3d t2 = new Vector3d();
			    foreach (Vector3d tt in srcpts)
			        t1 += tt;
			    foreach (Vector3d tt in tarpts)
			        t2 += tt;
			    t1 /= srcpts.Count;
			    t2 /= tarpts.Count;

			    trans = t2 - t1;

			    double mindis = double.MaxValue; int jj = -1;
			    for (int i = 0; i < 3; ++i)
			    {
			        Vector3d pt = tarpts[i];
			        double dis = (pt - t1).Length();
			        if (dis < mindis)
			        {
			            mindis = dis;
			            jj = i;
			        }
			    }

				Vector3d u1 = tarpts[jj] - tarpts[(jj + 2) % 3];
				Vector3d u2 = tarpts[jj] - tarpts[(jj + 1) % 3];
				double cos = u1.Normalize().Dot(u2.Normalize());
				double cos170 =  Math.Cos(17 / 18.0 * Math.PI);
				if (cos < cos170)
				{
					srcpts.RemoveAt(jj);
					tarpts.RemoveAt(jj);
				}
					
				this.FindBestTransform(srcpts, tarpts, out S, out T, out Q);

			}
			else // least-squaers solve
			{
				if (n == 2 && node2.isgroundcontact)
				{
					srcpts.Add(new Vector3d(srcpts[0].x, 0, srcpts[0].z));
					srcpts.Add(new Vector3d(srcpts[1].x, 0, srcpts[1].z));
					tarpts.Add(new Vector3d(tarpts[0].x, 0, tarpts[0].z));
					tarpts.Add(new Vector3d(tarpts[1].x, 0, tarpts[1].z));
					n = srcpts.Count;
				}

				this.FindBestTransform(srcpts,tarpts, out S, out T, out Q);

			//	Q = T * Matrix4d.TranslationMatrix(node2.pos) * S * Matrix4d.TranslationMatrix(new Vector3d() - node2.pos);

			}

			foreach (CuboidMeshComponent mc in cmc)
				mc.Transform(Q);
			if (cmc.Count > 1)
				node2.TransformSlots(Q);

			// symmetry
			if (symmcmc != null)
			{
				// get bounding box center
				Vector3d c1 = GetBoundingBoxCenterAndScale(cmc)[0];

			    Vector3d normal = (node1.pos - node1.symmetry.pos).Normalize();
			    Vector3d center = (node1.pos + node1.symmetry.pos)/2;
			    Vector3d cntr2 = HelperFunction.GetMirrorSymmetryPoint(c1, normal, center);
				
			    Vector3d cntr = GetBoundingBoxCenterAndScale(symmcmc)[0];
			    Vector3d tt = cntr2;
			    Q = Matrix4d.TranslationMatrix(tt) * S * Matrix4d.TranslationMatrix(new Vector3d() - cntr);
			    
				foreach (CuboidMeshComponent mc in symmcmc)
					mc.Transform(Q);
				if (symmcmc.Count > 1)
					node2.symmetry.TransformSlots(Q);
			   this.SetNodeReplace(node1.symmetry, node2.symmetry);
			}


			// assign and set tags
			this.SetNodeReplace(node1,node2);

			this.SetNodeTag(node1, true);
			this.SetNodeTag(node2, true);

			if (symmcmc != null)
			{
				this.SetNodeTag(node1.symmetry, true);
				this.SetNodeTag(node2.symmetry, true);
				
			}

		//	Program.OutputText("-------------------------------------------", true);

		}
		public void FindBestTransform(List<Vector3d> srcpts, List<Vector3d> tarpts, out Matrix4d S, out Matrix4d T, out Matrix4d Q)
		{
			int n = srcpts.Count;
			if (n == 2)
			{
				Vector3d c1 = (srcpts[0] + srcpts[1]) / 2;
				Vector3d c2 = (tarpts[0] + tarpts[1]) / 2;
				Vector3d v1 = srcpts[1] - srcpts[0];
				Vector3d v2 = tarpts[1] - tarpts[0];
				if (v1.Dot(v2) < 0) v1 = new Vector3d() - v1;
				double ss = v2.Length() / v1.Length();
				S = Matrix4d.ScalingMatrix(ss, ss, ss);
				Matrix4d R = Matrix4d.IdentityMatrix();
				double cos = v1.Normalize().Dot(v2.Normalize());
				if (cos < Math.Cos(1.0 / 18 * Math.PI))
				{
					Vector3d axis = v1.Cross(v2).Normalize();
					double theta = Math.Acos(cos);
					R = Matrix4d.RotationMatrix(axis, theta);
				}
				T = Matrix4d.TranslationMatrix(c2-c1);
				Q = Matrix4d.TranslationMatrix(c2) * R * S * Matrix4d.TranslationMatrix(new Vector3d() - c1);
			}
			else
			{
				Vector3d t1 = new Vector3d();
				Vector3d t2 = new Vector3d();
				foreach (Vector3d tt in srcpts)
					t1 += tt;
				foreach (Vector3d tt in tarpts)
					t2 += tt;
				t1 /= srcpts.Count;
				t2 /= tarpts.Count;

				Vector3d trans = t2 - t1;
				T = Matrix4d.TranslationMatrix(trans);

				// find the scales
				int k = srcpts.Count;
				double sx = 0, sy = 0, sz = 0;
				for (int i = 0; i < k; ++i)
				{
					Vector3d p1 = srcpts[i] - t1;
					Vector3d p2 = tarpts[i] - t2;
					sx += p2.x / p1.x;
					sy += p2.y / p1.y;
					sz += p2.z / p1.z;
				}
				sx /= k;
				sy /= k;
				sz /= k;

				if (double.IsNaN(sx) || double.IsInfinity(sx))
				{
					sx = 1.0;
				}
				if (double.IsNaN(sy) || double.IsInfinity(sy))
				{
					sy = 1.0;
				}
				if (double.IsNaN(sz) || double.IsInfinity(sz))
				{
					sz = 1.0;
				}

				// adjust scale 
				if (n > 3)
				{
					//// find the points plane
					//double[] points = new double[srcpts.Count * 3];
					//for (int i = 0, jj = 0; i < srcpts.Count; ++i, jj += 3)
					//{
					//    points[jj] = srcpts[i].x;
					//    points[jj + 1] = srcpts[i].y;
					//    points[jj + 2] = srcpts[i].z;
					//}
					//double[] plane = new double[4];
					//PlaneFitter.ZyyPlaneFitter _plfitter = new PlaneFitter.ZyyPlaneFitter();
					//fixed (double* _pts = points, _pl = plane)
					//    _plfitter.GetFittingPlane(_pts, srcpts.Count, _pl);
					//Vector3d normal = new Vector3d(plane, 0);
					//normal = normal.Normalize();

					// 4 permutations, 012, 123, 023, 013
					Vector3d aa = tarpts[0];
					Vector3d bb = tarpts[1];
					Vector3d cc = tarpts[2];
					Vector3d dd = tarpts[3];
					Vector3d[] nn = new Vector3d[4] {
						((aa - bb).Cross(bb - cc)).Normalize(),
						((bb - cc).Cross(cc - dd)).Normalize(),
						((aa - cc).Cross(cc - dd)).Normalize(),
						((aa - bb).Cross(bb - dd)).Normalize()
					};
					Vector3d nor = new Vector3d();
					for (int i = 0; i < 4; ++i)
					{
						if (!double.IsNaN(nn[i].x))
						{
							nor = nn[i];
							break;
						}
					}
					Vector3d normal = new Vector3d();
					int count = 0;
					for (int i = 0; i < 4; ++i)
					{
						if (!double.IsNaN(nn[i].x))
						{
							if (nn[i].Dot(nor) < 0)
							{
								nn[i] = new Vector3d() - nn[i];
							}
							normal += nn[i];
							count++;
						}
					}
					normal = normal.Normalize();

					if (double.IsNaN(normal.x)) throw new Exception();
					{
						if (Math.Abs(normal.x) > 0.5)
						{
							sx = 1.0;
						}
						else if (Math.Abs(normal.y) > 0.5)
						{
							sy = 1.0;
						}
						else if (Math.Abs(normal.z) > 0.5)
							sz = 1.0;
					}
				}

				Vector3d scale = new Vector3d(sx, sy, sz);

				if (double.IsNaN(scale.x) || double.IsNaN(trans.x)) throw new Exception();

				//Program.OutputText("trans: " + trans.x + " " + trans.y + " " + trans.z, true);
				//Program.OutputText("scale: " + scale.x + " " + scale.y + " " + scale.z, true);

				//// check error
				S = Matrix4d.ScalingMatrix(scale.x, scale.y, scale.z);
				Q = Matrix4d.TranslationMatrix(t2) * S * Matrix4d.TranslationMatrix(new Vector3d() - t1);
			}
		}
		public static Vector3d[] GetBoundingBoxCenterAndScale(List<CuboidMeshComponent> cmcs)
		{
			Vector3d maxCoord = new Vector3d(double.MinValue, double.MinValue, double.MinValue);
			Vector3d minCoord = new Vector3d(double.MaxValue, double.MaxValue, double.MaxValue);
			foreach (CuboidMeshComponent cm in cmcs)
			{
				foreach (Point3 pt in cm.boundingBox.points3)
				{
					Vector3d v = pt.pos;
					maxCoord = Vector3d.Max(maxCoord, v);
					minCoord = Vector3d.Min(minCoord, v);
				}
			}
			return new Vector3d[2] {(maxCoord + minCoord)/2, new Vector3d(
				Math.Abs(maxCoord.x - minCoord.x),
				Math.Abs(maxCoord.y - minCoord.y),
				Math.Abs(maxCoord.z - minCoord.z))
			};
		}
		public bool IsSelfCloning(ComponentGraph g)
		{
			foreach (ComponentGraphNode node in g.nodes)
			{
				if (node.replace != node)
					return false;
			}
			return true;
		}
		public bool IsCompleteReplacing(ComponentGraph g)
		{
			foreach (ComponentGraphNode node in g.nodes)
			{
				if (node.replace != null && node.replace.name == node.name)
					return false;
			}
			return true;
		}
		public Slot CreateCorrespondingSlot(ComponentGraphNode originalnode, ComponentGraphNode replacednode, Slot slot, bool usingICP)
		{
			if (slot == null) return null;

			// create in replacenode the corresponding node in the original node, slot is \in original node
			// 1. find the scale in two nodes
			Vector3d s1 = originalnode.scale, s2 = replacednode.scale;
			Vector3d t1 = originalnode.pos, t2 = replacednode.pos;
			if (originalnode.cuboidMeshComponent != null)
			{
				s1 = originalnode.cuboidMeshComponent.GetScale();
				t1 = originalnode.cuboidMeshComponent.boundingBox.center.pos;
			}
			if (replacednode.cuboidMeshComponent != null)
			{
				s2 = replacednode.cuboidMeshComponent.GetScale();
				t2 = replacednode.cuboidMeshComponent.boundingBox.center.pos;
			}

			Vector3d ss = new Vector3d(s2.x / s1.x, s2.y / s1.y, s2.z / s1.z);
			Matrix4d Q = Matrix4d.TranslationMatrix(t2) * Matrix4d.ScalingMatrix(ss.x, ss.y, ss.z) * Matrix4d.TranslationMatrix(new Vector3d()-t1);
			List<Vector3d> slotpoints = new List<Vector3d>();
			foreach (Vector3d t in slot.contactPoints)
			{
				Vector3d pt = (Q * new Vector4d(t, 1)).XYZ();
				slotpoints.Add(pt);
			}

			// if target cuboid components exist, find the cloest mesh point
			if (replacednode.cuboidMeshComponent != null && usingICP)
			{
				Mesh mesh = replacednode.cuboidMeshComponent.mesh;
				int N = mesh.VertexCount;
				List<Vector3d> newpoints = new List<Vector3d>();
				foreach (Vector3d pt in slotpoints)
				{
					double mindis = double.MaxValue; Vector3d target = pt;
					for (int i = 0, j = 0; i < N; ++i, j += 3)
					{
						Vector3d v = new Vector3d(mesh.VertexPos, j);
						double dis = (pt - v).Length();
						if (dis < mindis)
						{
							mindis = dis;
							target = v;
						}
					}
					for (int i = 0,j = 0; i < mesh.FaceCount; ++i,j+=3)
					{
					    int c0 = mesh.FaceIndex[j]*3;
					    int c1 = mesh.FaceIndex[j+1]*3;
					    int c2 = mesh.FaceIndex[j+2]*3;
					    Vector3d v0 = new Vector3d(mesh.VertexPos, c0);
					    Vector3d v1 = new Vector3d(mesh.VertexPos, c1);
					    Vector3d v2 = new Vector3d(mesh.VertexPos, c2);
					    Vector3d nn = new Vector3d(mesh.FaceNormal, j);
					    Vector3d pc = pt - v0;
					    Vector3d cc = pt - (pc).Dot(nn)*nn;
					    if (PointInTriangle(cc, v0, v1, v2))
					    {
					        double dis = Math.Abs(pc.Dot(nn));
					        if (dis < mindis)
					        {
					            mindis = dis;
					            target = cc;
					        }
					    }
					}
					newpoints.Add(target);
				}
				slotpoints = newpoints;
			}

			Slot another = new Slot(slot.start, slot.end, slotpoints, slot.index);
			return another;
		}
		public bool IsTopologicalCompatible(ComponentGraphNode node1, ComponentGraphNode node2)
		{
			if (node1 == node2) return true;
			if (node1.isgroundcontact != node2.isgroundcontact) return false;

			// edge counts for fuzzy correspondence
			int underneathEdgesCount1 = 0, uponEdgesCount1 = 0, symmetryEdgesCount1 = 0, contactEdgeCount1 = 0;
			int underneathEdgesCount2 = 0, uponEdgesCount2 = 0, symmetryEdgesCount2 = 0, contactEdgeCount2 = 0;
			foreach (ComponentGraphEdge e in node1.edges)
			{
				switch (e.type)
				{
					case ComponentEdgeType.Underneath:
						underneathEdgesCount1++;
						break;
					case ComponentEdgeType.Upon:
						uponEdgesCount1++;
						break;
					case ComponentEdgeType.Symmetry:
						symmetryEdgesCount1++;
						break;
					case ComponentEdgeType.Contact:
						contactEdgeCount1++;
						break;
				}
			}
			foreach (ComponentGraphEdge e in node2.edges)
			{
				switch (e.type)
				{
					case ComponentEdgeType.Underneath:
						underneathEdgesCount2++;
						break;
					case ComponentEdgeType.Upon:
						uponEdgesCount2++;
						break;
					case ComponentEdgeType.Symmetry:
						symmetryEdgesCount2++;
						break;
					case ComponentEdgeType.Contact:
						contactEdgeCount2++;
						break;
				}
			}

			return ((underneathEdgesCount1 == underneathEdgesCount2 && underneathEdgesCount1 > 0) || 
					(uponEdgesCount1 == uponEdgesCount2 && uponEdgesCount1 > 0) || 
					((underneathEdgesCount1 == underneathEdgesCount2 || uponEdgesCount1 == uponEdgesCount2) && 
						contactEdgeCount1 == contactEdgeCount2 && contactEdgeCount1 > 0)
					)
				&& symmetryEdgesCount1 == symmetryEdgesCount2;
		}
		public bool IsSlotCompatible(Slot s, Slot t)
		{
			return s.symmetrySupport == t.symmetrySupport;
			if (s.symmetrySupport == t.symmetrySupport && s.symmetrySupport != Functionality.none) return true;
			return (IsTopologicalCompatible(s.start, t.start) && IsTopologicalCompatible(s.end, t.end)) ||
				(IsTopologicalCompatible(s.end, t.start) && IsTopologicalCompatible(s.start, t.end));
		}
		public void ResetNodes()
		{
			// reset all tags
			foreach (ComponentGraph cg in spatialgraphs)
			{
				foreach (ComponentGraphNode nd in cg.nodes)
				{
					nd.tag = false;
					nd.replace = null;
					nd.cuboidMeshComponent = null;
					nd.compatibleNodes.Clear();
					nd.parent = null;
				}
			}
		}
		public void ResetNodesTags()
		{
			// reset all tags
			foreach (ComponentGraph cg in spatialgraphs)
			{
				foreach (ComponentGraphNode nd in cg.nodes)
				{
					nd.tag = false;
				}
			}
		}
		public void ResetNodesParent()
		{
			// reset all compatible nodes
			foreach (ComponentGraph cg in spatialgraphs)
			{
				foreach (ComponentGraphNode nd in cg.nodes)
				{
					nd.parent = null;
				}
			}
		}
		public void ResetNodesSlots()
		{
			// reset all compatible nodes
			foreach (ComponentGraph cg in spatialgraphs)
			{
				foreach (ComponentGraphNode nd in cg.nodes)
				{
					foreach (Slot s in nd.contactSlots)
						s.Reset();
				}
			}
		}


		// continued
		public bool IsNondeContained(ComponentGraphNode node, List<ComponentGraphNode> nodes)
		{
			foreach (ComponentGraphNode nd in nodes)
			{
				if (nd == node) return true;
				if (nd.kids != null && nd.kids.Contains(node)) return true;
			}
			return false;
		}
		public List<CuboidMeshComponent> CreateCuboidComponent(ComponentGraphNode node)
		{
			List<CuboidMeshComponent> ls = new List<CuboidMeshComponent>();
			if (node.kids == null) 
				ls.Add(new CuboidMeshComponent(node, node.meshfile));
			else
			foreach (ComponentGraphNode nd in node.kids)
			{
				ls.Add(new CuboidMeshComponent(nd, nd.meshfile));
			}
			return ls;
		}
		public void SetNodeTag(ComponentGraphNode node, bool tag)
		{
			node.tag = tag;
			if (node.kids != null)
			{
				foreach (ComponentGraphNode nd in node.kids)
					nd.tag = tag;
			}
		}
		public void SetNodeReplace(ComponentGraphNode node, ComponentGraphNode node2)
		{
			node.replace = node2;
			if (node.kids != null)
			{
				if (node == node2) // self replacing
				{
					foreach (ComponentGraphNode nd in node.kids)
						nd.replace = nd;
				}
				else
				foreach (ComponentGraphNode nd in node.kids)
					nd.replace = node2;
			}
		}
		public void TripletsGraphReshuffle(ComponentGraph g, List<ComponentGraph> targetgraphs)
		{
			if (g == null || targetgraphs.Count == 0) return;


			string targetfile = targetgraphs[0].meshfile;
			int idx = targetfile.LastIndexOf('\\');
			targetfile = targetfile.Substring(idx + 1, targetfile.Length - 5 - idx);


			// gets symmetry triplets...
			if (g.symmtryTriplets == null) g.GetSymmetryTriplets();
			foreach (ComponentGraph gg in targetgraphs)
			{
				if (gg.symmtryTriplets == null)
					gg.GetSymmetryTriplets();
			}

			// find compatible triplets for each triplet in g
			foreach (SymmetryTriplet tri in g.symmtryTriplets)
			{
				tri.compatibleTriplets.Clear();
				foreach (ComponentGraph gg in targetgraphs)
				{
					this.FindSymmetryNodesCompatible(g, gg);
					foreach (SymmetryTriplet tri2 in gg.symmtryTriplets)
					{
						if (IsCompatible(tri, tri2))
						{
							tri.compatibleTriplets.Add(tri2);
						}
					}
				}
			}

			// replace triplets in g with triplets in others, do we need an order?? -- ground touching first..
			List<SymmetryTriplet> first = new List<SymmetryTriplet>();
			List<SymmetryTriplet> second = new List<SymmetryTriplet>();
			foreach (SymmetryTriplet tri in g.symmtryTriplets)
			{
				if (tri.compatibleTriplets.Count == 0) continue;
				if (tri.isgroundtouching)
				{
					first.Add(tri);
				}
				else
					second.Add(tri);
			}
			first.AddRange(second);


			// handle super triplets
			List<SymmetryTriplet> _tris_remove = new List<SymmetryTriplet>();
			foreach (SymmetryTriplet tri in first)
			{
				if (tri.kidsTriplets != null)
				{
					bool _all_compatible = true;
					foreach (SymmetryTriplet tr in tri.kidsTriplets)
					{
						if (tr.compatibleTriplets.Count == 0)
						{
							_all_compatible = false;
							break;
						}
					}
					if (_all_compatible)
					{
						_tris_remove.Add(tri);
					}
				}
			}
			foreach (SymmetryTriplet tri in _tris_remove)
				first.Remove(tri);


			int K = first.Count;

			if (K == 0)
			{
				Program.OutputText(" -- not compatible --", true);
				return;
			}

			// load components which stay unchanged
			List<ComponentGraphNode> nodesstatic = new List<ComponentGraphNode>();
			List<ComponentGraphNode> nodes2replace = new List<ComponentGraphNode>();
			List<CuboidMeshComponent> cmcsstatic = new List<CuboidMeshComponent>();
			foreach (SymmetryTriplet tri in first)
			{
				foreach (ComponentGraphNode node in tri.nodes)
				{
					if (!nodes2replace.Contains(node))
					{
						nodes2replace.Add(node);
					}
				}
			}
			foreach (ComponentGraphNode node in g.nodes)
			{
				if (!this.IsNondeContained(node, nodes2replace))
				{
					CuboidMeshComponent cmc = new CuboidMeshComponent(node, node.meshfile);
					node.tag = true;
					this.SetNodeReplace(node, node);
					cmcsstatic.Add(cmc);
				}
			}


			// now replace, find all permutations --> replace
			List<int[]> perms = new List<int[]>();
			int[] A = new int[K];
			for (int i = 0; i < K; ++i)
			{
				A[i] = first[i].compatibleTriplets.Count;
			}
			int[] temp = new int[K];

			Perms(temp, A, K, 0, perms);

			// each triplet has 3 nodes, they also form permutations
			List<int[]> perms_nodes = new List<int[]>();
			A = new int[K];
			temp = new int[K];
			for (int i = 0; i < K; ++i)
			{
				A[i] = 4;
			}

			Perms(temp, A, K, 0, perms_nodes);


			List<ComponentGraphNode[]> _tmpModelList = new List<ComponentGraphNode[]>();

			int N = perms.Count; // ignore the last one, which is identical with the original model
			for (int i = 0; i < N; ++i)
			{
				int[] perm = perms[i];

				// perm for nodes in triplets
				for (int j = 0; j < perms_nodes.Count; ++j)
				{
					int[] perm2 = perms_nodes[j];

					List<CuboidMeshComponent> cmcsdynamic = new List<CuboidMeshComponent>();
					cmcsdynamic.AddRange(cmcsstatic);

					// reset nodes properties
					foreach (ComponentGraph gg in targetgraphs)
					{
						foreach (ComponentGraphNode nd in gg.nodes)
						{
							this.SetNodeTag(nd,false);
							this.SetNodeReplace(nd,null);
							nd.cuboidMeshComponent = null;
							foreach (Slot s in nd.contactSlots)
								s.Reset();
						}
						if (gg.supernodes != null)
							foreach (ComponentGraphNode nd in gg.supernodes)
							{
								this.SetNodeTag(nd,false);
								this.SetNodeReplace(nd,null);
								nd.cuboidMeshComponent = null;
								foreach (Slot s in nd.contactSlots)
									s.Reset();
							}
					}
					foreach (ComponentGraphNode nd in nodes2replace)
					{
						this.SetNodeTag(nd, false);
						this.SetNodeReplace(nd, null);
						nd.cuboidMeshComponent = null;
					}
					foreach (ComponentGraphNode nd in g.nodes)
					{
						foreach (Slot s in nd.contactSlots)
							s.Reset();
					}

					// handle unreplaceable nodes -- those have multiple functionalities
					int index = 0;
					foreach (SymmetryTriplet tri in first)
					{
						SymmetryTriplet tri2 = tri.compatibleTriplets[perm[index]];
						for (int jj = 0; jj < 3; ++jj)
						{
							if (!this.IsNodesFunctionalCompatible(tri.nodes[jj], tri2.nodes[jj]) ||
								!this.IsNodesSlotsFunctionalCompatible(tri.nodes[jj], tri2.nodes[jj],tri,tri2)
								)
							{
								ComponentGraphNode node = tri.nodes[jj];
								List<CuboidMeshComponent> lscmc = this.CreateCuboidComponent(node);
								cmcsdynamic.AddRange(lscmc);
								this.SetNodeTag(node, true);
								this.SetNodeReplace(node,node);
							}
						}
						index++;
					}


					// do node replacement
					List<ComponentGraphNode> _nodes_2_replace = new List<ComponentGraphNode>();
					Dictionary<ComponentGraphNode, ComponentGraphNode> _nodes_map = new Dictionary<ComponentGraphNode, ComponentGraphNode>();
					index = 0;
					foreach (SymmetryTriplet tri in first)
					{
						SymmetryTriplet tri2 = tri.compatibleTriplets[perm[index]];
						
						List<ComponentGraphNode> _unodes = new List<ComponentGraphNode>();
						int perm_index = perm2[index];
						switch (perm_index)
						{
							case 1: // 2 symmetry nodes stay unchanged
								_unodes.Add(tri.nodes[0]);
								_unodes.Add(tri.nodes[1]);
								break;
							case 2: // 1 connecting nodes stay unchanged
								_unodes.Add(tri.nodes[2]);
								break;
							case 3: // all stay unchanged
								_unodes.AddRange(tri.nodes);
								break;
						}
						foreach (ComponentGraphNode node in _unodes)
						{
							if (!node.tag && !_nodes_map.ContainsKey(node))
							{
								cmcsdynamic.AddRange(this.CreateCuboidComponent(node));
								this.SetNodeTag(node, true);
								this.SetNodeReplace(node,node);
							}
						}
						
						// now replace the triplet
					//	this.ReplaceTriplet(tri, tri2, cmcsdynamic);
						int jj = 0;
						foreach (ComponentGraphNode nd in tri.nodes)
						{
							if (!nd.tag && !_nodes_map.ContainsKey(nd))
							{
								_nodes_2_replace.Add(nd);
								_nodes_map.Add(nd, tri2.nodes[jj]);
							}
							++jj;
						}

						index++;
					}

					
					// perform nodes replacements at once
					//----------------------------------------------------------------
					this.ReplaceNodes(_nodes_2_replace, _nodes_map, cmcsdynamic);
					//----------------------------------------------------------------


					if (!this.IsCompleteReplacing(g) && !this.IsSelfCloning(g))
					{
						// check existance
						int M = g.nodes.Count;
						ComponentGraphNode[] _repNodes = new ComponentGraphNode[M];
						for (int jj = 0; jj < M; ++jj)
						{
							_repNodes[jj] = g.nodes[jj].replace;
						}
						bool exist = false;
						foreach (ComponentGraphNode[] _existNodes in _tmpModelList)
						{
							int jj = 0;
							for (; jj < M && _repNodes[jj] == _existNodes[jj]; ++jj) ;
							if (jj == M)
							{
								exist = true;
								break;
							}
						}
						if (!exist)
						{
							_tmpModelList.Add(_repNodes);
							Program.OutputText(" ... ", true);
							
							this.synthesizedCuboidMeshComponents.Add(cmcsdynamic);
							this.synthesisTargetsFiles.Add(targetfile);

						}
					}
				}

			//	break;
			}
		}
		public void SortNodesWithReplaceOrder(ref List<ComponentGraphNode> nodes)
		{
			if (nodes.Count == 0) return;

			List<ComponentGraphNode> _nodes_has_connection = new List<ComponentGraphNode>();
			foreach (ComponentGraphNode node in nodes)
			{
				foreach (ComponentGraphEdge e in node.edges)
				{
					if (e.type == ComponentEdgeType.Symmetry) continue;
					ComponentGraphNode other = e.start == node ? e.end : e.start;
					if (other.replace != null) // already replaed
					{
						_nodes_has_connection.Add(node);
						break;
					}
				}
			}

			if (_nodes_has_connection.Count == 0) // all need replace
			{
				_nodes_has_connection.Add(nodes[0]);
			}

			foreach (ComponentGraphNode node in _nodes_has_connection)
			{
				nodes.Remove(node);
			}

			// handle rest
			Queue<ComponentGraphNode> Q = new Queue<ComponentGraphNode>();
			foreach (ComponentGraphNode node in nodes)
				Q.Enqueue(node);

			foreach (ComponentGraphNode node in _nodes_has_connection) node.tag = true;
			while (Q.Count() > 0)
			{
				ComponentGraphNode node = Q.Dequeue();
				bool replaceable = false;
				foreach (ComponentGraphEdge e in node.edges)
				{
					if (e.type == ComponentEdgeType.Symmetry) continue;
					ComponentGraphNode other = e.start == node ? e.end : e.start;
					if (other.tag) // already replaed
					{
						_nodes_has_connection.Add(node);
						node.tag = true;
						replaceable = true;
						break;
					}
				}
				if (!replaceable)
				{
					Q.Enqueue(node);
				}
			}
			foreach (ComponentGraphNode node in _nodes_has_connection) node.tag = false;
			nodes = _nodes_has_connection;
		}
		public void ReplaceNodes(List<ComponentGraphNode> nodes, Dictionary<ComponentGraphNode, ComponentGraphNode> map, List<CuboidMeshComponent> cmcs)
		{
			if (nodes.Count == 0)
			{
		//		Program.OutputText("no nodes to repalce!", true);
				return;
			}

			this.SortNodesWithReplaceOrder(ref nodes);

			foreach (ComponentGraphNode node in nodes)
			{
				if (!node.tag)
				{
					List<CuboidMeshComponent> cmc = null;
					List<CuboidMeshComponent> symmcmc = null;
					this.ReplaceNode(node, map[node], out cmc, out symmcmc);
					if (cmc != null) cmcs.AddRange(cmc);
					if (symmcmc != null) cmcs.AddRange(symmcmc);
				}
			}

		}
		public void FindSymmetryNodesCompatible(ComponentGraph g1, ComponentGraph g2)
		{
			// find in g2 compatible nodes for nodes in g1, we restrict one correspondence at most
			foreach (ComponentGraphNode node in g1.nodes) node.tag = false;
			for (int i = 0; i < g1.nodes.Count; ++i)
			{
				ComponentGraphNode node = g1.nodes[i];
				
				if (node.tag || node.symmetry == null) continue;
				
				node.compatibleNodes.Clear();
				node.symmetry.compatibleNodes.Clear();

				Vector3d nor1 = (node.pos - node.symmetry.pos).Normalize();

				double mindis = double.MaxValue; ComponentGraphNode tarnode = null;
				foreach (ComponentGraphNode node2 in g2.nodes)
				{
					if (node2.symmetry == null) continue;
					Vector3d nor2 = (node2.pos - node2.symmetry.pos).Normalize();

					double cos = nor1.Dot(nor2);
					if (cos < 0.5) continue; // not the same reglectional plane, or not the same side
					
					double dis = (node2.pos - node.pos).Length();
					if (dis < mindis)
					{
						mindis = dis;
						tarnode = node2;
					}
				}
				if (tarnode != null)
				{
					node.compatibleNodes.Add(tarnode);
					node.symmetry.compatibleNodes.Add(tarnode.symmetry);
					node.tag = node.symmetry.tag = true;
				}
			}
			foreach (ComponentGraphNode node in g1.nodes) node.tag = false;
		}
		public bool IsCompatible(SymmetryTriplet t1, SymmetryTriplet t2)
		{
			if (t1.functionality != t2.functionality || t1.topple != t2.topple) return false;
			return true;
		}
		public bool IsNodesFunctionalCompatible(ComponentGraphNode node1, ComponentGraphNode node2)
		{
		//	if (node1.functionalities.Count > node2.functionalities.Count) return false;
			var intersect = node1.functionalities.Intersect(node2.functionalities);
			int n = intersect.Count();
			if (n == node1.functionalities.Count) return true;
			
			// weak cases, here we weaken the contact functionality
			if (intersect.Contains(Functionality.supportive) && intersect.Contains(Functionality.anti_supportive)) return true;
			if (intersect.Contains(Functionality.supportive) && !node1.functionalities.Contains(Functionality.anti_supportive)) return true;
			if (intersect.Contains(Functionality.anti_supportive) && !node1.functionalities.Contains(Functionality.supportive)) return true;
			
			return false;
		}
		public bool IsNodesSlotsFunctionalCompatible(ComponentGraphNode node1, ComponentGraphNode node2, SymmetryTriplet tri1, SymmetryTriplet tri2)
		{
			// find if the nodes in triplets are slot-compatible, prevent from mapping undesirable results, e.g., a bed - a pilar
			Slot s1 = null;
			foreach (Slot s in node1.contactSlots)
			{
				ComponentGraphNode other = s.start == node1 ? s.end : s.start;
				if (tri1.nodes.Contains(other))
				{
					s1 = s;
					break;
				}
			}
			Slot s2 = null;
			foreach (Slot s in node2.contactSlots)
			{
				ComponentGraphNode other = s.start == node2 ? s.end : s.start;
				if (tri2.nodes.Contains(other))
				{
					s2 = s;
					break;
				}
			}
			return s1.contactPoints.Count == s2.contactPoints.Count;
		}
		public double ComputeCompatibility(SymmetryTriplet t1, SymmetryTriplet t2)
		{
			ComponentGraphNode node1 = t1.nodes[0], node2 = t1.nodes[2];
			ComponentGraphNode node3 = t2.nodes[0], node4 = t2.nodes[2];
			double sx = Math.Abs(node1.scale.x - node3.scale.x) / (node1.scale.x + node3.scale.x);
			double sy = Math.Abs(node1.scale.y - node3.scale.y) / (node1.scale.y + node3.scale.y);
			double sz = Math.Abs(node1.scale.z - node3.scale.z) / (node1.scale.z + node3.scale.z);
			double res = (1+sx)*(1+sy)*(1+sz);

			sx = Math.Abs(node2.scale.x - node4.scale.x) / (node2.scale.x + node4.scale.x);
			sy = Math.Abs(node2.scale.y - node4.scale.y) / (node2.scale.y + node4.scale.y);
			sz = Math.Abs(node2.scale.z - node4.scale.z) / (node2.scale.z + node4.scale.z);
			res *= (1+sx)*(1+sy)*(1+sz);

			Vector3d u1 = node1.pos, u2 = node2.pos, u3 = t1.nodes[1].pos;
			Vector3d v1 = node3.pos, v2 = node4.pos, v3 = t2.nodes[1].pos;
			Vector3d u12 = u2 - u1, u13 = u3 - u1, u23 = u3 - u2;
			Vector3d v12 = v2 - v1, v13 = v3 - v1, v23 = v3 - v2;
			double c1 = u12.Length() * 2 + u12.Length();
			double c2 = v12.Length() * 2 + v12.Length();
			double cos1 = u12.Dot(u23) / (u12.Length() * u23.Length());
			double cos2 = v12.Dot(v23) / (v12.Length() * v23.Length());
			double theta1 = Math.Acos(cos1);
			double theta2 = Math.Acos(cos2);

			double cos = Math.Abs(theta1 - theta2) / (theta1 + theta2);
			double lengthRatio = Math.Acos(c1 - c2) / (c1 + c2);

			return res *(1 + cos) * (1 + lengthRatio);
		}


		// static members
        static void Swap(List<string> lsArray, int x, int y)
        {
            string t = lsArray[x];
            lsArray[x] = lsArray[y];
            lsArray[y] = t;
        }
		public static void Perms(int[] a, int[] bArray, int n, int temp, List<int[]> perms)
		{
			if (temp < n)
			{
				for (int j = 0; j < bArray[temp];)
				{
					Perms(a, bArray, n, temp+1, perms);

					a[temp] = ++j;
					if (a[temp] == bArray[temp])
						a[temp] = 0;

				}
			}
			else
			{
				perms.Add(a.Clone() as int[]);
			}
		}
		public static Dictionary<string, int> objfile2vocindexlookupmap = null;
		public static Dictionary<int, string> index2imgfilelookupmap = null;
		public static Dictionary<string, int> graphlookupmap = null;
		public static Dictionary<string, ComponentGraphNode> graphnodelookupmap = null;
		public static int[] topmatches = null;
		public static List<ComponentGraph> spatialgraphs = null;
		public static void LoadSpatialGraphs()
		{
			string file = @".\data\spatialGraphs.txt";

			if (!File.Exists(file)) return;

			char[] delimiters = { ' ', '\t' };
			string s = "";


			spatialgraphs = new List<ComponentGraph>();
			graphlookupmap = new Dictionary<string, int>();
			graphnodelookupmap = new Dictionary<string, ComponentGraphNode>();

			StreamReader sr = new StreamReader(file);
			ComponentGraph g = null; int graphindex = 0;
			while (sr.Peek() > -1)
			{
				s = sr.ReadLine();

				int index = 0;

				string[] tokens = s.Split(delimiters);
				switch (tokens[0].ToLower())
				{
					case "#":
						graphlookupmap.Add(tokens[1], graphindex++);
						if (g != null)
						{
							foreach (ComponentGraphEdge e in g.edges)
							{
								if (e.type != ComponentEdgeType.Underneath && e.type != ComponentEdgeType.Upon)
								{
									e.start.edges.Add(e);
									e.end.edges.Add(e);
								}
								else
									e.start.edges.Add(e);
							}
							g.classname = g.nodes[0].name;
							spatialgraphs.Add(g);
						}

						g = new ComponentGraph();
						g.meshfile = tokens[1];
						index = 0;

						break;
					case "n":
						{
							ComponentGraphNode node = new ComponentGraphNode(index++);
							node.meshfile = tokens[1];
							graphnodelookupmap.Add(node.meshfile, node);
							node.name = tokens[2];
							node.componentIndex = int.Parse(tokens[3]);
							node.pos = new Vector3d(double.Parse(tokens[4]),
									double.Parse(tokens[5]),
									double.Parse(tokens[6])
								);
							node.scale = new Vector3d(double.Parse(tokens[7]),
									double.Parse(tokens[8]),
									double.Parse(tokens[9])
								);
							g.nodes.Add(node);
						}
						break;
					case "e":
						{
							int r = int.Parse(tokens[1]);
							int t = int.Parse(tokens[2]);
							bool dirt = bool.Parse(tokens[4]);
							ComponentGraphEdge edge = new ComponentGraphEdge(g.nodes[r],g.nodes[t],dirt);
							switch (tokens[3])
							{
								case "Contact":
									edge.type = ComponentEdgeType.Contact;
									break;
								case "Symmetry":
									edge.type = ComponentEdgeType.Symmetry;
									break;
								case "Underneath":
									edge.type = ComponentEdgeType.Underneath;
									break;
								case "Upon":
									edge.type = ComponentEdgeType.Upon;
									break;
							}
							g.edges.Add(edge);
						}
						break;
				}
			}

			// last graph
			if (g != null)
			{
				foreach (ComponentGraphEdge e in g.edges)
				{
					if (e.type != ComponentEdgeType.Underneath && e.type != ComponentEdgeType.Upon)
					{
						e.start.edges.Add(e);
						e.end.edges.Add(e);
					}
					else
						e.start.edges.Add(e);
				}
				g.classname = g.nodes[0].name;
				spatialgraphs.Add(g);
			}

			foreach (ComponentGraph gg in spatialgraphs)
			{
				foreach (ComponentGraphNode node in gg.nodes)
				{
					node.getSymmetry();
				}
			}

			sr.Close();
		}
		public static void LocateComponentSlotPoints()
		{
			if (spatialgraphs == null) return;

			Program.OutputText("compute contacts points ...",true);

			foreach (ComponentGraph g in spatialgraphs)
			{
				// load the meshes
				List<CuboidMeshComponent> meshes = new List<CuboidMeshComponent>();
				foreach (ComponentGraphNode node in g.nodes)
				{
					StreamReader sr = new StreamReader(node.meshfile);
					Mesh m = new Mesh(sr, false);
					CuboidMeshComponent cmc = new CuboidMeshComponent(node, m);
					meshes.Add(cmc);
					sr.Close();
				}

				// ground contact slots
				double miny = double.MaxValue; Polygon plg = null;
				foreach (CuboidMeshComponent cmc in meshes)
				{
					foreach (Polygon p in cmc.boundingBox.polygons)
					{
						double yy = p.center.pos.y;
						if (yy < miny)
						{
							miny = yy;
							plg = p;
						}
					}
				}

				double thres2 = 1.0 / 64;
				foreach (CuboidMeshComponent cmc in meshes)
				{
					foreach (Polygon p in cmc.boundingBox.polygons)
					{
						double yy = p.center.pos.y;
						if (Math.Abs(yy - miny) < thres2)
						{
							cmc.node.isgroundcontact = true;
						}
					}
				}


				int slotindex = 0;
				foreach (ComponentGraphEdge e in g.edges)
				{
					if (e.type == ComponentEdgeType.Underneath || e.type == ComponentEdgeType.Contact)
					{
						CuboidMeshComponent m1 = meshes[e.start.componentIndex], m2 = meshes[e.end.componentIndex];
						
						// find contact faces
						Polygon face = null; Mesh m = null;
						Polygon face1 = null, face2 = null; double mindis = double.MaxValue;
						foreach (Polygon plg1 in m1.boundingBox.polygons)
						{
							Vector3d c1 = plg1.center.pos; Vector3d n1 = plg1.normal;
							foreach (Polygon plg2 in m2.boundingBox.polygons)
							{
								Vector3d c2 = plg2.center.pos; Vector3d n2 = plg2.normal;
								if (Math.Abs(n1.Dot(n2)) < 0.3) // not parallel
									continue;
								if (PointProjectionInPolygon(c1, plg2) || PointProjectionInPolygon(c2, plg1))
								{
									double dis = Math.Abs((c2 - c1).Dot(n1));
									if (dis < mindis)
									{
										mindis = dis;
										face1 = plg1;
										face2 = plg2;
									}
								}
							}
						}

						if (face1 == null && face2 == null)
						{
							// not found!!!
							double min = double.MaxValue;
							foreach (Polygon plg1 in m1.boundingBox.polygons)
							{
								Vector3d c1 = plg1.center.pos;
								foreach (Point3 pt in m2.boundingBox.points3)
								{
									double dis = (pt.pos - c1).Length();
									if (dis < min)
									{
										m = m1.mesh;
										face = plg1;
										min = dis;
									}
								}
							}
							foreach (Polygon plg1 in m2.boundingBox.polygons)
							{
								Vector3d c1 = plg1.center.pos;
								foreach (Point3 pt in m1.boundingBox.points3)
								{
									double dis = (pt.pos - c1).Length();
									if (dis < min)
									{
										m = m2.mesh;
										face = plg1;
										min = dis;
									}
								}
							}
						}
						else
						{
							double s1 = (face1.points3[1].pos - face1.points3[0].pos).Length() *
								(face1.points3[2].pos - face1.points3[1].pos).Length();
							double s2 = (face2.points3[1].pos - face2.points3[0].pos).Length() *
								(face2.points3[2].pos - face2.points3[1].pos).Length();
							if (s1 < s2)
							{
								face = face1;
								m = m1.mesh;
							}
							else
							{
								face = face2;
								m = m2.mesh;
							}
						}


						double thres = 1.0 / 64;
						List<Vector3d> slotpoints = new List<Vector3d>();
						while (slotpoints.Count == 0)
						{
							for (int i = 0, j = 0; i < m.VertexCount; ++i, j += 3)
							{
								Vector3d pos = new Vector3d(m.VertexPos, j);
								double dis = Math.Abs((pos - face.center.pos).Dot(face.normal));
								if (dis < thres)
								{
									slotpoints.Add(pos);
								}
							}
							thres *= 2;
						}

						Slot slot = new Slot(e.start, e.end, slotpoints, slotindex++);
	
						e.start.contactSlots.Add(slot);
						e.end.contactSlots.Add(slot);
					}
				}

				GC.Collect();
			}

			Program.OutputText("finished ...",true);

		}
		public static void LoadComponentSlots()
		{
			if (spatialgraphs == null) return;

			string file = @".\data\slotcontacts.txt";

			if (!File.Exists(file))
				return;

			char[] delimiters = { ' ', '\t' };
			string s = "";

			StreamReader sr = new StreamReader(file);
			while (sr.Peek() > -1)
			{
				s = sr.ReadLine();
				string[] tokens = s.Split(delimiters);
				if (tokens[0] == "#")
				{
					int N = int.Parse(tokens[3]);
					N = graphlookupmap[tokens[1]];
					ComponentGraph g = spatialgraphs[N];
					int slotCount = int.Parse(tokens[2]);

					List<Slot> slots = new List<Slot>();
					for (int it = 0; it < slotCount; ++it)
					{
						s = sr.ReadLine();
						tokens = s.Split(delimiters);

						int slotindex = int.Parse(tokens[0]);
						int nodeIndex1 = int.Parse(tokens[1]);
						int nodeIndex2 = int.Parse(tokens[2]);

						ComponentGraphNode node1 = g.nodes[nodeIndex1];
						ComponentGraphNode node2 = g.nodes[nodeIndex2];

						List<double> pos = new List<double>();
						for (int j = 3; j < tokens.Length; ++j)
						{
							if (tokens[j] != "")
							{
								pos.Add(double.Parse(tokens[j]));
							}
						}
						List<Vector3d> points = new List<Vector3d>();
						double[] vpos = pos.ToArray();
						for (int j = 0; j < vpos.Length; j += 3)
						{
							Vector3d v = new Vector3d(vpos, j);
							points.Add(v);
						}

						List<Vector3d> pts = new List<Vector3d>();
						pts.AddRange(points);
						Slot slot1 = new Slot(node1, node2, points, slotindex);
						Slot slot2 = new Slot(node1, node2, pts, slotindex);
						slot1.clonned = slot2;
						slot2.clonned = slot1;
						node1.contactSlots.Add(slot1);
						node2.contactSlots.Add(slot2);
					}
				}
			}

			// handle symmetry
			foreach (ComponentGraph g in spatialgraphs)
			{
				foreach (ComponentGraphNode node in g.nodes)
				{
					if (node.symmetry != null)
					{
						foreach (Slot ss in node.contactSlots)
						{
							if (ss.symmetry != null) continue;
							foreach (Slot tt in node.symmetry.contactSlots)
							{
								if (tt.symmetry != null) continue;
								if ((ss.start.symmetry == tt.start && ss.end.symmetry == tt.end) ||
									(ss.start.symmetry == tt.start && ss.end == tt.end) ||
									(ss.start.symmetry == tt.end && ss.end.symmetry == tt.start) ||
									(ss.start.symmetry == tt.end && ss.end == tt.start) ||
									(ss.start == tt.start && ss.end.symmetry == tt.end) ||
									(ss.start == tt.end && ss.end.symmetry == tt.start)
								)
								{
									ss.symmetry = tt;
									tt.symmetry = ss;
									ss.clonned.symmetry = tt.clonned;
									tt.clonned.symmetry = ss.clonned;
								}
							}
						}
					}
				}
			}

			sr.Close();
		}
		public static void LoadSymmetryTriplets()
		{
			if (spatialgraphs == null) return;

			string file = @".\data\symmetryTriplets.txt";

			if (!File.Exists(file))
				return;

			char[] delimiters = { ' ', '\t' };
			string s = "";

			StreamReader sr = new StreamReader(file);
			while (sr.Peek() > -1)
			{
				s = sr.ReadLine();
				string[] tokens = s.Split(delimiters);
				if (tokens[0] == "#")
				{
					int N = int.Parse(tokens[2]);
					N = graphlookupmap[tokens[1]];
					ComponentGraph g = spatialgraphs[N];
					g.symmtryTriplets = new List<SymmetryTriplet>();
					int n = int.Parse(tokens[3]);
					for (int i = 0; i < n; ++i)
					{
						s = sr.ReadLine();
						tokens = s.Split(delimiters);

						int index1 = int.Parse(tokens[0]);
						int index2 = int.Parse(tokens[1]);
						int index3 = int.Parse(tokens[2]);

						ComponentGraphNode node1 = g.nodes[index1];
						ComponentGraphNode node2 = g.nodes[index2];
						ComponentGraphNode node3 = g.nodes[index3];

						Functionality type = Functionality.none;
						switch (tokens[3].ToLower())
						{
							case "contacting":
								type = Functionality.contacting;
								break;
							case "supportive":
								type = Functionality.supportive;
								break;
							case "anti_supportive":
								type = Functionality.anti_supportive;
								break;
						}

						SymmetryTriplet tri = new SymmetryTriplet(node1, node2, node3, type);
						switch (tokens[4].ToLower())
						{
							case "true":
								tri.isgroundtouching = true;
								break;
							case "false":
								tri.isgroundtouching = false;
								break;
						}
						switch (tokens[5].ToLower())
						{
							case "true":
								tri.topple = true;
								break;
							case "false":
								tri.topple = false;
								break;
						}

						g.symmtryTriplets.Add(tri);
					}
				}
			}

			sr.Close();

		}
		public static void ComputeGraphNodeFunctionalities()
		{
			if (spatialgraphs == null) return;
			foreach (ComponentGraph g in spatialgraphs)
			{
				if (g.symmtryTriplets == null)
					g.GetSymmetryTriplets();
				g.ComputeNodeFunctionalities();
			}
		}
		public static void LoadGroundContactSlots()
		{
			if (spatialgraphs == null) return;

			string file = @".\data\slotgroundcontacts.txt";

			if (!File.Exists(file))
				return;

			char[] delimiters = { ' ', '\t' };
			string s = "";

			StreamReader sr = new StreamReader(file);
			while (sr.Peek() > -1)
			{
				s = sr.ReadLine();
				string[] tokens = s.Split(delimiters);
				if (tokens[0] == "#")
				{
					int N = int.Parse(tokens[2]);
					N = graphlookupmap[tokens[1]];
					ComponentGraph g = spatialgraphs[N];
					int n = int.Parse(tokens[3]);
					for (int i = 0; i < n; ++i)
					{
						ComponentGraphNode node = g.nodes[i];
						s = sr.ReadLine();
						tokens = s.Split(delimiters);
						switch (tokens[0].ToLower())
						{
							case "true":
								node.isgroundcontact = true;
								break;
							case "false":
								node.isgroundcontact = false;
								break;
						}
					}
				}
			}

			sr.Close();
		}
		public static bool PointProjectionInPolygon(Vector3d p, Polygon plg)
		{
			Vector3d c = plg.center.pos; Vector3d n1 = plg.normal;
			Vector3d e1 = plg.points3[1].pos - plg.points3[0].pos;
			Vector3d e2 = plg.points3[2].pos - plg.points3[1].pos;
			double l1 = e1.Length() / 2;
			double l2 = e2.Length() / 2;
			Vector3d pc = p - c;
			return Math.Abs(pc.Dot(e1.Normalize())) < l1 &&
				Math.Abs(pc.Dot(e2.Normalize())) < l2;
		}
		public static bool PointInTriangle(Vector3d p, Vector3d a, Vector3d b, Vector3d c)
		{
			Vector3d pa = p - a;
			Vector3d pb = p - b;
			Vector3d pc = p - c;
			Vector3d v1 = pa.Cross(pb);
			Vector3d v2 = pb.Cross(pc);
			Vector3d v3 = pc.Cross(pa);
			double sign1 = v1.Dot(v2);
			double sign2 = v2.Dot(v3);
			double sign3 = v3.Dot(v1);
			return (sign1 <= 0 && sign2 <= 0 && sign3 <= 0) || 
				(sign1 >= 0 && sign2 >= 0 && sign3 >= 0);
		}

	}
}
