/Main/src/DynamicDataDisplay/Charts/Isolines/IsolineBuilder.cs

using System;

using System.Linq;

using System.Diagnostics;

using System.Runtime.Serialization;

using System.Windows;

using System.Windows.Media;

using Microsoft.Research.DynamicDataDisplay.Common.Auxiliary;

using Microsoft.Research.DynamicDataDisplay.DataSources;

using System.Collections.Generic;



namespace Microsoft.Research.DynamicDataDisplay.Charts.Isolines

{

	/// <summary>

	/// Generates geometric object for isolines of the input 2d scalar field.

	/// </summary>

	public sealed class IsolineBuilder

	{

		/// <summary>

		/// The density of isolines means the number of levels to draw.

		/// </summary>

		private int density = 12;



		private bool[,] processed;



		/// <summary>Number to be treated as missing value. NaN if no missing value is specified</summary>

		private double missingValue = Double.NaN;



		static IsolineBuilder()

		{

			SetCellDictionaries();

		}



		/// <summary>

		/// Initializes a new instance of the <see cref="IsolineBuilder"/> class.

		/// </summary>

		public IsolineBuilder() { }



		/// <summary>

		/// Initializes a new instance of the <see cref="IsolineBuilder"/> class for specified 2d scalar data source.

		/// </summary>

		/// <param name="dataSource">The data source with 2d scalar data.</param>

		public IsolineBuilder(IDataSource2D<double> dataSource)

		{

			DataSource = dataSource;

		}



		public double MissingValue

		{

			get

			{

				return missingValue;

			}

			set

			{

				missingValue = value;

			}

		}



		#region Private methods



		private static Dictionary<int, Dictionary<int, Edge>> dictChooser = new Dictionary<int, Dictionary<int, Edge>>();

		private static void SetCellDictionaries()

		{

			var bottomDict = new Dictionary<int, Edge>();

			bottomDict.Add((int)CellBitmask.RightBottom, Edge.Right);

			bottomDict.Add(Edge.Left,

				CellBitmask.LeftTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom | CellBitmask.RightTop,

				CellBitmask.LeftTop | CellBitmask.RightBottom | CellBitmask.RightTop,

				CellBitmask.LeftBottom);

			bottomDict.Add(Edge.Right,

				CellBitmask.RightTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom | CellBitmask.LeftTop,

				CellBitmask.LeftBottom | CellBitmask.LeftTop | CellBitmask.RightTop);

			bottomDict.Add(Edge.Top,

				CellBitmask.RightBottom | CellBitmask.RightTop,

				CellBitmask.LeftBottom | CellBitmask.LeftTop);



			var leftDict = new Dictionary<int, Edge>();

			leftDict.Add(Edge.Top,

				CellBitmask.LeftTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom | CellBitmask.RightTop);

			leftDict.Add(Edge.Right,

				CellBitmask.LeftTop | CellBitmask.RightTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom);

			leftDict.Add(Edge.Bottom,

				CellBitmask.RightBottom | CellBitmask.RightTop | CellBitmask.LeftTop,

				CellBitmask.LeftBottom);



			var topDict = new Dictionary<int, Edge>();

			topDict.Add(Edge.Right,

				CellBitmask.RightTop,

				CellBitmask.LeftTop | CellBitmask.LeftBottom | CellBitmask.RightBottom);

			topDict.Add(Edge.Right,

				CellBitmask.RightBottom,

				CellBitmask.LeftTop | CellBitmask.LeftBottom | CellBitmask.RightTop);

			topDict.Add(Edge.Left,

				CellBitmask.RightBottom | CellBitmask.RightTop | CellBitmask.LeftTop,

				CellBitmask.LeftBottom,

				CellBitmask.LeftTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom | CellBitmask.RightTop);

			topDict.Add(Edge.Bottom,

				CellBitmask.RightBottom | CellBitmask.RightTop,

				CellBitmask.LeftTop | CellBitmask.LeftBottom);



			var rightDict = new Dictionary<int, Edge>();

			rightDict.Add(Edge.Top,

				CellBitmask.RightTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom | CellBitmask.LeftTop);

			rightDict.Add(Edge.Left,

				CellBitmask.LeftTop | CellBitmask.RightTop,

				CellBitmask.LeftBottom | CellBitmask.RightBottom);

			rightDict.Add(Edge.Bottom,

				CellBitmask.RightBottom,

				CellBitmask.LeftTop | CellBitmask.LeftBottom | CellBitmask.RightTop);



			dictChooser.Add((int)Edge.Left, leftDict);

			dictChooser.Add((int)Edge.Right, rightDict);

			dictChooser.Add((int)Edge.Bottom, bottomDict);

			dictChooser.Add((int)Edge.Top, topDict);

		}



		private Edge GetOutEdge(Edge inEdge, ValuesInCell cv, IrregularCell rect, double value)

		{

			// value smaller than all values in corners or 

			// value greater than all values in corners

			if (!cv.ValueBelongTo(value))

			{

				throw new IsolineGenerationException(Strings.Exceptions.IsolinesValueIsOutOfCell);

			}



			CellBitmask cellVal = cv.GetCellValue(value);

			var dict = dictChooser[(int)inEdge];

			if (dict.ContainsKey((int)cellVal))

			{

				Edge result = dict[(int)cellVal];

				switch (result)

				{

					case Edge.Left:

						if (cv.LeftTop.IsNaN() || cv.LeftBottom.IsNaN())

							result = Edge.None;

						break;

					case Edge.Right:

						if (cv.RightTop.IsNaN() || cv.RightBottom.IsNaN())

							result = Edge.None;

						break;

					case Edge.Top:

						if (cv.RightTop.IsNaN() || cv.LeftTop.IsNaN())

							result = Edge.None;

						break;

					case Edge.Bottom:

						if (cv.LeftBottom.IsNaN() || cv.RightBottom.IsNaN())

							result = Edge.None;

						break;

				}

				return result;

			}

			else if (cellVal.IsDiagonal())

			{

				return GetOutForOpposite(inEdge, cellVal, value, cv, rect);

			}



			const double near_zero = 0.0001;

			const double near_one = 1 - near_zero;



			double lt = cv.LeftTop;

			double rt = cv.RightTop;

			double rb = cv.RightBottom;

			double lb = cv.LeftBottom;



			switch (inEdge)

			{

				case Edge.Left:

					if (value == lt)

						value = near_one * lt + near_zero * lb;

					else if (value == lb)

						value = near_one * lb + near_zero * lt;

					else

						return Edge.None;

					// Now this is possible because of missing value

					//throw new IsolineGenerationException(Strings.Exceptions.IsolinesUnsupportedCase);

					break;

				case Edge.Top:

					if (value == rt)

						value = near_one * rt + near_zero * lt;

					else if (value == lt)

						value = near_one * lt + near_zero * rt;

					else

						return Edge.None;

					// Now this is possibe because of missing value

					//throw new IsolineGenerationException(Strings.Exceptions.IsolinesUnsupportedCase);

					break;

				case Edge.Right:

					if (value == rb)

						value = near_one * rb + near_zero * rt;

					else if (value == rt)

						value = near_one * rt + near_zero * rb;

					else

						return Edge.None;

					// Now this is possibe because of missing value

					//throw new IsolineGenerationException(Strings.Exceptions.IsolinesUnsupportedCase);

					break;

				case Edge.Bottom:

					if (value == rb)

						value = near_one * rb + near_zero * lb;

					else if (value == lb)

						value = near_one * lb + near_zero * rb;

					else

						return Edge.None;

					// Now this is possibe because of missing value

					//throw new IsolineGenerationException(Strings.Exceptions.IsolinesUnsupportedCase);

					break;

			}



			// Recursion?

			//return GetOutEdge(inEdge, cv, rect, value);



			return Edge.None;

		}



		private Edge GetOutForOpposite(Edge inEdge, CellBitmask cellVal, double value, ValuesInCell cellValues, IrregularCell rect)

		{

			Edge outEdge;



			SubCell subCell = GetSubCell(inEdge, value, cellValues);



			int iters = 1000; // max number of iterations

			do

			{

				ValuesInCell subValues = cellValues.GetSubCell(subCell);

				IrregularCell subRect = rect.GetSubRect(subCell);

				outEdge = GetOutEdge(inEdge, subValues, subRect, value);

				if (outEdge == Edge.None)

					return Edge.None;

				bool isAppropriate = subCell.IsAppropriate(outEdge);

				if (isAppropriate)

				{

					ValuesInCell sValues = subValues.GetSubCell(subCell);



					Point point = GetPointXY(outEdge, value, subValues, subRect);

					segments.AddPoint(point);

					return outEdge;

				}

				else

				{

					subCell = GetAdjacentEdge(subCell, outEdge);

				}



				byte e = (byte)outEdge;

				inEdge = (Edge)((e > 2) ? (e >> 2) : (e << 2));

				iters--;

			} while (iters >= 0);



			throw new IsolineGenerationException(Strings.Exceptions.IsolinesDataIsUndetailized);

		}



		private static SubCell GetAdjacentEdge(SubCell sub, Edge edge)

		{

			SubCell res = SubCell.LeftBottom;



			switch (sub)

			{

				case SubCell.LeftBottom:

					res = edge == Edge.Top ? SubCell.LeftTop : SubCell.RightBottom;

					break;

				case SubCell.LeftTop:

					res = edge == Edge.Bottom ? SubCell.LeftBottom : SubCell.RightTop;

					break;

				case SubCell.RightBottom:

					res = edge == Edge.Top ? SubCell.RightTop : SubCell.LeftBottom;

					break;

				case SubCell.RightTop:

				default:

					res = edge == Edge.Bottom ? SubCell.RightBottom : SubCell.LeftTop;

					break;

			}



			return res;

		}



		private static SubCell GetSubCell(Edge inEdge, double value, ValuesInCell vc)

		{

			double lb = vc.LeftBottom;

			double rb = vc.RightBottom;

			double rt = vc.RightTop;

			double lt = vc.LeftTop;



			SubCell res = SubCell.LeftBottom;

			switch (inEdge)

			{

				case Edge.Left:

					res = (Math.Abs(value - lb) < Math.Abs(value - lt)) ? SubCell.LeftBottom : SubCell.LeftTop;

					break;

				case Edge.Top:

					res = (Math.Abs(value - lt) < Math.Abs(value - rt)) ? SubCell.LeftTop : SubCell.RightTop;

					break;

				case Edge.Right:

					res = (Math.Abs(value - rb) < Math.Abs(value - rt)) ? SubCell.RightBottom : SubCell.RightTop;

					break;

				case Edge.Bottom:

				default:

					res = (Math.Abs(value - lb) < Math.Abs(value - rb)) ? SubCell.LeftBottom : SubCell.RightBottom;

					break;

			}



			ValuesInCell subValues = vc.GetSubCell(res);

			bool valueInside = subValues.ValueBelongTo(value);

			if (!valueInside)

			{

				throw new IsolineGenerationException(Strings.Exceptions.IsolinesDataIsUndetailized);

			}



			return res;

		}



		private static Point GetPoint(double value, double a1, double a2, Vector v1, Vector v2)

		{

			double ratio = (value - a1) / (a2 - a1);



			Verify.IsTrue(0 <= ratio && ratio <= 1);



			Vector r = (1 - ratio) * v1 + ratio * v2;

			return new Point(r.X, r.Y);

		}



		private Point GetPointXY(Edge edge, double value, ValuesInCell vc, IrregularCell rect)

		{

			double lt = vc.LeftTop;

			double lb = vc.LeftBottom;

			double rb = vc.RightBottom;

			double rt = vc.RightTop;



			switch (edge)

			{

				case Edge.Left:

					return GetPoint(value, lb, lt, rect.LeftBottom, rect.LeftTop);

				case Edge.Top:

					return GetPoint(value, lt, rt, rect.LeftTop, rect.RightTop);

				case Edge.Right:

					return GetPoint(value, rb, rt, rect.RightBottom, rect.RightTop);

				case Edge.Bottom:

					return GetPoint(value, lb, rb, rect.LeftBottom, rect.RightBottom);

				default:

					throw new InvalidOperationException();

			}

		}



		private bool BelongsToEdge(double value, double edgeValue1, double edgeValue2, bool onBoundary)

		{

			if (!Double.IsNaN(missingValue) && (edgeValue1 == missingValue || edgeValue2 == missingValue))

				return false;



			if (onBoundary)

			{

				return (edgeValue1 <= value && value < edgeValue2) ||

				(edgeValue2 <= value && value < edgeValue1);

			}

			else

			{

				return (edgeValue1 < value && value < edgeValue2) ||

					(edgeValue2 < value && value < edgeValue1);

			}

		}



		private bool IsPassed(Edge edge, int i, int j, byte[,] edges)

		{

			switch (edge)

			{

				case Edge.Left:

					return (i == 0) || (edges[i, j] & (byte)edge) != 0;

				case Edge.Bottom:

					return (j == 0) || (edges[i, j] & (byte)edge) != 0;

				case Edge.Top:

					return (j == edges.GetLength(1) - 2) || (edges[i, j + 1] & (byte)Edge.Bottom) != 0;

				case Edge.Right:

					return (i == edges.GetLength(0) - 2) || (edges[i + 1, j] & (byte)Edge.Left) != 0;

				default:

					throw new InvalidOperationException();

			}

		}



		private void MakeEdgePassed(Edge edge, int i, int j)

		{

			switch (edge)

			{

				case Edge.Left:

				case Edge.Bottom:

					edges[i, j] |= (byte)edge;

					break;

				case Edge.Top:

					edges[i, j + 1] |= (byte)Edge.Bottom;

					break;

				case Edge.Right:

					edges[i + 1, j] |= (byte)Edge.Left;

					break;

				default:

					throw new InvalidOperationException();

			}

		}



		private Edge TrackLine(Edge inEdge, double value, ref int x, ref int y, out double newX, out double newY)

		{

			// Getting output edge

			ValuesInCell vc = (missingValue.IsNaN()) ?

				(new ValuesInCell(values[x, y],

					values[x + 1, y],

					values[x + 1, y + 1],

					values[x, y + 1])) :

				(new ValuesInCell(values[x, y],

					values[x + 1, y],

					values[x + 1, y + 1],

					values[x, y + 1],

					missingValue));



			IrregularCell rect = new IrregularCell(

				grid[x, y],

				grid[x + 1, y],

				grid[x + 1, y + 1],

				grid[x, y + 1]);



			Edge outEdge = GetOutEdge(inEdge, vc, rect, value);

			if (outEdge == Edge.None)

			{

				newX = newY = -1; // Impossible cell indices

				return Edge.None;

			}



			// Drawing new segment

			Point point = GetPointXY(outEdge, value, vc, rect);

			newX = point.X;

			newY = point.Y;

			segments.AddPoint(point);

			processed[x, y] = true;



			// Whether out-edge already was passed?

			if (IsPassed(outEdge, x, y, edges)) // line is closed

			{

				//MakeEdgePassed(outEdge, x, y); // boundaries should be marked as passed too

				return Edge.None;

			}



			// Make this edge passed

			MakeEdgePassed(outEdge, x, y);



			// Getting next cell's indices

			switch (outEdge)

			{

				case Edge.Left:

					x--;

					return Edge.Right;

				case Edge.Top:

					y++;

					return Edge.Bottom;

				case Edge.Right:

					x++;

					return Edge.Left;

				case Edge.Bottom:

					y--;

					return Edge.Top;

				default:

					throw new InvalidOperationException();

			}

		}



		private void TrackLineNonRecursive(Edge inEdge, double value, int x, int y)

		{

			int s = x, t = y;



			ValuesInCell vc = (missingValue.IsNaN()) ?

				(new ValuesInCell(values[x, y],

					values[x + 1, y],

					values[x + 1, y + 1],

					values[x, y + 1])) :

				(new ValuesInCell(values[x, y],

					values[x + 1, y],

					values[x + 1, y + 1],

					values[x, y + 1],

					missingValue));



			IrregularCell rect = new IrregularCell(

				grid[x, y],

				grid[x + 1, y],

				grid[x + 1, y + 1],

				grid[x, y + 1]);



			Point point = GetPointXY(inEdge, value, vc, rect);



			segments.StartLine(point, (value - minMax.Min) / (minMax.Max - minMax.Min), value);



			MakeEdgePassed(inEdge, x, y);



			//processed[x, y] = true;



			double x2, y2;

			do

			{

				inEdge = TrackLine(inEdge, value, ref s, ref t, out x2, out y2);

			} while (inEdge != Edge.None);

		}



		#endregion



		private bool HasIsoline(int x, int y)

		{

			return (edges[x, y] != 0 &&

				((x < edges.GetLength(0) - 1 && edges[x + 1, y] != 0) ||

				 (y < edges.GetLength(1) - 1 && edges[x, y + 1] != 0)));

		}



		/// <summary>Finds isoline for specified reference value</summary>

		/// <param name="value">Reference value</param>

		private void PrepareCells(double value)

		{

			double currentRatio = (value - minMax.Min) / (minMax.Max - minMax.Min);



			if (currentRatio < 0 || currentRatio > 1)

				return; // No contour lines for such value



			int xSize = dataSource.Width;

			int ySize = dataSource.Height;

			int x, y;

			for (x = 0; x < xSize; x++)

				for (y = 0; y < ySize; y++)

					edges[x, y] = 0;



			processed = new bool[xSize, ySize];



			// Looking in boundaries.

			// left

			for (y = 1; y < ySize; y++)

			{

				if (BelongsToEdge(value, values[0, y - 1], values[0, y], true) &&

					(edges[0, y - 1] & (byte)Edge.Left) == 0)

				{

					TrackLineNonRecursive(Edge.Left, value, 0, y - 1);

				}

			}



			// bottom

			for (x = 0; x < xSize - 1; x++)

			{

				if (BelongsToEdge(value, values[x, 0], values[x + 1, 0], true)

					&& (edges[x, 0] & (byte)Edge.Bottom) == 0)

				{

					TrackLineNonRecursive(Edge.Bottom, value, x, 0);

				};

			}



			// right

			x = xSize - 1;

			for (y = 1; y < ySize; y++)

			{

				// Is this correct?

				//if (BelongsToEdge(value, values[0, y - 1], values[0, y], true) &&

				//    (edges[0, y - 1] & (byte)Edge.Left) == 0)

				//{

				//    TrackLineNonRecursive(Edge.Left, value, 0, y - 1);

				//};



				if (BelongsToEdge(value, values[x, y - 1], values[x, y], true) &&

					(edges[x, y - 1] & (byte)Edge.Left) == 0)

				{

					TrackLineNonRecursive(Edge.Right, value, x - 1, y - 1);

				};

			}



			// horizontals

			for (x = 1; x < xSize - 1; x++)

				for (y = 1; y < ySize - 1; y++)

				{

					if ((edges[x, y] & (byte)Edge.Bottom) == 0 &&

						BelongsToEdge(value, values[x, y], values[x + 1, y], false) &&

						!processed[x, y - 1])

					{

						TrackLineNonRecursive(Edge.Top, value, x, y - 1);

					}

					if ((edges[x, y] & (byte)Edge.Bottom) == 0 &&

						BelongsToEdge(value, values[x, y], values[x + 1, y], false) &&

						!processed[x, y])

					{

						TrackLineNonRecursive(Edge.Bottom, value, x, y);

					}

					if ((edges[x, y] & (byte)Edge.Left) == 0 &&

						BelongsToEdge(value, values[x, y], values[x, y - 1], false) &&

						!processed[x - 1, y - 1])

					{

						TrackLineNonRecursive(Edge.Right, value, x - 1, y - 1);

					}

					if ((edges[x, y] & (byte)Edge.Left) == 0 &&

						BelongsToEdge(value, values[x, y], values[x, y - 1], false) &&

						!processed[x, y - 1])

					{

						TrackLineNonRecursive(Edge.Left, value, x, y - 1);

					}

				}

		}



		/// <summary>

		/// Builds isoline data for 2d scalar field contained in data source.

		/// </summary>

		/// <returns>Collection of data describing built isolines.</returns>

		public IsolineCollection BuildIsoline()

		{

			VerifyDataSource();



			segments = new IsolineCollection();



			// Cannot draw isolines for fields with one dimension lesser than 2

			if (dataSource.Width < 2 || dataSource.Height < 2)

				return segments;



			Init();



			if (!minMax.IsEmpty)

			{

				values = dataSource.Data;

				double[] levels = GetLevelsForIsolines();



				foreach (double level in levels)

				{

					PrepareCells(level);

				}



				if (segments.Lines.Count > 0 && segments.Lines[segments.Lines.Count - 1].OtherPoints.Count == 0)

					segments.Lines.RemoveAt(segments.Lines.Count - 1);

			}

			return segments;

		}



		private void Init()

		{

			if (dataSource.Range.HasValue)

				minMax = dataSource.Range.Value;

			else

				minMax = (Double.IsNaN(missingValue) ? dataSource.GetMinMax() : dataSource.GetMinMax(missingValue));



			if (dataSource.MissingValue.HasValue)

				missingValue = dataSource.MissingValue.Value;



			segments.Min = minMax.Min;

			segments.Max = minMax.Max;

		}



		/// <summary>

		/// Builds isoline data for the specified level in 2d scalar field.

		/// </summary>

		/// <param name="level">The level.</param>

		/// <returns></returns>

		public IsolineCollection BuildIsoline(double level)

		{

			VerifyDataSource();



			segments = new IsolineCollection();



			Init();

			

			if (!minMax.IsEmpty)

			{

				values = dataSource.Data;





				PrepareCells(level);



				if (segments.Lines.Count > 0 && segments.Lines[segments.Lines.Count - 1].OtherPoints.Count == 0)

					segments.Lines.RemoveAt(segments.Lines.Count - 1);

			}

			return segments;

		}



		private void VerifyDataSource()

		{

			if (dataSource == null)

				throw new InvalidOperationException(Strings.Exceptions.IsolinesDataSourceShouldBeSet);

		}



		IsolineCollection segments;



		private double[,] values;

		private byte[,] edges;

		private Point[,] grid;



		private Range<double> minMax;

		private IDataSource2D<double> dataSource;

		/// <summary>

		/// Gets or sets the data source - 2d scalar field.

		/// </summary>

		/// <value>The data source.</value>

		public IDataSource2D<double> DataSource

		{

			get { return dataSource; }

			set

			{

				if (dataSource != value)

				{

					value.VerifyNotNull("value");



					dataSource = value;

					grid = dataSource.Grid;

					edges = new byte[dataSource.Width, dataSource.Height];

				}

			}

		}



		private const double shiftPercent = 0.05;

		private double[] GetLevelsForIsolines()

		{

			double[] levels;

			double min = minMax.Min;

			double max = minMax.Max;



			double step = (max - min) / (density - 1);

			double delta = (max - min);



			levels = new double[density];

			levels[0] = min + delta * shiftPercent;

			levels[levels.Length - 1] = max - delta * shiftPercent;



			for (int i = 1; i < levels.Length - 1; i++)

				levels[i] = min + i * step;



			return levels;

		}

	}

}