#region Written by Tamas Budavari (budavari@jhu.edu)
/**
 * SphericalLib.Halfspace
 * 
 * SphericalLib is a .NET library to deal with geometry on the unit sphere.
 * It is intented to be used for complex survey footprint descriptions in
 * standalone applications as well as in web services and perhaps Yukon.
 * 
 * Written by Tamas Budavari (budavari@jhu.edu)
 * See bottom of file for revision history
 * 
 * Current revision:
 *   ID:          $Id: Halfspace.cs,v 1.1.1.1 2006/11/14 02:56:39 gupchup Exp $
 *   Revision:    $Revision: 1.1.1.1 $
 *   Date:        $Date: 2006/11/14 02:56:39 $
 */
using System;
using System.Collections.Generic;
using System.Xml.Serialization;
#endregion

namespace Spherical
{
	/// <summary>
	/// Struct to hold a Halfspace that consists of a 
	/// center point and the cosine of the angle.
    /// <remarks>It is a mutable value type.</remarks>
	/// </summary>
    [Serializable, XmlRoot(Namespace = Const.NameSpace, IsNullable=false)]
    [CLSCompliant(true)]
	public struct Halfspace
    {
        #region Private members

        private Cartesian vector;
        private double cos0;

        #endregion
        #region Public properties

        public Cartesian Vector
        {
            get { return vector; }
            set { vector = value; }
        }

        [XmlAttribute]
        public double Offset
        {
            get { return cos0; }
            set { cos0 = value; }
        }

        [XmlIgnore]
        public double Cos0
        {
            get { return cos0; }
			set { cos0 = value; }
        }

        /// <summary>
        /// Extended sign property
        /// </summary>
        public ESign ESign
        {
            get
            {
                if (this.cos0 > Const.Epsilon) return ESign.Positive;
                else if (this.cos0 < -Const.Epsilon) return ESign.Negative;
                else return ESign.Zero;
            }
        }

        /// <summary>
        /// Radius of circle in radians
        /// </summary>
        public double RadiusInRadian
        {
            get
            {
                double rad;
                if (cos0 >= 1) rad = 0;
                else if (cos0 <= -1) rad = Math.PI;
                else rad = Math.Acos(cos0);
                return rad;
            }
        }

        #endregion
        #region Constructors of various kinds

        /// <summary>
        /// Create halfspace from xyz and cosine of angle.
        /// Normalize if requested.
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="z"></param>
        /// <param name="normalize"></param>
        /// <param name="cos0"></param>
		public Halfspace(double x, double y, double z, bool normalize, double cos0)
		{
            this.vector = new Cartesian(x, y, z, normalize);
			this.cos0 = cos0;
		}

		/// <summary>
		/// Create Halfspace from a center point and cosine of angle
		/// </summary>
		/// <param name="center"></param>
		/// <param name="d"></param>
		public Halfspace(Cartesian center, double cos0)
		{
            this.vector = center;
			this.cos0 = cos0;
		}

		/// <summary>
		/// Create Halfspace from center (RA,Dec) and radius in arcmin
		/// </summary>
		/// <param name="ra"></param>
		/// <param name="dec"></param>
		/// <param name="arcmin"></param>
		public Halfspace(double ra, double dec, double arcmin)
		{
			this.vector = new Cartesian(ra,dec);
			this.cos0 = Math.Cos(arcmin*Const.Arcmin2Radian);
		}

        /// <summary>
        /// Smallest circle that goes through 2 points
        /// </summary>
        /// <param name="p0"></param>
        /// <param name="p1"></param>
        public Halfspace(Cartesian p0, Cartesian p1)
        {
            if (p0.Equals(p1))
                throw new ArgumentException
                    ("Halfspace.ctor(p0,p1): same points");
            
            if (p0.Same(p1,Const.Tolerance))
                throw new ArgumentException
                    ("Halfspace.ctor(p0,p1): same points within default tolerance");

            Cartesian c = p0.MiddlePoint(p1, false);
            double norm = c.Normalize();
            if (double.IsNaN(norm))
            {
                norm = Math.Sqrt(p0.X * p0.X + p0.Y * p0.Y);
                if (!double.IsNaN(norm))
                {
                    this.vector = new Cartesian(p0.Y / norm, -p0.X / norm, 0, false);
                }
                else
                {
                    norm = Math.Sqrt(p0.Y * p0.Y + p0.Z * p0.Z);
                    this.vector = new Cartesian(0, -p0.Z / norm, p0.Z / norm, false);
                }
                this.cos0 = 0;
            }
            else
            {
                this.vector = c;
                this.cos0 = c.Dot(p0);
            }
        }

		/// <summary>
		/// Small circle constraint that goes through 3 points.
		/// </summary>
		/// <param name="p0"></param>
		/// <param name="p1"></param>
		/// <param name="p2"></param>
		public Halfspace(Cartesian p0, Cartesian p1, Cartesian p2)
		{
            if (p0 == p1) throw new ArgumentException
                    ("Halfspace.ctor(p0,p1,p2): points 0 & 1 are the same");
            if (p0 == p2) throw new ArgumentException
                    ("Halfspace.ctor(p0,p1,p2): points 0 & 2 are the same");
            if (p1 == p2) throw new ArgumentException
                    ("Halfspace.ctor(p0,p1,p2): points 1 & 2 are the same");

			Cartesian v1 = new Cartesian (p1.X-p0.X, p1.Y-p0.Y, p1.Z-p0.Z, false);
			Cartesian v2 = new Cartesian (p2.X-p0.X, p2.Y-p0.Y, p2.Z-p0.Z, false);
            this.vector = v1.Cross(v2, true);
			this.cos0 = this.vector.Dot(p0);
			if (this.cos0 < 0) 
				this.Invert();
		}

        #endregion
        #region Morphological operations

        /// <summary>
		/// Invert this halfspace (complementer on sphere)
		/// Operation: center *= -1 and cos0 *= -1
		/// </summary>
		public void Invert()
		{
            this.vector.Mirror();
            this.cos0 *= -1;
		}

		/// <summary>
		/// Return inverse (complementer) of this halfspace
		/// </summary>
		/// <returns></returns>
		public Halfspace Inverse()
		{
            return new Halfspace(-vector.X, -vector.Y, -vector.Z, false, -cos0);
		}

        /// <summary>
        /// Increase the angle by the given amount, and return this
        /// </summary>
        /// <param name="arcmin"></param>
        public Halfspace Grow(double arcmin)
        {
            double radius = Math.Acos(this.cos0) + arcmin * Const.Arcmin2Radian;
            if (radius > Math.PI) throw new ArgumentException
                ("Halfspace.Grow(): radius bigger than Pi");
            if (radius < 0) throw new ArgumentException
                ("Halfspace.Grow(): radius less than 0");
            this.cos0 = Math.Cos(radius);
            return this;
        }

		#endregion
		#region Geometrical relation to other objects

		/// <summary>
		/// Check if Halfspace covers a point (x,y,z) unit vector with some tolerance
		/// </summary>
		/// <param name="x"></param>
		/// <param name="y"></param>
		/// <param name="z"></param>
		/// <returns></returns>
		public bool Covers(double x, double y, double z, double tolerance)
		{
            double dot = x*vector.X + y*vector.Y + z*vector.Z;
            double lim = cos0 - tolerance;
            bool res = dot > lim;
            return res;
		}

		/// <summary>
		/// Check if Halfspace covers a point within given tolerance limit
		/// </summary>
		/// <param name="p"></param>
		/// <returns></returns>
		public bool Covers(Cartesian p, double tolerance)
		{
            return Covers(p.X, p.Y, p.Z, tolerance);
		}

		/// <summary>
		/// Check if this halfspace covers a list of points
		/// </summary>
		/// <param name="list"></param>
		/// <returns></returns>
		public bool Covers(IEnumerable<Cartesian> list, double tolerance)
		{
            foreach (Cartesian p in list)
                if (!this.Covers(p.X, p.Y, p.Z, tolerance))
                    return false;
			return true;
		}

		/// <summary>
		/// Count how many points are covered
		/// </summary>
		/// <param name="arr"></param>
		/// <returns></returns>
		public int CoversHowMany(IEnumerable<Cartesian> arr, double tolerance)
		{
			int n = 0;
            foreach (Cartesian p in arr)
                if (Covers(p.X, p.Y, p.Z, tolerance))
                    n++;
			return n;
		}

		/// <summary>
		/// Count how many points are excluded by this constraint
		/// </summary>
		/// <param name="arr"></param>
		/// <returns></returns>
		public int MasksHowMany(IEnumerable<Cartesian> arr, double tolerance)
		{
			int n = 0;
            foreach (Cartesian p in arr)
                if (!Covers(p.X, p.Y, p.Z, tolerance))
                    n++;
			return n;
		}

		/// <summary>
		/// Count how many roots are excluded by this constraint
		/// </summary>
		/// <param name="arr"></param>
		/// <returns></returns>
		internal int MasksHowMany(IEnumerable<Root> arr, double tolerance)
		{
			int n = 0;
            foreach (Root r in arr)
                if (!Covers(r.Point, tolerance))
                    n++;
			return n;
		}

		/// <summary>
		/// Check if Halfspace is identical to another within a given tolerance
		/// </summary>
		/// <param name="h"></param>
		/// <param name="tolerance"></param>
		/// <returns></returns>
		public bool Same(Halfspace h, double epsilon)
		{
			bool eq = false;
            if (Math.Abs(cos0 - h.cos0) < epsilon
                    && Math.Abs(vector.X - h.vector.X) < epsilon
                    && Math.Abs(vector.Y - h.vector.Y) < epsilon
                    && Math.Abs(vector.Z - h.vector.Z) < epsilon)
                eq = true;
			return eq;
		}

		/// <summary>
		/// Check if Halfspace covers another completely
		/// </summary>
		/// <param name="h"></param>
		/// <returns></returns>
		public bool Covers(Halfspace h, double tolerance)
		{
			double angle = vector.AngleInRadian(h.vector);
			return angle - tolerance < this.RadiusInRadian - h.RadiusInRadian;
		}

		/// <summary>
		/// Check if Halfspace is disjoint from another
		/// </summary>
		/// <param name="h"></param>
		/// <returns></returns>
		public bool Disjoint(Halfspace h, double tolerance)
		{
			double angle = Math.Acos(vector.Dot(h.vector));
			return angle > h.RadiusInRadian + this.RadiusInRadian - tolerance;
		}

        /// <summary>
        /// Point on circle to the west from the center
        /// </summary>
        public Cartesian PointWest()
        {
            double sin0 = Math.Sin(this.RadiusInRadian);
            double ra = vector.RA;
            double sinRA = Math.Sin(ra * Const.Degree2Radian);
            double cosRA = Math.Cos(ra * Const.Degree2Radian);
            double x = vector.X * cos0 + sinRA * sin0;
            double y = vector.Y * cos0 - cosRA * sin0;
            double z = vector.Z * cos0;
            Cartesian w = new Cartesian(x, y, z, false);
            return w;
        }

		#endregion
		#region Root finding

		/// <summary>
		/// Return a directed root or PointNaN if invalid. Call Roots() for better performance!
		/// </summary>
		/// <param name="h"></param>
		/// <returns></returns>
		public Cartesian Root (Halfspace h)
		{
			double u,v,w, g,d,b;
			
			g = vector.Dot(h.vector);
			if (g > 1 || g < -1) 
				return Cartesian.NaN;
			
			d = 1 - g*g;
			b = cos0*cos0 + h.cos0*h.cos0 - 2*g*cos0*h.cos0;			
			if (d < b) 
				return Cartesian.NaN;
			
			u = (cos0-h.cos0*g ) / d;
			v = (h.cos0-cos0*g ) / d;
			w = Math.Sqrt((d-b)/(d*d));

			Cartesian c = vector.Cross(h.vector,false);

			// calc point
			double x = u*vector.X + v*h.vector.X + w*c.X;
			double y = u*vector.Y + v*h.vector.Y + w*c.Y;
			double z = u*vector.Z + v*h.vector.Z + w*c.Z;
            Cartesian r = new Cartesian(x, y, z, true);
			
            return r;
		}

		/// <summary>
		/// Find both roots in one go and return a status flag, where negative 
        /// values signal a problem (points are set to Cartesian.NaN)
		/// </summary>
		/// <param name="h"></param>
		/// <param name="pos"></param>
		/// <param name="neg"></param>
		/// <returns></returns>
		public int Roots (Halfspace h, out Cartesian pos, out Cartesian neg)
		{
			double u,v,w, g,d,b;
			
			g = vector.Dot(h.vector);
			if (g > 1 || g < -1) 
			{
				pos = neg = Cartesian.NaN;
				return -1;
			}
			
			d = 1 - g*g;
			b = cos0*cos0 + h.cos0*h.cos0 - 2*g*cos0*h.cos0;			
			if (d < b) 
			{
				pos = neg = Cartesian.NaN;
				return -2;
			}

			u = (cos0-h.cos0*g ) / d;
			v = (h.cos0-cos0*g ) / d;
			w = Math.Sqrt((d-b)/(d*d));

			Cartesian c = vector.Cross(h.vector,false);

			// calc points
            double x, y, z;

			x = u*vector.X + v*h.vector.X + w*c.X;
			y = u*vector.Y + v*h.vector.Y + w*c.Y;
			z = u*vector.Z + v*h.vector.Z + w*c.Z;
            pos = new Cartesian(x, y, z, true);
            
			x = u*vector.X + v*h.vector.X - w*c.X;
			y = u*vector.Y + v*h.vector.Y - w*c.Y;
			z = u*vector.Z + v*h.vector.Z - w*c.Z;
            neg = new Cartesian(x, y, z, true);

			return 2;
		}

		#endregion
        #region Equality related methods...

        public static bool operator ==(Halfspace left, Halfspace right)
        {
            return (left.vector == right.vector) && (left.cos0 == right.cos0);
        }

        public static bool operator !=(Halfspace left, Halfspace right)
        {
            return (left.vector != right.vector) || (left.cos0 != right.cos0);
        }

        public override bool Equals(object right)
        {
            if (right == null)
                return false;

            if (object.ReferenceEquals(this, right))
                return true;

            if (this.GetType() != right.GetType())
                return false;

            Halfspace h = (Halfspace)right;
            return (vector == h.vector) && (cos0 == h.cos0);
        }


        public override int GetHashCode()
        {
            return base.GetHashCode();
        }

        #endregion
        #region Parsing, ToString, revision and radius comparison

        /// <summary>
		/// String representation is: x y z cos0
		/// </summary>
		/// <returns></returns>
		public override string ToString()
		{
			return String.Format("{0} {1:R}", vector, cos0);
		}
	
		/// <summary>
		/// Parse a string of "x y z d" into a Halfspace. 
		/// </summary>
		/// <param name="repr"></param>
		/// <returns></returns>
		public static Halfspace Parse(string repr, bool normalize)
		{
			char[] sep = new char[] {' '};

			// only spaces left after this (hopefully)
			string proc = repr.Trim().Replace(',',' ').Replace("\r\n"," ").Replace('\n',' ').Replace('\t',' ');
			
			// replace duplicate spaces with a single one
			string trim = "";
			while (trim != proc) 
			{
				trim = proc;
				proc = trim.Replace("  "," ");
			}

			string[] t = trim.Split(sep,4);
			double x = double.Parse(t[0]);
			double y = double.Parse(t[1]);
			double z = double.Parse(t[2]);
			double cos0 = double.Parse(t[3]);
            return new Halfspace(x, y, z, normalize, cos0);
        }

        /// <summary>
        /// Comparison by angle to be used for sorting
        /// </summary>
        /// <param name="h"></param>
        /// <param name="g"></param>
        /// <returns></returns>
        public static int ComparisonRadiusAscending(Halfspace h, Halfspace g)
        {
            return g.cos0.CompareTo(h.cos0);
        }

		/// <summary>
		/// Revision from CVS
		/// </summary>
		public static string Revision
		{
			get { return "$Revision: 1.1.1.1 $"; }
		}

        #endregion
    }
}

#region Revision History
/* Revision History

        $Log: Halfspace.cs,v $
        Revision 1.1.1.1  2006/11/14 02:56:39  gupchup
        no message

        Revision 1.10  2006/05/22 18:55:54  budavari
        BinaryFormatter now works with Region and co

        Revision 1.9  2006/04/19 16:51:01  budavari
        EVerything is now CLS compliant

        Revision 1.8  2006/04/10 19:58:26  budavari
        Explicit not nullable attrib added

        Revision 1.7  2006/03/09 20:30:07  budavari
        Equality checks for value types

        Revision 1.6  2006/03/07 14:54:37  budavari
        Added set to public properties...

        Revision 1.5  2006/03/03 21:22:16  gyuri
        Made Cos0 settable not just gettable

        Revision 1.4  2006/02/28 16:25:28  budavari
        Major revision: struct for speed

        Revision 1.3  2006/01/30 19:09:50  budavari
        MinimalEnclosingCircle w/ 3 points and Halfspace.ctor w/ 2 points

        Revision 1.2  2006/01/25 19:06:14  budavari
        Fixed serialization attribute

        Revision 1.1  2006/01/24 15:56:05  budavari
        Re-structured directory layout to accommodate more projects

        Revision 1.8  2006/01/20 15:19:13  budavari
        ESign.Negative bug fixed (-epsilon)

        Revision 1.7  2006/01/19 18:09:38  budavari
        Comparison for sorting by ascending radius

        Revision 1.6  2006/01/19 17:27:10  budavari
        Extended sign enum added

        Revision 1.5  2006/01/09 21:11:02  gyuri
        Added an extra space at the end of ToString() to make it work better for output

        Revision 1.4  2005/12/29 14:23:47  budavari
        Normalize in root finding

        Revision 1.3  2005/12/20 02:26:56  budavari
        New checks in Grow(), which now also returns itself

        Revision 1.2  2005/12/02 15:37:49  budavari
        Static Cartesian.TripleProduct() added

        Revision 1.1  2005/12/02 02:47:58  budavari
        HalfSpace to Halfspace

        Revision 1.16  2005/12/01 22:34:57  budavari
        Area struct added and fully functional

        Revision 1.15  2005/11/30 20:50:27  budavari
        Cosmetics... #region specifications etc added

        Revision 1.14  2005/11/30 20:01:53  budavari
        Dilate() renamed to Gro() per Alex's request

        Revision 1.13  2005/11/30 19:50:37  budavari
        Dilate() methods added. Also explicit tolerance required everywhere for clarity

        Revision 1.12  2005/11/30 16:03:53  budavari
        Doc added

        Revision 1.11  2005/11/30 15:32:37  budavari
        Header changes

        Revision 1.10  2005/11/24 16:27:16  budavari
        Region areas compared to SQL code and releated fixes of degenerate arcs.
        Looking good except maybe for one patch (convex 0) of the 'footprint',
        which has 2 degenerate triangles and patch area changes when shifting
        pivot point.

        Revision 1.9  2005/11/18 05:58:09  budavari
        Working on boolean operations on regions...

        Revision 1.8  2005/11/17 01:58:54  budavari
        Done with the convex stuff... Hurray!!! It can even deal with the "orange
        slices" to some extent: the problem is the area - coming up next...

        Revision 1.7  2005/11/16 18:56:59  budavari
        Only degenerate roots are missing... Also added sample convex static class.

        Revision 1.6  2005/11/15 22:20:20  budavari
        Using Generics - almost there

        Revision 1.5  2005/11/15 15:25:36  budavari
        Everything is a class, things are coming together for area computation, etc. Will moved to C# 2.0

        Revision 1.4  2005/11/05 17:47:44  budavari
        Fixed a bug in Convex caused by renaming things and did some tidying up...

        Revision 1.3  2005/11/04 22:49:36  budavari
        Re-write to have Point and Halfspace as struct again and other tinkerings. XML is still fine...

        Revision 1.2  2005/11/03 15:22:01  budavari
        Fairly major revision to make sure the objects serialize to XML nicely:
          - Added separate containers of Halfspaces and convexes
          - Customized XML output for compactness and friendly proxy interface
          - Implemented simple .Parse() methods taking cartesion representation

        Revision 1.1.1.1  2005/10/28 14:13:43  budavari
        Initial import of the SphericalLib library. It is pretty complete in terms
        of the classes but lacking the simplication and other advanced methods.
		
*/
#endregion