context: #region Copyright
  /*
  * Copyright (c) 2005-2010, OpenMI Association
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
  *     * Redistributions of source code must retain the above copyright
  *       notice, this list of conditions and the following disclaimer.
  *     * Redistributions in binary form must reproduce the above copyright
  *       notice, this list of conditions and the following disclaimer in the
  *       documentation and/or other materials provided with the distribution.
  *     * Neither the name of the OpenMI Association nor the
  *       names of its contributors may be used to endorse or promote products
  *       derived from this software without specific prior written permission.
  *
  * THIS SOFTWARE IS PROVIDED BY "OpenMI Association" ``AS IS'' AND ANY
  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  * DISCLAIMED. IN NO EVENT SHALL "OpenMI Association" BE LIABLE FOR ANY
  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
  #endregion
  #region Copyright
  
  ///////////////////////////////////////////////////////////
  //
  // namespace: org.OpenMI.Utilities.Spatial 
  // purpose: Utility for spatial data-operations within the LinkableComponents wrappers
  // file: XYGeometryTools.cs
  //
  ///////////////////////////////////////////////////////////
  //
  //    Copyright (C) 2006 OpenMI Association
  //
  //    This library is free software; you can redistribute it and/or
  //    modify it under the terms of the GNU Lesser General Public
  //    License as published by the Free Software Foundation; either
  //    version 2.1 of the License, or (at your option) any later version.
  //
  //    This library 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
  //    Lesser General Public License for more details.
  //
  //    You should have received a copy of the GNU Lesser General Public
  //    License along with this library; if not, write to the Free Software
  //    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  //    or look at URL www.gnu.org/licenses/lgpl.html
  //
  //    Contact info: 
  //      URL: www.openmi.org
  //	Email: sourcecode@openmi.org
  //	Discussion forum available at www.sourceforge.net
  //
  //      Coordinator: Roger Moore, CEH Wallingford, Wallingford, Oxon, UK
  //
  ///////////////////////////////////////////////////////////
  //
  //  Original author: Peter Sinding & Jan B. Gregersen, DHI - Water & Environment, Horsholm, Denmark
  //  Created on:      6 April 2005
  //  Version:         1.0.0 
  //
  //  Modification history:
  //  
  //
  ///////////////////////////////////////////////////////////
  
  #endregion
  
  using System;
  using System.Collections;
  
  namespace Oatc.UpwardsComp.Spatial
  {
  	/// <summary>
  	/// The XYGeometryTools class is a collection of general geometry functions. All functions are 
  	/// static methods that performs calculations on input given as parameters and returns a result.
  	/// 
  	/// The parameters passed to the XYGeometryTools methods are typically of type XYPoint, XYLine, 
  	/// XYPolyline or XYPolygon.
  	/// </summary>
  	public class XYGeometryTools
  	{
  		private const double EPSILON = 1e-5;
  
  		//===============================================================================================
  		//CalculatePointToPointDistance(XYPoint p1, XYPoint p2) : double
  		//===============================================================================================
  		/// <summary>
  		/// Returns the distance between the two points.
  		/// </summary>
  		/// <param name="p1">Point</param>
  		/// <param name="p2">Point</param>
  		/// <returns>Point to point distance</returns>
  		public static double CalculatePointToPointDistance(XYPoint p1, XYPoint p2)
  		{
  			return Math.Sqrt((p1.X - p2.X)*(p1.X - p2.X) + (p1.Y - p2.Y)*(p1.Y - p2.Y));
  		}
  
  		//===============================================================================================
  		//DoLineSegmentsIntersect(double x1, double y1, double x2,  double y2, double x3,	double y3, double x4, double y4) : bool
  		//===============================================================================================
  		/// <summary>
  		/// Returns true if two line segments intersects. The lines are said to intersect if the lines
  		/// axctually crosses and not if they only share a point. 
  		/// </summary>
  		/// <param name="x1">x-coordiante for first point in first line segment</param>
  		/// <param name="y1">y-coordinate for first point in first line segment </param>
  		/// <param name="x2">x-cooedinate for second point in first line segment</param>
  		/// <param name="y2">y-coordinate for second point in first line segment</param>
  		/// <param name="x3">x-coordinate for the first point in second line segment</param>
  		/// <param name="y3">y-coordinate for the first point in second line segment</param>
  		/// <param name="x4">x-coordinate for the second point in the second line segment</param>
  		/// <param name="y4">y-coordinate for the second point in the second line segment</param>
  		/// <returns>True if the line segments intersects otherwise false.</returns>
  		public static bool DoLineSegmentsIntersect(double x1, double y1, double x2, double y2, double x3, double y3, double x4,
  		                                           double y4)
  		{
  			double detP1P2P3, detP1P2P4, detP3P4P1, detP3P4P2;
  			bool intersect = false;
  
  			detP1P2P3 = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
  			detP1P2P4 = (x2 - x1)*(y4 - y1) - (x4 - x1)*(y2 - y1);
  			detP3P4P1 = (x3 - x1)*(y4 - y1) - (x4 - x1)*(y3 - y1);
  			detP3P4P2 = detP1P2P3 - detP1P2P4 + detP3P4P1;
  
  			if ((detP1P2P3*detP1P2P4 < 0) && (detP3P4P1*detP3P4P2 < 0))
  			{
  				intersect = true;
  			}
  			return intersect;
  		}
  
  		//===============================================================================================
  		//DoLineSegmentsIntersect(XYPoint p1, XYPoint p2, XYPoint p3, XYPoint p4) : bool
  		//===============================================================================================
  		/// <summary>
  		/// OverLoad of DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4)
  		/// </summary>
  		/// <param name="p1">First point in first line</param>
  		/// <param name="p2">Second point in first line</param>
  		/// <param name="p3">First point in second line</param>
  		/// <param name="p4">Second point in second line</param>
  		/// <returns>true if the line segmenst intersects otherwise false</returns>
  		public static bool DoLineSegmentsIntersect(XYPoint p1, XYPoint p2, XYPoint p3, XYPoint p4)
  		{
  			return DoLineSegmentsIntersect(p1.X, p1.Y, p2.X, p2.Y, p3.X, p3.Y, p4.X, p4.Y);
  		}
  
  		//===============================================================================================
  		//DoLineSegmentsIntersect(XYLine line1, XYLine line2) : bool
  		//===============================================================================================
  		/// <summary>
  		/// OverLoad of DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4)
  		/// </summary>
  		/// <param name="line1">First line</param>
  		/// <param name="line2">Second line</param>
  		/// <returns>true if the line segmenst intersects otherwise false</returns>
  		public static bool DoLineSegmentsIntersect(XYLine line1, XYLine line2)
  		{
  			return
  				DoLineSegmentsIntersect(line1.P1.X, line1.P1.Y, line1.P2.X, line1.P2.Y, line2.P1.X, line2.P1.Y, line2.P2.X,
  				                        line2.P2.Y);
  		}
  
  		//=========================================================================================
  		//CalculateIntersectionPoint(XYPoint p1, XYPoint p2, XYPoint p3, XYPoint p4) : XYPoint
  		//=========================================================================================
  		/// <summary>
  		/// Calculate intersection point between two line segments.
  		/// </summary>
  		/// <param name="p1">First point in first line</param>
  		/// <param name="p2">Second point in first line</param>
  		/// <param name="p3">First point in second line</param>
  		/// <param name="p4">Second point in second line</param>
  		/// <returns>Intersection point</returns>
  		public static XYPoint CalculateIntersectionPoint(XYPoint p1, XYPoint p2, XYPoint p3, XYPoint p4)
  		{
  			if (!DoLineSegmentsIntersect(p1, p2, p3, p4))
  			{
  				throw new Exception(
  					"Attempt to calculate intersection point between non intersecting lines. CalculateIntersectionPoint failed.");
  			}
  
  			XYPoint interSectionPoint = new XYPoint();
  
  			double a = p1.X*p2.Y - p2.X*p1.Y;
  			double b = p3.X*p4.Y - p4.X*p3.Y;
  			double c = (p1.X - p2.X)*(p3.Y - p4.Y) - (p3.X - p4.X)*(p1.Y - p2.Y);
  
  			interSectionPoint.X = (a*(p3.X - p4.X) - (b*(p1.X - p2.X)))/c;
  			interSectionPoint.Y = (a*(p3.Y - p4.Y) - (b*(p1.Y - p2.Y)))/c;
  
  			return interSectionPoint;
  		}
  
  		//=========================================================================================
  		//CalculateIntersectionPoint(XYLine line1, XYLine line2) : XYPoint
  		//=========================================================================================
  		/// <summary>
  		/// OverLoad of CalculateIntersectionPoint(XYPoint p1, XYPoint p2, XYPoint p3, XYPoint p4)
  		/// </summary>
  		/// <param name="line1">First line</param>
  		/// <param name="line2">Second line</param>
  		/// <returns>Intersection point</returns>
  		public static XYPoint CalculateIntersectionPoint(XYLine line1, XYLine line2)
  		{
  			return CalculateIntersectionPoint(line1.P1, line1.P2, line2.P1, line2.P2);
  		}
  
  		//=========================================================================================
  		//CalculateLengthOfPolylineInsidePolygon(XYPolyline polyline, XYPolygon polygon) : double
  		//=========================================================================================
  		/// <summary>
  		/// Calculates the length of polyline inside polygon. Lines segments on the edges of 
  		/// polygons are included with half their length.
  		/// </summary>
  		/// <param name="polyline">Polyline</param>
  		/// <param name="polygon">Polygon</param>
  		/// <returns>
  		/// Length of polyline inside polygon.
  		/// </returns>
  		public static double CalculateLengthOfPolylineInsidePolygon(XYPolyline polyline, XYPolygon polygon)
  		{
  			double lengthInside = 0;
  			int numberOfLineSegments = polyline.Points.Count - 1;
  			for (int i = 0; i < numberOfLineSegments; i++)
  			{
  				XYLine line = new XYLine(polyline.GetLine(i));
  				lengthInside += CalculateLengthOfLineInsidePolygon(line, polygon);
  			}
  			return lengthInside;
  		}
  
  		//=========================================================================================
  		//CalculateSharedLength(XYLine lineA, XYLine lineB) : double
  		//=========================================================================================
  		/// <summary>
  		/// Calculates the length that two lines overlap.
  		/// </summary>
  		/// <param name="lineA">Line</param>
  		/// <param name="lineB">Line</param>
  		/// <returns>
  		/// Length of shared line segment.
  		/// </returns>
  		protected static double CalculateSharedLength(XYLine lineA, XYLine lineB)
  		{
  			if (Math.Abs(lineA.P2.X - lineA.P1.X) < EPSILON && Math.Abs(lineB.P2.X - lineB.P1.X) < EPSILON &&
  			    Math.Abs(lineA.P1.X - lineB.P1.X) < EPSILON)
  			{
  				double YP1A = Math.Min(lineA.P1.Y, lineA.P2.Y);
  				double YP2A = Math.Max(lineA.P1.Y, lineA.P2.Y);
  				double YP1B = Math.Min(lineB.P1.Y, lineB.P2.Y);
  
gt: 				double YP2B = Math.Max(lineB.P1.Y, lineB.P2.Y);