source: branches/HeuristicLab.Problems.GrammaticalOptimization/DynamicDataDisplay/Charts/Isolines/IsolineBuilder.cs @ 13777

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;
    }
  }
}