source: branches/1265_HeuristicLab.Visualization/HeuristicLab.Visualization/3.3/Primitives/PrimitiveUtil.cs @ 17717

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;

namespace HeuristicLab.Visualization {
  public static class PrimitiveUtil {
    public static List<PointD> ComputeIntersect(this Rectangle rect, Line line) {
      var a = rect.LowerLeft;
      var c = rect.UpperRight;
      var b = new PointD(c.X, a.Y);
      var d = new PointD(a.X, c.Y);
      var ls = line.Start;
      var le = line.End;

      var intersectionPoints = new List<PointD>();
      var p = PointD.Empty;
      if (LineIntersect(ls, le, a, b, ref p)) intersectionPoints.Add(p);
      if (LineIntersect(ls, le, b, c, ref p)) intersectionPoints.Add(p);
      if (LineIntersect(ls, le, c, d, ref p)) intersectionPoints.Add(p);
      if (LineIntersect(ls, le, d, a, ref p)) intersectionPoints.Add(p);

      // return the intersection points in the order of the distance to the start of the line
      intersectionPoints.Sort((p1, p2) => EuclideanDistance(p1, line.Start).CompareTo(EuclideanDistance(p2, line.Start)));

      return intersectionPoints;
    }

    public static List<PointD> ComputeIntersect(this Ellipse ell, Line line) {
      var intersectionPoints = new List<PointD>();

      var x1 = line.Start.X;
      var y1 = line.Start.Y;
      var x2 = line.End.X;
      var y2 = line.End.Y;
      var midX = ell.Center().X;
      var midY = ell.Center().Y;
      var h = ell.Size.Width / 2;
      var v = ell.Size.Height / 2;

      // translate coordinate system to be the center of the ellipse
      x1 -= midX;
      y1 -= midY;

      x2 -= midX;
      y2 -= midY;

      if (x1.IsAlmost(x2)) {
        var y = (v / h) * Math.Sqrt(h * h - x1 * x1);
        if (Math.Min(y1, y2) <= y && y <= Math.Max(y1, y2)) {
          intersectionPoints.Add(new PointD(x1 + midX, y + midY));
        }
        if (Math.Min(y1, y2) <= -y && -y <= Math.Max(y1, y2)) {
          intersectionPoints.Add(new PointD(x1 + midX, -y + midY));
        }
        return intersectionPoints;
      }

      var a = (y2 - y1) / (x2 - x1);
      var b = y1 - a * x1;
      var h2 = h * h;
      var v2 = v * v;
      var r = a * a * h2 + v2;
      var s = 2 * a * b * h2;
      var t = h2 * (b * b - v2);
      var d = s * s - 4 * r * t;

      if (d > 0) {
        var xi1 = (-s + Math.Sqrt(d)) / (2 * r);
        var xi2 = (-s - Math.Sqrt(d)) / (2 * r);

        var yi1 = a * xi1 + b;
        var yi2 = a * xi2 + b;

        if (isPointInLine(x1, x2, y1, y2, xi1, yi1)) {
          intersectionPoints.Add(new PointD(xi1 + midX, yi1 + midY));
        }
        if (isPointInLine(x1, x2, y1, y2, xi2, yi2)) {
          intersectionPoints.Add(new PointD(xi2 + midX, yi2 + midY));
        }
      } else if (d.IsAlmost(0)) {
        var xi = -s / (2 * r);
        var yi = a * xi + b;

        if (isPointInLine(x1, x2, y1, y2, xi, yi)) {
          intersectionPoints.Add(new PointD(xi + midX, yi + midY));
        }
      }
      intersectionPoints.Sort((p1, p2) => EuclideanDistance(p1, line.Start).CompareTo(EuclideanDistance(p2, line.Start)));
      return intersectionPoints;
    }

    public static PointD Center(this RectangularPrimitiveBase primitive) {
      return new PointD((primitive.LowerLeft.X + primitive.UpperRight.X) / 2, (primitive.LowerLeft.Y + primitive.UpperRight.Y) / 2);
    }

    #region lowlevel
    public static double EuclideanDistance(PointD a, PointD b) {
      var x = a.X - b.X;
      var y = a.Y - b.Y;
      return Math.Sqrt(x * x + y * y);
    }

    public static bool isPointInLine(double x1, double x2, double y1, double y2, double px, double py) {
      double xMin = Math.Min(x1, x2);
      double xMax = Math.Max(x1, x2);

      double yMin = Math.Min(y1, y2);
      double yMax = Math.Max(y1, y2);

      return xMin <= px && px <= xMax && yMin <= py && py <= yMax;
    }

    private static bool IsAlmost(this double x, double y, double eps = 1e-12) {
      return Math.Abs(x - y) < eps;
    }
    // finds the intersection point between line segments AB and CD
    private static bool LineIntersect(PointD a, PointD b, PointD c, PointD d, ref PointD intersectionPoint) {
      // false if either segment is zero-length
      if (a == b || c == d)
        return false;

      // false if the segments share a common endpoint
      if (a == c || b == c || a == d || b == d)
        return false;

      // translate the reference system so that point A is the new origin
      var bx = b.X - a.X;
      var by = b.Y - a.Y;
      var cx = c.X - a.X;
      var cy = c.Y - a.Y;
      var dx = d.X - a.X;
      var dy = d.Y - a.Y;

      var ab = Math.Sqrt(bx * bx + by * by);

      var sin = by / ab;
      var cos = bx / ab;

      // rotate the system so that point B is on the positive X axis
      var tmp = cx * cos + cy * sin;
      cy = cy * cos - cx * sin;
      cx = tmp;

      tmp = dx * cos + dy * sin;
      dy = dy * cos - dx * sin;
      dx = tmp;

      // return false if CD doesn't cross AB
      if (cy < 0 && dy < 0 || cy >= 0 && dy >= 0)
        return false;

      var xIntersect = dx + (cx - dx) * dy / (dy - cy);

      if (xIntersect < 0 || xIntersect > ab)
        return false;

      intersectionPoint = new PointD(a.X + xIntersect * cos, a.Y + xIntersect * sin);
      return true;
    }
    #endregion
  }
}