Scalar Data Visualization: Part 3

//-----------------------------------------------------------------------------
// File: ISOsurface.cs
//
// DemolitiOn is My MissiOn
//
// Copyright (c) Ibrahim Hamed. All rights reserved.
//-----------------------------------------------------------------------------

using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.DirectX;
using Microsoft.DirectX.Direct3D;

namespace ScoOteRLoAdeR
{
	class ISOsurface
	{
		protected ContourManager Parent;
		protected VertexBuffer VB;
		protected CustomVertex.PositionNormalColored[] Verts;
		protected List<Point> points;
		protected int TriangleCount;

		public int Color;
		public double Value;

		public ISOsurface(ContourManager parent, double value, int color)
		{
			Parent = parent;
			Value = value;
			Color = color;
			points = new List<Point>();
		}

		public virtual bool CreateSurface()
		{
			return false;
		}

		public bool SetVertexBuffer()
		{
			if (points.Count == 0)
				return false;

			VB = new VertexBuffer(typeof(CustomVertex.PositionNormalColored), points.Count /* 2*/, Parent.Engine.CurrentDevice,
				Usage.WriteOnly, CustomVertex.PositionNormalColored.Format, Pool.Managed);
			Verts = new CustomVertex.PositionNormalColored[points.Count];// * 2];
			TriangleCount = points.Count / 3;
			for (int i = 0; i < points.Count; i++)
			{
				Verts[i] = new CustomVertex.PositionNormalColored();
				Verts[i].X = (float)points[i].Data[0];
				Verts[i].Y = (float)points[i].Data[1];
				Verts[i].Z = (float)points[i].Data[2];// +0.05f;
				Verts[i].Nx = points[i].NormalX;
				Verts[i].Ny = points[i].NormalY;
				Verts[i].Nz = points[i].NormalZ;
				Verts[i].Color = Color;
			}
			//for (int i = 0; i < points.Count; i++)
			//{
			//    Verts[i + points.Count] = new CustomVertex.PositionNormalColored();
			//    Verts[i + points.Count].X = (float)points[i].Data[0];
			//    Verts[i + points.Count].Y = (float)points[i].Data[1];
			//    Verts[i + points.Count].Z = (float)points[i].Data[2];// +0.05f;
			//    Verts[i + points.Count].Nx = -points[i].NormalX;
			//    Verts[i + points.Count].Ny = -points[i].NormalY;
			//    Verts[i + points.Count].Nz = -points[i].NormalZ;
			//    Verts[i + points.Count].Color = Color;
			//}
			VB.SetData(Verts, 0, 0);
			points.Clear();
			Verts = null;
			return true;
		}

		public void Render()
		{
			Parent.Engine.CurrentDevice.SetStreamSource(0, VB, 0);
			Parent.Engine.CurrentDevice.VertexFormat = CustomVertex.PositionNormalColored.Format;
			Parent.Engine.CurrentDevice.DrawPrimitives(PrimitiveType.TriangleList, 0, TriangleCount);
		}
	}

	class MarchingCubes : ISOsurface
	{
		public MarchingCubes(ContourManager parent, double value, int color)
			: base(parent, value, color)
		{ }

		public override bool CreateSurface()
		{
			foreach (MarchedCube c in Parent.Cubes)
			{
				uint TableIndex = 0;
				if (Parent.Points[c.Vertex1].Data[Parent.VariableIndex] < Value)
					TableIndex |= 1;
				if (Parent.Points[c.Vertex2].Data[Parent.VariableIndex] < Value)
					TableIndex |= 2;
				if (Parent.Points[c.Vertex3].Data[Parent.VariableIndex] < Value)
					TableIndex |= 4;
				if (Parent.Points[c.Vertex4].Data[Parent.VariableIndex] < Value)
					TableIndex |= 8;
				if (Parent.Points[c.Vertex5].Data[Parent.VariableIndex] < Value)
					TableIndex |= 16;
				if (Parent.Points[c.Vertex6].Data[Parent.VariableIndex] < Value)
					TableIndex |= 32;
				if (Parent.Points[c.Vertex7].Data[Parent.VariableIndex] < Value)
					TableIndex |= 64;
				if (Parent.Points[c.Vertex8].Data[Parent.VariableIndex] < Value)
					TableIndex |= 128;
				if (Parent.EdgeTable[TableIndex] != 0)
				{
					for (int i = 0; Parent.TriTable[TableIndex,i] != -1; i += 3)
					{
						int p1 = CreatePoint(Parent.TriTable[TableIndex,i], Value, c);
						int p2 = CreatePoint(Parent.TriTable[TableIndex,i + 1], Value, c);
						int p3 = CreatePoint(Parent.TriTable[TableIndex,i + 2], Value, c);
						CreateFaceNormals(p1, p2, p3);
					}
				}
			}
			return SetVertexBuffer();
		}

		private int CreatePoint(int EdgeNum, double Value, MarchedCube c)
		{
			int [] pt = GetVerts(EdgeNum,c);

			double diff = (Parent.Points[pt[1]].Data[Parent.VariableIndex] - Parent.Points[pt[0]].Data[Parent.VariableIndex]);
			double t = 0.0;
			if (diff != 0)
				t = (Value - Parent.Points[pt[0]].Data[Parent.VariableIndex]) / diff;

			Point p = new Point(13);
			double DeltaX = (Parent.Points[pt[1]].Data[0] - Parent.Points[pt[0]].Data[0]);
			double DeltaY = (Parent.Points[pt[1]].Data[1] - Parent.Points[pt[0]].Data[1]);
			double DeltaZ = (Parent.Points[pt[1]].Data[2] - Parent.Points[pt[0]].Data[2]);
			p.Data.Add(Parent.Points[pt[0]].Data[0] + t * DeltaX);
			p.Data.Add(Parent.Points[pt[0]].Data[1] + t * DeltaY);
			p.Data.Add(Parent.Points[pt[0]].Data[2] + t * DeltaZ);
			points.Add(p);

			return points.Count - 1;
		}

		private int[] GetVerts(int EdgeNum, MarchedCube c)
		{
			int[] pts = new int[2];
			double[] sorter = new double[2];
			switch (EdgeNum)
			{
				case 0:
					pts[0] = c.Vertex1;
					pts[1] = c.Vertex2;
					break;
				case 1:
					pts[0] = c.Vertex2;
					pts[1] = c.Vertex3;
					break;
				case 2:
					pts[0] = c.Vertex3;
					pts[1] = c.Vertex4;
					break;
				case 3:
					pts[0] = c.Vertex4;
					pts[1] = c.Vertex1;
					break;
				case 4:
					pts[0] = c.Vertex5;
					pts[1] = c.Vertex6;
					break;
				case 5:
					pts[0] = c.Vertex6;
					pts[1] = c.Vertex7;
					break;
				case 6:
					pts[0] = c.Vertex7;
					pts[1] = c.Vertex8;
					break;
				case 7:
					pts[0] = c.Vertex8;
					pts[1] = c.Vertex5;
					break;
				case 8:
					pts[0] = c.Vertex1;
					pts[1] = c.Vertex5;
					break;
				case 9:
					pts[0] = c.Vertex6;
					pts[1] = c.Vertex2;
					break;
				case 10:
					pts[0] = c.Vertex7;
					pts[1] = c.Vertex3;
					break;
				case 11:
					pts[0] = c.Vertex8;
					pts[1] = c.Vertex4;
					break;
			}
			sorter[0] = Parent.Points[pts[0]].Data[Parent.VariableIndex];
			sorter[1] = Parent.Points[pts[1]].Data[Parent.VariableIndex];
			Array.Sort(sorter, pts);
			return pts;
		}

		private void CreateFaceNormals(int P1, int P2, int P3)
		{
			Point p1 = points[P1];
			Point p2 = points[P2];
			Point p3 = points[P3];

			Vector3 v1 = new Vector3((float)(p2.Data[0] - p1.Data[0]),
									(float)(p2.Data[1] - p1.Data[1]),
									(float)(p2.Data[2] - p1.Data[2]));

			Vector3 v2 = new Vector3((float)(p3.Data[0] - p1.Data[0]),
									(float)(p3.Data[1] - p1.Data[1]),
									(float)(p3.Data[2] - p1.Data[2]));

			v1 = Vector3.Cross(v2, v1);
			v1.Normalize();
			p1.NormalX = p2.NormalX = p3.NormalX = v1.X;
			p1.NormalY = p2.NormalY = p3.NormalY = v1.Y;
			p1.NormalZ = p2.NormalZ = p3.NormalZ = v1.Z;

			points[P1] = p1;
			points[P2] = p2;
			points[P3] = p3;
		}
	}
}