#region Written by Tamas Budavari (budavari@jhu.edu)
/**
 * SphericalLib.PatchFactory
 * 
 * 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: PatchFactory.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;
using System.Collections.Generic;
using System.Collections.ObjectModel;
#endregion

namespace Spherical
{
	/// <summary>
	/// PatchFactory has a set of roots and arclets plus related info
	/// </summary>
	internal class PatchFactory
    {
        #region Arclet struct

        /// <summary>
        /// Internally point to roots and halfspaces that make up an arc
        /// </summary>
        private struct Arclet
        {
            public int iCircle, iRoot1, iRoot2;
            public Arclet(int iCircle, int iRoot1, int iRoot2)
            {
                this.iCircle = iCircle;
                this.iRoot1 = iRoot1;
                this.iRoot2 = iRoot2;
            }
            public override string ToString()
            {
                return String.Format("{0} : {1} {2}", iCircle, iRoot1, iRoot2);
            }
        }

        #endregion
        #region Private members

        // convex and extra info about its halfspaces
		private Convex convex;
        private bool[] hasRoot;
        private Cartesian[] west;

        // list of roots and (alter)egos
        private List<Root> roots;
        private int[] ego;

        // list of arclets
        private List<Arclet> arclets;

        #endregion
        #region Public properties

        public ReadOnlyCollection<Root> Roots
        {
            get { return roots.AsReadOnly(); }
        }

        #endregion
        #region Constructors

		/// <summary>
		/// Derive roots of convex, save 'west point' of halfspaces, 
        /// merge roots and set their (alter)ego
		/// </summary>
		/// <returns></returns>
		public PatchFactory (Convex convex, double tolerance)
		{
            int num = convex.HalfspaceList.Count;
            this.west = new Cartesian[num];
            this.hasRoot = new bool[num];
            this.convex = convex;
            
            // derive roots
            this.roots = new List<Root>(num * (num - 1));
            Cartesian p1, p2;
			int nroots;
            // all the pairs of halfspaces
			for (int i=0; i<num-1; i++) 
			{
                west[i] = convex.HalfspaceList[i].PointWest();
				for (int j=i+1; j<num; j++)
				{
					nroots = convex.HalfspaceList[i].Roots 
                        (convex.HalfspaceList[j], out p1, out p2);
                    // # of roots is 2 - otherwise no roots
                    if (nroots == 2)
                    {
                        roots.Add(new Root(p1, convex, tolerance, i, j));
                        roots.Add(new Root(p2, convex, tolerance, j, i));
                        hasRoot[i] = true;
                        hasRoot[j] = true;
                    }
				}
			}
            // set the last west point as well
            west[num - 1] = convex.HalfspaceList[num - 1].PointWest();

            // add west point for those who have no roots
            for (int i = 0; i < hasRoot.Length; i++)
            {
                if (!hasRoot[i])
                    roots.Add(new Root(west[i], convex, tolerance, i, i));
            }
            
            // set ego of roots
            ego = new int[roots.Count];
            for (int i = 0; i < ego.Length; i++)
            {
                ego[i] = i;
            }

            // alter-ego (degeneracy)
            // set identical point for roots with the same location
            Cartesian Pi;
            Root Rj;
            for (int i = 0; i < roots.Count - 1; i++)
            {
                Pi = roots[i].Point;
                for (int j = i + 1; j < roots.Count; j++)
                {
                    Rj = roots[j];
                    if (Pi.Same(Rj.Point, tolerance))
                    {
                        Rj.Point = Pi;
                        ego[j] = ego[i];
                    }
                }
            }
            // setup empty arclet list
            arclets = new List<Arclet>(2*roots.Count);
        }
        #endregion
        #region Select this and that...

        /// <summary>
        /// Select roots with their status set to something
        /// </summary>
        /// <returns></returns>
        private List<int> SelectRoots(RootStatus status)
		{
			List<int> idx = new List<int>();
			for (int i=0; i<roots.Count; i++)
				if (status==RootStatus.Any || status==roots[i].Status)
					idx.Add(i);
			return idx;
		}

		/// <summary>
		/// Select roots with a certain 'good' flag value
		/// </summary>
		/// <param name="goodness"></param>
		/// <returns></returns>
        private List<int> SelectRoots(int goodness)
		{
            List<int> idx = new List<int>();
            for (int i = 0; i < roots.Count; i++)
				if (goodness == roots[i].Good)	
					idx.Add(i);
			return idx;
		}

		/// <summary>
		/// Return indices of circles that have roots in specified set
		/// </summary>
		/// <returns></returns>
        private List<int> SelectCirclesWithRoots(List<int> iRoots)
		{
			List<int> res = new List<int>();
            int p;
			foreach (int i in iRoots)
			{
                p = roots[i].Parent1;
                if (!res.Contains(p)) res.Add(p);
                p = roots[i].Parent2;
                if (!res.Contains(p)) res.Add(p);
			}				
			return res;
		}

		/// <summary>
		/// Return indices of roots on a circle specified by halfspace index.
		/// The list is sorted such that the roots have counterclockwise ordering.
		/// </summary>
		/// <param name="c"></param>
		/// <returns></returns>
        private int[] SelectRootsOnCircle(int c)
		{
			List<int> r = new List<int>();
			// loop over roots and include if circle 'c' is a parent
			for (int i=0; i<roots.Count; i++)
				if (c==roots[i].Parent1 || c==roots[i].Parent2)
					r.Add(i);
		
			// projection to halfspace's plane
			Halfspace h = convex.HalfspaceList[c];
            Cartesian w = west[c].Sub(h.Vector.Scaled(h.Cos0));
            w.Normalize(); // west in tangent plane
            Cartesian u = w.Cross(h.Vector, true); // north in tangent plane

			// calculate winding angle
			double[] phi = new double[r.Count];
			for (int j=0; j<r.Count; j++)
			{
				Cartesian x = roots[r[j]].Point;  // root on the sky
				phi[j] = Math.Atan2 (w.Dot(x), u.Dot(x));
			}

            // order copy of 'r' by increasing phi
            int[] ir = r.ToArray();
            Array.Sort(phi, ir);
            
			return ir;
		}

        /// <summary>
		/// Return all arclets of all circles
		/// </summary>
		/// <returns></returns>
        private void BuildArclets()
		{
			for (int c=0; c<this.convex.HalfspaceList.Count; c++)
			{
				// ordered roots of circle 'c'
				int[] idx = this.SelectRootsOnCircle(c);
				if (idx.Length==0) continue;
				// loop on pairs
				for (int k=1; k<idx.Length; k++)
				{
					arclets.Add (new Arclet(c, idx[k-1], idx[k]));
				}
				// last with first
				arclets.Add (new Arclet(c, idx[idx.Length-1],idx[0]));
			}
		}

        /// <summary>
        /// Return indices of visible arclets on given circles
        /// </summary>
        /// <returns></returns>
        private List<int> SelectArcletsVisible(RootStatus status, double tolerance)
        {
            List<int> iVis = new List<int>();
            Cartesian test;
            for (int i = 0; i < arclets.Count; i++)
            {
                Arclet a = arclets[i];

                // if endpoints have different status than requested, continue
                if (status != RootStatus.Any)
                    if (roots[a.iRoot1].Status != status || roots[a.iRoot2].Status != status) 
                        continue;

                // if circle is not the relevant parent of the roots...
                if (a.iCircle != roots[a.iRoot1].Parent2) continue;
                if (a.iCircle != roots[a.iRoot2].Parent1) continue;

                // materialize arc
                Arc arc = this.MakeArc(a);

                // if two points are identical (roots merged in constructor)
                if (arc.Point1.Equals(arc.Point2))
                {
                    // if different endpoints -> infinitesimal length
                    if (a.iRoot1 != a.iRoot2)
                        continue;

                    test = arc.Point1;
                }
                else
                {
                    test = arc.Middle;
                }

                if (convex.Covers(test, tolerance))
                    iVis.Add(i);
            }
            return iVis;
        }

        /// <summary>
        /// Materialize an arc from an arclet
        /// </summary>
        /// <param name="arclet"></param>
        /// <returns></returns>
        private Arc MakeArc(Arclet arclet)
        {
            Arc arc = new Arc(this.convex.HalfspaceList[arclet.iCircle],
                    roots[arclet.iRoot1].Point, roots[arclet.iRoot2].Point);
            return arc;
        }

		#endregion
		#region Masking circles...

		/// <summary>
		/// Initialize the 'good' flag of roots like this
		///		good =  1 : root is inside convex
		///		good =  0 : root's parents are 'goodcircles' (in S) but masked out by other 'goodcircles'
		///		good = -1 : root's parents are NOT both 'goodcircles'
		///		good = -2 : root's parents ARE 'goodcircles' (in S) but masked by at least one 'not goodcircle'
		/// </summary>
		/// <param name="S"></param>
		private void InitRootGoodness (List<int> S)
		{
			foreach (Root r in roots)
			{
                // if inside, then good
				if (r.Status == RootStatus.Inside) 
				{
					r.Good = 1;	
				}
                // if outside and parent's are NOT both goodcircles
				else if (!S.Contains(r.Parent1) || !S.Contains(r.Parent2)) 
				{
					r.Good = -1;
				}
                // if outside and both parents ARE goodcircles
				else 
				{
					// if there is anyone in S who already masks root
                    if (r.MaskedByExists(S.Contains))
						r.Good = 0;
					else
						r.Good = -2;
				}
			}
		}

		/// <summary>
		/// Return indices of halfspaces that are needed to represent the convex (still not minimal)
		/// and which one is visible.
		/// </summary>
        private void SelectCompleteCircles(double tolerance, int maxiter, 
            ref List<int> S, out List<int> Visible)
		{
			// everybody in S is visible
			Visible = new List<int>(S);
            
            int iter;
			for (iter=0; iter<maxiter; iter++)
			{
                // set 'good' flag of roots
                this.InitRootGoodness(S);

                // find halfspace that covers most roots with good=-2
				List<int> MinusTwo = SelectRoots(-2);
				if (MinusTwo.Count==0) break;
                Cartesian[] pts = Cartesian.SelectPoints(roots, MinusTwo);
                int nmax = -1, imax = 99999999, n;
				for (int i=0; i<convex.HalfspaceList.Count; i++)
				{
					if (S.Contains(i)) continue;
					// # of roots covered by halfspace 'i'
					n = convex.HalfspaceList[i].MasksHowMany(pts,tolerance);
					if (n > nmax)
					{
						nmax = n;
						imax = i;
					}
				}
				// throw up if nobody from the rest covers any of the good=-2 roots
				if (nmax==-1)
					throw new ApplicationException("SelectCompleteCircles(): nmax=-1 should never happen!");

                // add it to minimal set and invisible
				S.Add(imax);
			}

			// throw up if hit maxiter limit - return list is probably not right
			if (iter==maxiter)
				throw new ApplicationException
					(String.Format("SelectCompleteCircles(): exceeded maximum number of iterations ({0})",
                     maxiter));
        }

		#endregion
        #region Arclets and company...

        private List<int> SelectCirclesWithArclets(List<int> iGoodArclets)
        {
            List<int> iCirc = new List<int>(iGoodArclets.Count);
            int circ;
            foreach (int i in iGoodArclets)
            {
                circ = arclets[i].iCircle;
                if (!iCirc.Contains(circ))
                    iCirc.Add(circ);
            }
            return iCirc;
        }

        private List<int> SelectRootsInArclets(List<int> iArclets)
        {
            List<int> iRoot = new List<int>(iArclets.Count);
            int p;
            foreach (int i in iArclets)
            {
                p = arclets[i].iRoot1;
                if (!iRoot.Contains(p)) iRoot.Add(p);
                p = arclets[i].iRoot2;
                if (!iRoot.Contains(p)) iRoot.Add(p);
            }
            return iRoot;
        }

        /// <summary>
        /// Connect arclets into patch. Picks the first arclet index from the
        /// given list of arclets and return the indices of that make a patch.
        /// Also removes the used indices from the passed list of arclets.
        /// </summary>
        /// <param name="iArclet"></param>
        /// <returns></returns>
        private List<int> ConnectArclets (ref List<int> iArclet)
        {
            // first element of good arcs is where we start 
            // and its start point is where we want to arrive
            int first = iArclet[0];
            int first1 = arclets[first].iRoot1;

            // return array of arclets and remove the ones used
            List<int> iPatch = new List<int>(1);
            iPatch.Add(first);
            iArclet.RemoveAt(0);

            int prev2 = arclets[first].iRoot2;
            bool done = false;

            // done if its endpoint is the first
            if (ego[prev2] == ego[first1])
            {
                return iPatch;
            }

            // assemble a patch
            Arclet a;
            List<int> taken = new List<int>();
            while (!done)
            {
                bool found = false;
                foreach (int i in iArclet)
                {
                    a = arclets[i];
                    // if matching part, add
                    if (ego[a.iRoot1] == ego[prev2])
                    {
                        // if same circle as before p1, reject
                        if (roots[a.iRoot1].Parent1 == a.iCircle)
                        {
                            taken.Add(i);
                            found = true;
                            break;
                        }
                        iPatch.Add(i);
                        prev2 = arclets[i].iRoot2;
                        found = true;
                        taken.Add(i);
                        break;
                    }
                }
                if (!found)
                    throw new ApplicationException("Arclet.Connect(): no matching arclet found");

                if (ego[prev2] == ego[first1])
                    done = true;

                // remove 'taken' arcs
                foreach (int i in taken)
                    iArclet.Remove(i);
                taken.Clear();
            }

            return iPatch;
        }
        #endregion
        #region Derive patches

        /// <summary>
        /// Derive all patches of its convex and simplify it if requested
        /// </summary>
        /// <param name="maxiter"></param>
        /// <param name="simplify_convex"></param>
        /// <returns></returns>
        public List<Patch> DerivePatches(double tolerance, int maxiter, bool simplify_convex)
        {
            // make all possible arcs (with proper ordering)
            this.BuildArclets();

            // arclets with middle point inside convex
            List<int> iGoodArclet = this.SelectArcletsVisible(RootStatus.Inside, tolerance);

            // prune good circles/roots
            List<int> iGoodCirc = this.SelectCirclesWithArclets(iGoodArclet);
            //List<int> iGoodRoot = this.SelectRootsInArclets(iGoodArclet);

            // add halfspaces needed to mask outside roots
            List<int> iVisible;
            this.SelectCompleteCircles(tolerance, maxiter, ref iGoodCirc, out iVisible);

            // assemble patches one by one and update iGoodArclet list until it's empty
            List<Patch> patches = new List<Patch>();
            while (iGoodArclet.Count > 0)
            {
                List<int> iPatch = this.ConnectArclets(ref iGoodArclet);
                Arc[] arcs = new Arc[iPatch.Count];
                for (int j = 0; j < iPatch.Count; j++)
                {
                    int i = iPatch[j];
                    arcs[j] = this.MakeArc(arclets[i]);
                }
                patches.Add(new Patch(arcs,tolerance));
            }

            // simplify convex
            if (simplify_convex)
            {
                for (int i = this.convex.HalfspaceList.Count - 1; i >= 0; i--)
                    if (!iGoodCirc.Contains(i))
                        this.convex.RemoveAt(i);
            }

            return patches;
        }

        /// <summary>
        /// Test function prints info to stdout while deriving patches...
        /// </summary>
        public void Test()
        {
            Convex c = this.convex;
            Title(c.Name);
            Console.WriteLine("# of halfspaces: {0}", c.HalfspaceList.Count);
            Console.WriteLine(c);

            this.BuildArclets();
            Console.WriteLine(this);
            Console.WriteLine();

            // arclets with middle point inside convex (needed to assemble patches)
            List<int> iGoodArclet = this.SelectArcletsVisible(RootStatus.Inside, Const.Tolerance);

            // prune good circles/roots
            List<int> iGoodCirc = this.SelectCirclesWithArclets(iGoodArclet);
            List<int> iGoodRoot = this.SelectRootsInArclets(iGoodArclet);

            // add halfspaces needed to mask outside roots
            List<int> iVisible;
            this.SelectCompleteCircles(Const.Tolerance, Const.MaximumIteration, 
                ref iGoodCirc, out iVisible);

            // assemble a patch and update iGoodArclet list
            Console.WriteLine("Good arclets:");
            foreach (int i in iGoodArclet) Console.Out.Write(" {0}", i);
            Console.Out.WriteLine();

            List<int> iPatch;
            int counter = 0;
            while (iGoodArclet.Count > 0)
            {
                iPatch = this.ConnectArclets(ref iGoodArclet);
                Console.Out.WriteLine("Patch {0}", counter++);
                Arc[] arcs = new Arc[iPatch.Count];
                for (int j=0; j<arcs.Length; j++)
                {
                    int i = iPatch[j];
                    arcs[j] = this.MakeArc(arclets[i]);
                    Console.Out.Write("{0}\n{1}\n", i, arcs[j]);
                }
                Patch p = new Patch(arcs,Const.Tolerance);
                Console.Out.WriteLine("Patch area: {0}", p.Area);
                Console.Out.WriteLine("Patch bounding circle: {0}", p.Mec);
                Console.Out.WriteLine("      w/ center: {0} {1}", 
                    p.Mec.Vector.RA, p.Mec.Vector.Dec);
                Console.Out.WriteLine();
                Console.Out.WriteLine("# of more good arclets: {0}", iGoodArclet.Count);
            }

            // simplify convex
            for (int i = c.HalfspaceList.Count - 1; i >= 0; i--)
                if (!iGoodCirc.Contains(i))
                    c.RemoveAt(i);

            Console.WriteLine("# of halfspaces: {0}", c.HalfspaceList.Count);
            Console.WriteLine(c);
        }
        private static void Title(string title)
        { 
            string line = "---------------------------------------------------------------------------------------";
            Console.Out.WriteLine(line);
            Console.Out.WriteLine("          {0}", title);
            Console.Out.WriteLine(line);
        }

        #endregion
        #region ToString and revision

        /// <summary>
		/// String representation (roots and arclets)
		/// </summary>
		/// <returns></returns> 
		public override string ToString()
		{
			string r = "Roots:";
            for (int i = 0; i < roots.Count; i++)
                r += String.Format("\n\t{0}. {1} -> {2}", i, roots[i], ego[i]);
            if (arclets.Count > 0)
            {
                r += "\nArclets:";
                for (int i = 0; i < arclets.Count; i++)
                    r += String.Format("\n\t{0}. {1} ", i, arclets[i]);
            }
			return r;
		}

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

        #endregion
    }
}

#region Revision History
/* Revision History

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

        Revision 1.2  2006/02/28 16:26:59  budavari
        Major revision: public properties

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

        Revision 1.17  2005/12/29 14:27:05  budavari
        Bugfix for covering circles

        Revision 1.16  2005/12/20 02:22:17  budavari
        Bug fix for masking circles

        Revision 1.15  2005/12/04 16:50:55  budavari
        Small features and fixes. Also added Outline class - not function, yet

        Revision 1.14  2005/12/02 21:32:30  budavari
        member variables do start with small letter ;-) XML still the same

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

        Revision 1.12  2005/12/02 15:02:57  budavari
        RootStatus.Any and Const.Epsilon=1e-15 added

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

        Revision 1.10  2005/11/30 20:55:39  budavari
        Cosmetics... #region specifications etc added

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

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

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

        Revision 1.6  2005/11/29 15:19:21  budavari
        All public fields in Patch for now

        Revision 1.5  2005/11/29 14:00:37  budavari
        Arc visibility test now checks parents of roots, as well

        Revision 1.4  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.3  2005/11/18 05:58:09  budavari
        Working on boolean operations on regions...

        Revision 1.2  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.1  2005/11/16 19:12:37  budavari
        PatchFactory replaces RootArray - does most of the magic

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

        Revision 1.1  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

*/
#endregion