using System;
using System.Collections.Generic;
using System.Collections;
using System.Text;
using Spherical;

/*=====================================================================

  File:      Cover.cs
  Summary:   Creates HTMids or HTMid ranges from Regions
  Date:	     Jan 08, 2006
  Author:    Gyorgy Fekete

---------------------------------------------------------------------
  This spatial search library was developed by Alexander Szalay, Robert Brunner, 
  Peter Kunszt, and Gyorgy Fekete at the Department of Physics and Astronomy, 
  The Johns Hopkins University, in collaboration with Jim Gray, Microsoft Research.
  Details of the algorithms can be found at http://www.sdss.jhu.edu/htm/.
 
======================================================= */

namespace Spherical.HTM
{
	/// <summary>
	/// Generate a covermap for a region object
    /// Ways to call the static
    /// Cover.HidRange(reg);
    /// Cover.HidRange(reg, false, fullonly, dir);
    /// Cover.HidAugRange(reg, false, null);
	/// </summary>
	public class Cover {
		/// <summary>
		/// list of hid ranges are collected and manipulated here
		/// </summary>
		private Interval _range;
		/// <summary>
		/// The Region whose cover sill be computed
		/// </summary>
		private Region _reg;

        private int _deepestFullLevel;
        private int _deepestFullIndex;
        private int _fullIndex;
        private int _previousNodeLevel;
        private static  readonly int _minRanges_ = 12;
        private static  readonly int _maxRanges_ = 20;
		/// <summary>
		/// Cover will make list of level 20 HtmID ranges if this flag is false
		/// 
		/// If set, trixels are saved at their natural level.
		/// 
		/// </summary>
		bool _makeTrixels = false;
        /// <summary>
        /// If this flag is set, then fulls and partials are treated separately
        /// </summary>
        bool _fulleqpartial = true;
        bool _collectFulls = true;
        bool _collectPartials = true;

		/// <summary>
		/// used to limit recursion
		/// </summary>
		int _maxlevel;
		/// <summary>
		/// subdivide trixels to at least this level
		/// 
		/// even if th
		/// </summary>
		int _minlevel = 0;  
		/// <summary>
		/// the covermap hids or hid ranges go here
		/// </summary>
        /// savePair gets cleared exactly in one place, modified
        /// in exactly one place and returned in exactly one place
		List<Int64Pair> _savePairs;
        List<Int64AugPair> _saveAugPairs;
        List<Node> _AList;
        List<Node> _BList;
        List<Node> _CurrentList; // references to A/B above
        List<Node> _PreviousList;// references to A/B above
        List<Node> _FullList; // Once full, children of full need not be re-examined
		/// <summary>
		/// When set, tracing causes files to be created in the 
		/// _directory
		/// </summary>
		bool _dotrace = true;
		/// <summary>
		/// directory is the root of the trace text files
		/// 
		/// there are about 4 files generated here, if tracing is enabled
		/// </summary>
		string _directory = null;
		// 
		// These are invoked only when debugging the class
		// or when these four files are necessary
		// _tracefile gets the actual trace of the intersect algorithm
		// _savedfile gets the HIDS saved as a result
		// _halfspacesfile gets a description of the halfspaces
		// _extendedfile gets the level 20 ranges. isomorphic
		// to _savedfile
		/// <summary>
		/// trace object for algorithm, saves all steps
		/// </summary>
		private Trace _tracefile = null;
		/// <summary>
		/// trace object, saves all trixels
		/// </summary>
		private Trace _savedfile = null;
		/// <summary>
		/// trace object, saves all HID ranges
		/// </summary>
		private Trace _extendedfile = null;
		/// <summary>
		/// trace file, saves all halfspaces
		/// </summary>
		private Trace _halfspacesfile = null;
		/// <summary>
		/// mode of operation of this object controls whether or not traces are generated
		/// </summary>
        private int _nr;
        private int _nrOld;
        private int _nrDel;
        private Mode _mode = Mode.Normal;
        /// <summary>
        /// Compute the list of HtmIDs or level 20 HtmID ranges that cover the given region
        /// </summary>
        /// <param name="reg">The Region in question</param>
        /// <param name="varlen">if set HtmIDs are computed, ranges otherwise</param>
        /// <param name="directory">If not null, name of directory where algorithm trace information is written</param>
        /// <returns></returns>
        [CLSCompliant(true)]
        public static List<Int64Pair> HidRange(Region reg, bool varlen, bool fullonly, string directory) {

            List<Int64Pair> returnResult;
            Cover cover = new Cover(reg, varlen, directory);
            cover._collectFulls = true;
            cover._collectPartials = !fullonly;
            returnResult = cover.Xintersect();
            return returnResult;
        }
		/// <summary>
		/// Return a hid range list
		/// </summary>
		/// <param name="reg"></param>
		/// <returns></returns>
		[CLSCompliant(true)]
		public static List<Int64Pair> HidRange(Region reg) {
			return HidRange(reg, false, false, null);

		}
        /// <summary>
        /// Return a hid range list augmented with partial/full information
        /// </summary>
        /// <param name="reg"></param>
        /// <returns></returns>
        public static List<Int64AugPair> HidAugRange(Region reg, bool varlen, string directory) {

            List<Int64AugPair> returnResult;
            Cover cover = new Cover(reg, varlen, directory);
            returnResult = cover.Xintersect(false);
            return returnResult;
        }
		/// <summary>
		/// Return a hidlist of the covermap
		/// </summary>
		/// <param name="reg"></param>
		/// <returns>Only the "lo" component is significant in the list of pairs</returns>
		[CLSCompliant(true)]
		public static List<Int64Pair> HidList(Region reg) {
			return HidRange(reg, true, false, null);
		}

		/// <summary>
		/// Create an instance with optional tracing
		///
        /// <param name="reg">Input region</param>
        /// <param name="varlen">If set, save trixel HtmID's, and not ranges</param>
        /// <param name="dir">optional directory for algorithm tracing</param>
		/// </summary>
		private Cover(Region reg, bool varlen, string dir) {
			_range = new Interval();
			_savePairs = new List<Int64Pair>();
            _saveAugPairs = new List<Int64AugPair>();
            _AList = new List<Node>(50);
            _BList = new List<Node>(50);
            _FullList = new List<Node>(50);
			_dotrace = true;

			_reg = reg;
			_directory = dir;
			_makeTrixels = varlen;
            _previousNodeLevel = -1;

			if (_directory == null) {
				_mode = Mode.Normal;
			} else {
				_mode = Mode.File;
			}			
		}
        /// <summary>
        /// Save a trixel on the "where" list
        /// </summary>
        /// <param name="q"></param>
        /// <param name="reason"></param>
        /// <param name="where"></param>
        /// <param name="force">save without checking if it or ancestor is already on list</param>
        private void NodeSave(Node q, List<Node> where, bool force) {
            if(force) {
                where.Add(q);
            } else {
                foreach(Node n in where) {
                    if(Trixel.IsAncestor(n.hid, q.hid)) { // already in list, quit
                        return;
                    }
                }
                where.Add(q);
            }
        }
        private int NRangesInList(List<Node> where) {
            Interval it = new Interval();
            foreach(Node q in where) {
                Int64Pair range = Trixel.Extend(q.hid, 20);
                it.AddRange(range);
            }
            it.sort();
            it.compact();
            it.check();
            return it.pairList.Count; //
        }
        /// <summary>
        /// Compute angle between points A and B in the context of the given circle.
        /// 
        /// There are many ways to get to point B from A. The shortest distance is
        /// along a great circle segment, but sometimes the path is constrained to a small
        /// circle. This function computes the angular separation of points A and B
        /// constrained to a circle. Both A and B are presumed to be on the given
        /// circle, otherwise the return value may not be meaningful. No check is performed
        /// that ensures that A and B are in fact on the circel.
        /// 
        /// </summary>
        /// <param name="A">First point</param>
        /// <param name="B">Second point</param>
        /// <param name="circle">The circle expressed as a Halfspac</param>
        /// <returns>angles in degrees</returns>
        [CLSCompliant(true)]
   		public static double angleInContext(Cartesian A, Cartesian B, Halfspace circle){
			double ang;
			// double dot = A.Dot(B);
			bool negative = circle.Cos0 < 0.0;
			double D = Math.Abs(circle.Cos0);
			Cartesian axis;
			if (negative) {
				axis = new Cartesian(-circle.Vector.X, -circle.Vector.Y, -circle.Vector.Z, true);
			} else {
				axis = new Cartesian(circle.Vector.X, circle.Vector.Y, circle.Vector.Z, true);
			}
			// NEW: 2/20/2006 added D* ... 
			Cartesian Ap = new Cartesian(A.X - D* axis.X, A.Y - D * axis.Y, A.Z - D * axis.Z, true);
			Cartesian Bp = new Cartesian(B.X - D* axis.X, B.Y - D * axis.Y, B.Z - D * axis.Z, true);
			double dot = Ap.Dot(Bp);
			if (Math.Abs(-1.0 - dot) < Const.Tolerance) {
				ang = 180.0;
			} else if (Math.Abs(1.0 - dot) < Const.Tolerance) {
				ang = 0.0;
			} else {
				double radians;
				radians = Math.Acos(dot);
				ang = radians * Const.Radian2Degree;
				Cartesian AB = Ap.Cross(Bp, true);
				double sign = AB.Dot(axis);
				if (sign < 0.0 && !negative) {
					ang = 360.0 - ang;
				}
			}
			return ang;

		}

		/// <summary>
		/// See if triangle intersects patch or contains patch
		/// </summary>
		/// <param name="v0"></param>
		/// <param name="v1"></param>
		/// <param name="v2"></param>
		/// <param name="patch"></param>
		/// <returns></returns>
		private bool patchIntersect(Cartesian v0, Cartesian v1, Cartesian v2, Patch patch) {
			double width, th1, th2;
			int nroots;
			Cartesian root1, root2;
			Cartesian beg, end;
			Cartesian[] vertices = new Cartesian[3];
			Halfspace[] edges = new Halfspace[3];
			vertices[0] = v0;
			vertices[1] = v1;
			vertices[2] = v2;
			beg = vertices[2];

			edges[0] = new Halfspace(v2.Cross(v0, true), 0.0);
			edges[1] = new Halfspace(v0.Cross(v1, true), 0.0);
			edges[2] = new Halfspace(v1.Cross(v2, true), 0.0);

			if (edges[0].Covers(patch.ArcList[0].Point1, Const.Epsilon)) {
				if (edges[1].Covers(patch.ArcList[0].Point1, Const.Epsilon)) {
					if (edges[2].Covers(patch.ArcList[0].Point1, Const.Epsilon)) {
						return true;
					}
				}
			}
			for (int i = 0; i < 3; i++) {
				end = vertices[i];
				// Halfspace trixelEdge = new Halfspace(beg.Cross(end, true), 0.0);
				if (patch.ArcList.Count == 1) {
					// Just look at roots...
					nroots = edges[i].Roots(patch.ArcList[0].Circle, out root1, out root2);
					if (nroots > 0) { //check if root is on the trixel
						width = beg.AngleInDegree(end);
						th1 = angleInContext(beg, root1, edges[i]);
						th2 = angleInContext(root1, end, edges[i]);
						if (th1 + th2 <= width + Const.Tolerance) {
							return true;
						}
						th1 = angleInContext(beg, root2, edges[i]);
						th2 = angleInContext(root2, end, edges[i]);
						if (th1 + th2 <= width + Const.Tolerance) {
							return true;
						}
					}
				} else {
					foreach (Arc arc in patch.ArcList) { // does the arc intersect an edge? 
						nroots = edges[i].Roots(arc.Circle, out root1, out root2);
						// Edge from v0 to v1, and first arc in second patch:
						// +		root1	{0.989910532471754 -0.117043851642381 0.0798615958656336}
						// +		root2	{0.987586804864801 -0.150073961190551 0.0463692681580186}
						if (nroots > 0) {
							// See where the roots are on the arc first, then on the trixel
							// They must be on both
							// 
							// width = arc.length; // In degrees: 	31.115205244568092
							width = angleInContext(arc.Point1, arc.Point2, arc.Circle);
							th1 = angleInContext(arc.Point1, root1, arc.Circle);
							th2 = angleInContext(root1, arc.Point2, arc.Circle);
							if (th1 + th2 <= width + Const.Tolerance) {
								width = beg.AngleInDegree(end);
								th1 = angleInContext(beg, root1, edges[i]);
								th2 = angleInContext(root1, end, edges[i]);
								if (th1 + th2 <= width + Const.Tolerance) {
									return true;
								}
							}
							width = angleInContext(arc.Point1, arc.Point2, arc.Circle);
							th1 = angleInContext(arc.Point1, root2, arc.Circle);
							th2 = angleInContext(root2, arc.Point2, arc.Circle);
							if (th1 + th2 <= width + Const.Tolerance) {
								width = beg.AngleInDegree(end);
								th1 = angleInContext(beg, root2, edges[i]);
								th2 = angleInContext(root2, end, edges[i]);
								if (th1 + th2 <= width + Const.Tolerance) {
									return true;
								}
							}
						}
					}
				}
				beg = end;
			}
			// edges[0-2] already computed, see if a patchpoint is inside trixel
			//
			//if (Cartesian.TripleProduct(v1, v2, p) < -Const.Epsilon)
			//    return false;
			//if (Cartesian.TripleProduct(v2, v3, p) < -Const.Epsilon)
			//    return false;
			//if (Cartesian.TripleProduct(v3, v1, p) < -Const.Epsilon)
			//    return false;
			//return true;
			return false;
		}
		/// <summary>
		/// Either a patchpoint is the trixel or a path intersects the trixel
		/// </summary>
		/// <param name="v0"></param>
		/// <param name="v1"></param>
		/// <param name="v2"></param>
		/// <returns></returns>
		private bool testCombined(Cartesian v0, Cartesian v1, Cartesian v2){
			foreach (Convex c in _reg.ConvexList) {
				foreach (Patch p in c.PatchList) {
					if (patchIntersect(v0, v1, v2, p)) {
						return true;
					}
				}
			}
			return false;
		}

		//private bool testSomePatchPoint(Cartesian v0, Cartesian v1, Cartesian v2) {
		//    foreach (SuperConvex c in _sreg.convexList) {
		//        foreach (Patch p in c.patchList) {
		//            foreach (Arc a in p.arcList) {
		//                if (Trixel.Contains(a.p1, v0, v1, v2)){
		//                    return true;
		//                } else if (Trixel.Contains(a.p2, v0, v1, v2)) {
		//                    return true;
		//                }
		//            }
		//        }
		//    }
		//    return false;
		//}


		//private bool testAllPatch(Cartesian v0, Cartesian v1, Cartesian v2) {
		//    foreach (SuperConvex c in _sreg.convexList){
		//        foreach (Patch p in c.patchList){
		//            if (patchIntersect(v0, v1, v2, p)){
		//                return true;
		//            }
		//        }
		//    }
		//    return false;
		//}
		private Patch findPatch(Cartesian v, Halfspace bc,
			out bool reject) {
			bool vinsidep;
			bool alldisjoint = true; // see if at least one patch BC is not disjoint with trixel BC
			reject = true;
			foreach (Convex c in _reg.ConvexList){
				foreach (Patch p in c.PatchList){
					vinsidep = true; // Assume yes
					if (!p.Mec.Disjoint(bc, Const.Tolerance)){
						alldisjoint = false;
						reject = false;
						// is v inside the patch?
						foreach (Arc arc in p.ArcList){
							if (!arc.Circle.Covers(v, Const.Tolerance)){ // v is not in there!
								vinsidep = false;
								break; // move on to next patch
							}
						}
						if (vinsidep) { // inside this patch!
							return p;
						}
					}
					// Not inside this patch, keep looking
				}
			}
			reject = alldisjoint;
			return null;
		}
        private  Markup testNode(Node q) {
            return testNode(q.v0, q.v1, q.v2, q.hid);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="v0"></param>
        /// <param name="v1"></param>
        /// <param name="v2"></param>
        /// <param name="hid"></param>
        /// <returns></returns>
		private Markup testNode(
			Cartesian v0, Cartesian v1, Cartesian v2, Int64 hid) {
			//
			// Test this node against all patches, which are disjoint.
			// If BC of trixel is disjoint from BC of all patches, then clearly reject

			//
			// Find a patch that contains at least one vertex.
			//      Are all vertices is this patch?
			//                 Yes
			// Then see if trixel is full contained in this patch
			// If so, the trixel is 'FULL'
			// Else it is partial
			//
			// A vertex can be in at most one patch.
			Halfspace BC = new Halfspace(v0, v1, v2);
			Patch p0, p1, p2;
			bool alldisjoint;

			Markup mark = Markup.Dontknow;
			while (hid == 10L){ // 3085, 233907352L) {
				break;
			}

			//
			// Sets alldisjoint, if trixel BC is clear of all patch BC's
			p0 = findPatch(v0, BC, out alldisjoint);
			if (alldisjoint) {
				mark = Markup.Reject;
			} else {
				p1 = findPatch(v1, BC, out alldisjoint);
				if (alldisjoint)
					throw new ApplicationException("Internal error testNode:1");
				p2 = findPatch(v2, BC, out alldisjoint);
				if (alldisjoint)
					throw new ApplicationException("Internal error testNode:2");

				// If we get here, then not all BC's are disjoint...

				if (p0 != null) {
					if (p0 == p1 && p1 == p2) { // All three are inside same patch
						mark = Markup.Full;
					} else {
						mark = Markup.Partial;
					}
				} else { // p0 is null
					if (p1 == null && p2 == null) {
						// All three vertices of triangle are outside of patch
						// 1) is there a patch point inside? yes -> partial
						//

						if (testCombined(v0, v1, v2)) {
							mark = Markup.Partial;
						} 
						//// 1) is a patchpoint inside the triangle?
						//// 
						//if (testSomePatchPoint(v0, v1, v2)) {
						//    mark = Markup.Partial;
						//}
						////
						//// 2) is there an intersection?
						////
						//else if (testAllPatch(v0, v1, v2)) { // There was an intersection 
						//    mark = Markup.Partial;
						//} else {
						else {
							mark = Markup.Reject;
						}
					} else {
						mark = Markup.Partial;
					}
				}
			}
			if (_dotrace && _mode == Mode.File) {
				this._tracefile.dump("testHID {0} {1}", hid, mark);
			}
			if (mark == Markup.Dontknow) {
				throw new ApplicationException("Internal error testNode:3");
			}
			return mark;
		}
        internal enum Traversal {
            DepthFirst,
            BreadthFirst
        };
        internal class Node {
            public Int64 hid;
            public Cartesian v0;
            public Cartesian v1;
            public Cartesian v2;
            public bool terminal;
            public int level;
            public Savecode savecode;
            public override string ToString(){
                return Trixel.ToString(this.hid);
            }
            public Node(Int64 hid, Cartesian a, Cartesian b, Cartesian c, int level) {
                v0 = a;
                v1 = b;
                v2 = c;
                this.level = level;
                this.hid = hid;
                this.terminal = false;
            }
            /// <summary>
            /// Expand this node, put children on given queue
            /// </summary>
            /// <param name="queue"></param>
            public void Expand(Queue<Node> queue) {
                Node qc0, qc1, qc2, qc3;
                if(Trixel.IsValid(hid)) {
                    Int64 id0;
                    int newlevel = level + 1;
                    Cartesian w0 = v1.Add(v2); w0.Normalize();
                    Cartesian w1 = v2.Add(v0); w1.Normalize();
                    Cartesian w2 = v0.Add(v1); w2.Normalize();

                    id0 = hid << 2;
                    qc0 = new Node(id0  , v0, w2, w1, newlevel);
                    qc1 = new Node(id0+1, v1, w0, w2, newlevel);
                    qc2 = new Node(id0+2, v2, w1, w0, newlevel);
                    qc3 = new Node(id0+3, w0, w1, w2, newlevel);
                    queue.Enqueue(qc0);
                    queue.Enqueue(qc1);
                    queue.Enqueue(qc2);
                    queue.Enqueue(qc3);
                } 
            }
        }
        /// <summary>
        /// Creates and fills savePairs and saveTriples from supplied
        /// list of full and partial nodes
        /// </summary>
        /// <param name="partialNodes"></param>
        /// <param name="fullNodes"></param>
        private void SaveList(List<Node> partialNodes, List<Node> fullNodes) {
            //
            // Combine the two
            // for the range
            // but trace them separately
            // 
            Interval partials = new Interval();
            Interval fulls = new Interval();

            
            // ///////////////////////////////////////
            // Save Partials
            foreach(Node qt in partialNodes) {
                Int64Pair range = Trixel.Extend(qt.hid, 20);
                partials.AddRange(range);
            }
            partials.sort();
            partials.compact();
            partials.check();

            // ///////////////////////////////////////
            // Save Fulls
            foreach(Node qt in fullNodes) {
                Int64Pair range = Trixel.Extend(qt.hid, 20);
                fulls.AddRange(range);
            }
            fulls.sort();
            fulls.compact();
            fulls.check();
            //
            // INTERVALS fulls and partials must be compact and checked before continuing
            //
            if(_dotrace && _mode == Mode.File) {
                foreach(Node q in fullNodes) {
                    _savedfile.dump("{0,20} F", q.hid);
                    _tracefile.dump("Saving {0} F why {1} (evel {2})",
                        q.hid, Savecode.Full, q.level);
                }
                foreach(Node q in partialNodes) {
                    _savedfile.dump("{0,20} P", q.hid);
                    _tracefile.dump("Saving {0} P why {1} (evel {2})",
                        q.hid, Savecode.Partial, q.level);
                }
            }
            // ///////////////////////////////////////
            // COMBINE PARTIAL AND FULL into ALLIST
            if(_fulleqpartial) {
                List<Int64Pair> allList = Interval.combine(partials.pairList, fulls.pairList);

                Interval.compact(allList);
                Interval.IsWellFormed(allList);
                if(_makeTrixels) {
                    foreach(Node qt in fullNodes) {
                        _savePairs.Add(new Int64Pair(qt.hid, qt.hid));
                    }
                    foreach(Node qt in partialNodes) {
                        _savePairs.Add(new Int64Pair(qt.hid, qt.hid));
                    }

                } else {
                    foreach(Int64Pair lohi in allList) {
                        if(_dotrace && _mode == Mode.File) {
                            _extendedfile.dump("{0,20}  {1,20}", lohi.lo, lohi.hi);
                        }
                        _savePairs.Add(lohi);
                    }
                }
            } else { // Keep fulls and partials separate

                if(_makeTrixels) {
                    foreach(Node qt in fullNodes) {
                        _saveAugPairs.Add(new Int64AugPair(qt.hid, qt.hid, true));
                    }
                    foreach(Node qt in partialNodes) {
                        _saveAugPairs.Add(new Int64AugPair(qt.hid, qt.hid, false));
                    }

                } else {
                    ///
                    /// fulls and partials are compacted and checked...
                    /// 
                    foreach(Int64Pair lohi in fulls.pairList) {
                        if(_dotrace && _mode == Mode.File) {
                            _extendedfile.dump("{0,20}  {1,20}", lohi.lo, lohi.hi);
                        }
                        _saveAugPairs.Add(new Int64AugPair (lohi, true));
                    }
                    foreach(Int64Pair lohi in partials.pairList) {
                        if(_dotrace && _mode == Mode.File) {
                            _extendedfile.dump("{0,20}  {1,20}", lohi.lo, lohi.hi);
                        }
                        _saveAugPairs.Add(new Int64AugPair(lohi, false));
                    }

                }
            }
        }
        private void Visit(Node q) {
            //
            // compute markups, save,reject or partial
            // Node is terminal if any of the follolowing is true
            // (a) budget is exhausted
            // (b) level has been reached
            // We save partials and fulls on a temporary list
            // so that when we complete a level, we can examine the budget.
            // We shgould at least save two levels, so that if the Current level
            // is too pricey in terms of number of ranges, we can output the
            // previous level's range.
            Markup m;
            m = testNode(q);
            
            switch(m) {
                case Markup.Full:
                    if(_collectFulls) {
                        q.savecode = Savecode.Full;
                        NodeSave(q, _FullList, false);
                        _fullIndex++; // index of last slot taken
                        if(_deepestFullLevel < 0) {
                            _deepestFullLevel = q.level;
                        }
                        if(q.level > _deepestFullLevel) {
                            _deepestFullIndex = _fullIndex; // last slot taken is at a new level
                            _deepestFullLevel = q.level;
                        }
                        // Keep track of level transitions. you may have to delete the highest level
                        // nodes...

                        q.terminal = true;
                    }
                    break;
                case Markup.Reject:
                    q.terminal = true;
                    break;
                case Markup.Partial:
                    if(_collectPartials) {
                        q.savecode = Savecode.Partial;
                        NodeSave(q, _CurrentList, true);
                        // Let traverse decide whether or not its a terminal node
                    }
                    break;
            }
        }
        Queue<Node> theq = new Queue<Node>(100);

        private void traverse() {
            bool bailout = false;
            bool usecurrent = true;
            _CurrentList = _AList;
            _PreviousList = _BList;
            _CurrentList.Clear();
            _PreviousList.Clear(); // we like a fresh start
            _deepestFullLevel = -1;
            _deepestFullIndex = -1;
            _fullIndex = -1;
            //
            // _Current and _Previous simply pointers/references
            // to A and B list
            while(theq.Count > 0) {

                Node q = theq.Peek();
                if(q.level != _previousNodeLevel) {
                    // A new level is about to be explored.
                    // Before we rush into it, see if can stop.
                    // Current points to most recently finished level
                    // and previous is the one before it.
                    // First time special case, dont wait and see                        
                    if(_previousNodeLevel < 0) {
                        _previousNodeLevel = q.level; // Should be zero.
                        _nrOld = 0;
                    } else {
                        // Watch change in number of ranges (nr)
                        // (1) If the budget was exceeded, output stuff from previous
                        // level. 
                        // (2) If nr does not change, and maxlevel
                        // is exceeded, then stop (so yes, we do need maxlevel).
                        // output is in currentlist.
                        // (3)If the budget was exceeded by only a bit,
                        // output stuff from current level.
                        //
                        // Otherwise current becomes previous, and new current 
                        // should be cleared.
                        // CAREFUL!
                        // We really want the number of ranges after the current partial is combined
                        // with the list of full nodes.

                        Interval currents = new Interval();
                        Interval fulls = new Interval();
                        foreach(Node qt in _CurrentList) {
                            Int64Pair range = Trixel.Extend(qt.hid, 20);
                            currents.AddRange(range);
                        }
                        currents.sort();
                        currents.compact();
                        currents.check();

                        foreach(Node qt in _FullList) {
                            Int64Pair range = Trixel.Extend(qt.hid, 20);
                            fulls.AddRange(range);
                        }
                        fulls.sort();
                        fulls.compact();
                        fulls.check();
                        //
                        // Each level has its share of full and partial (currents)

                        List<Int64Pair> allList = Interval.combine(currents.pairList, fulls.pairList);
                        currents.Clear();
                        currents.pairList = allList;
                        currents.compact();
                        currents.check();
                        //
                        // Currents is now combined with 
                        // partial and full
                        //
                        _nr = currents.pairList.Count;
                        //
                        // was: _nr = NRangesInList(_CurrentList);
                        //
                        _nrDel = _nr - _nrOld;
                        if(_nr > _maxRanges_) { // (1)
                            theq.Clear();
                            bailout = true;
                            usecurrent = false; // USE PREVIOUS
                        } else if(q.level > _maxlevel) { // (2)
                            theq.Clear();
                            bailout = true;
                            usecurrent = true;// USE CURRENT
                        } else if(_nr > _minRanges_ && _nr <= _maxRanges_) { // (3)
                            theq.Clear();
                            bailout = true;
                            usecurrent = true;
                        } else {
                            List<Node> temp;
                            temp = _CurrentList;
                            _CurrentList = _PreviousList;
                            _PreviousList = temp;
                            _CurrentList.Clear();
                        }
                    }
                    //
                    //
                    _previousNodeLevel = q.level;
                } // if level changed
                //
                // NOTE: When bailout is set, theq is also cleared to prevent further
                //       iteration.
                if(!bailout){
                    theq.Dequeue();
                    Visit(q);

                    // visitation should build the binary tree of ranges
                    // can be sloppy, only playing with less than 1000 ranges
                    // The tree of intervals (TOI) can be traversed and cleared
                    // and inserted. Insertion also carries with it sorting and merging 
                    // intervals and nodes. For the future, now we hack the Intervals

                    if(!q.terminal) {
                        q.Expand(theq);
                    }
                }
            } // While
            if(usecurrent) {
                SaveList(_CurrentList, _FullList);
            } else {
                // See if Fullist is too large, like in the case of collecting only
                // fulls, and no partials...
                if(_FullList.Count > _maxRanges_) {
                    //
                    // trim it back...
                    // don't consider the highest level
                    // It will be organized like this:
                    // lower level nodes first, followed by higher levels.
                    // look from the back and find the level transition. delete the last
                    // batch of iso-level nodes.
                    // But we anticipated this, and kept track ...
                    //
                    _FullList.RemoveRange(_deepestFullIndex, _FullList.Count - _deepestFullIndex);
                }
                SaveList(_PreviousList, _FullList);
            }

        } // Method
        private void IntersectSetup() {
            _range.Clear();
            // ///////////////////////////////////////////////////
            double theta;
            double featuresize = 1.0;
            double a;
            double twopi = Math.PI * 2.0;
            double minArea = Const.WholeSphereInSquareDegree, maxArea = 0.0;
            foreach(Convex con in _reg.ConvexList) {
                foreach(Patch p in con.PatchList) {
                    a = p.Area;
                    if(a < 0.0 + Const.Tolerance) {
                        a += Const.WholeSphereInSquareDegree;
                    }
                    if(minArea > a) minArea = a;
                    if(maxArea < a) maxArea = a;
                }
            }
            a = maxArea / Const.SquareRadian2SquareDegree;
            theta = Math.Acos(1.0 - (a / twopi));
            featuresize = 2.0 * theta;
            _minlevel = HtmState.Instance.getLevel(featuresize);

            a = minArea / Const.SquareRadian2SquareDegree;
            theta = Math.Acos(1.0 - (a / twopi));
            featuresize = 2.0 * theta;
            _maxlevel = HtmState.Instance.getLevel(featuresize);

            if(_maxlevel < _minlevel + HtmState.Instance.hdelta) {
                _maxlevel += HtmState.Instance.hdelta / 2;
                _minlevel -= HtmState.Instance.hdelta / 2;
            }
            //
            // Min must never be below 3
            // 
            int nudgeup = _minlevel - 3;
            if(nudgeup < 0) {
                _minlevel -= nudgeup;
                _maxlevel -= nudgeup;
            }

            _maxlevel += HtmState.Instance.deltalevel;
            _minlevel += HtmState.Instance.deltalevel;

            HtmState.Instance.maxlevel = _maxlevel;
            HtmState.Instance.minlevel = _minlevel;
            HtmState.Instance.appendTrace = false;
            bool appendflag = HtmState.Instance.appendTrace;
            if(_dotrace && _mode == Mode.File) {
                if(_directory == null) {
                    throw (new Exception("Cover.IntersectSetup():Undefined debug directory"));
                }
                this._tracefile = new Trace(String.Concat(_directory, "ix.txt"), appendflag);
                this._savedfile = new Trace(String.Concat(_directory, "i"), appendflag);
                this._extendedfile = new Trace(String.Concat(_directory, "e"), appendflag);
                this._halfspacesfile = new Trace(String.Concat(_directory, "c"), appendflag);
                foreach(Convex con in _reg.ConvexList) {
                    foreach(Halfspace hs in con.HalfspaceList) {
                        _halfspacesfile.dump(hs.Vector.X, hs.Vector.Y, hs.Vector.Z, hs.Cos0, 1);
                    }
                }
            }
            _savePairs.Clear();
            _saveAugPairs.Clear();
            if(!_collectPartials) {
                _maxlevel += 2;
            }
        }
        private void IntersectCleanup(bool rangeMergeFlag) {
            if(_dotrace && _mode == Mode.File) {
                _tracefile.close();
                _savedfile.close();
                _extendedfile.close();
                _halfspacesfile.close();
            }
        }
        private List<Int64AugPair> Xintersect(bool fulleqpartial) {
            IntersectSetup();
            _fulleqpartial = fulleqpartial;
            for(int i = 0; i < 8; i++) {
                // Put them on the queue

                Int64 hid;
                Cartesian v0, v1, v2;
                Cartesian fa, fb, fc;

                // For breadth first, this is replaced
                // witgh creat

                hid = 8 + i;

                Trixel.TopVertices(i, out fa, out fb, out fc);
                v0 = new Cartesian(fa, false);
                v1 = new Cartesian(fb, false);
                v2 = new Cartesian(fc, false);
                theq.Enqueue(new Node(hid, v0, v1, v2, 0));
            }
            traverse();
            IntersectCleanup(false);
            return _saveAugPairs;

        }

        private List<Int64Pair> Xintersect() {

            IntersectSetup();

            for(int i = 0; i < 8; i++) {
                // Put them on the queue

                Int64 hid;
                Cartesian v0, v1, v2;
                Cartesian fa, fb, fc;

                // For breadth first, this is replaced
                // witgh creat

                hid = 8 + i;

                Trixel.TopVertices(i, out fa, out fb, out fc);
                v0 = new Cartesian(fa, false);
                v1 = new Cartesian(fb, false);
                v2 = new Cartesian(fc, false);
                theq.Enqueue(new Node(hid, v0, v1, v2, 0));
            }
            traverse();
            IntersectCleanup(false);
            return _savePairs;
        }
		// -------------------------------------------------------------------------------------- END
		/// <summary>
		/// Check if a hid is in the given list of ranges
		/// 
		/// Used for diagnostic purposes
		/// </summary>
		/// <param name="hid"></param>
		/// <param name="coverlist"></param>
		/// <returns></returns>
		private static bool IsHidInCoverlist(Int64 hid, List<Int64Pair> coverlist) {
			if (!Trixel.IsValid(hid))
				return false;

			if (coverlist == null)
				return false;

			if (coverlist.Count < 1)
				return false;

			Int64Pair hid20 = Trixel.Extend(hid, 20);
			foreach (Int64Pair range in coverlist) {
				//
				// ORDER OF the next two tests is important!
				//
				if (hid20.hi < range.lo)
					return false;
				if (hid20.lo <= range.hi && hid20.hi >= range.lo)
					return true;
			}
			return false;
		}
	}
    /// <summary>
    /// Mode controls whether or not trace information is generated into text files.
    /// 
    /// For debugging purposes only
    /// </summary>
    public enum Mode {
        /// <summary>Normal operation without file i/o</summary>
        Normal,
        /// <summary>Write trace files</summary>
        File,
        /// <summary>Not used</summary>
        NullMode
    };
    /// <summary>
    /// Markup is a classification of a trixel with respect to a convex
    /// </summary>
    public enum Markup {
        /// </summary>trixel is completely inside convex</summary>
        Full,
        /// <summary>trixel is completely outside convex</summary>
        Reject,
        /// <summary>trixel non-trivially intersects convex</summary>
        Partial,
        /// <summary>trixel's status is not known</summary>
        Dontknow
    };
    /// <summary>
    /// Savecode is a local condition that
    /// guides the heuristic algorithm that performs trixel-convex
    /// intersection.
    /// </summary>
    public enum Savecode {
        /// do not save this trixel
        Nosave,
        /// reached maximum depth
        Reachedlevel,
        /// 4 children are partial
        Peq4,
        /// at least 2 children are full			
        Fge2,
        /// 1 child full, AND 3 are partial			
        Peq3Feq1,
        /// at least 1 child is partial,
        /// and parent had three partial children 
        Pgt1Ppreveq3,
        /// <summary>
        /// Full, completely contained
        /// </summary>
        Full,
        Partial,
    };
}
