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
{

	public class HeuristicCover
	{
		/// <summary>
		/// list of hid ranges are collected and manipulated here
		/// </summary>
		private Interval _range;
		/// <summary>
		/// The Region whose cover sill be computed
		/// </summary>
		Region _reg;


		/// <summary>
		/// Variable length (level) indicator
		/// 
		/// If set, trixels are saved at their natural level.
		/// Otherwise hid ranges at level 20 are used.
		/// </summary>
		bool _varlen = false;

		/// <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>
		List<Int64Pair> _savePairs;
		/// <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 Mode _mode = Mode.Normal;
        
        public static List<Int64Pair> HidRange(Region reg, bool varlen, string directory) {

            List<Int64Pair> returnResult;
            HeuristicCover cover = new HeuristicCover(reg, varlen, directory);
            returnResult = cover.Rintersect();
            return returnResult;
        }
		
		public static List<Int64Pair> HidRange(Region reg, bool varlen, string directory, bool useR) {

			List<Int64Pair> returnResult;
			HeuristicCover cover = new HeuristicCover(reg, varlen, directory);
            if(useR) {
                returnResult = cover.Rintersect();
            } else {
                returnResult = cover.Rintersect();
            }
			return returnResult;
		}
		/// <summary>
		/// Return a hid range list
		/// </summary>
		/// <param name="reg"></param>
		/// <returns></returns>
		public static List<Int64Pair> HidRange(Region reg) {
			return HidRange(reg, false, null);

		}
		/// <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>
		public static List<Int64Pair> HidList(Region reg) {
			return HidRange(reg, true, null);
		}

		/// <summary>
		/// Create an instance with optional tracing
		///
		/// </summary>
		/// <param name="directory">where debug trace text files are created</param>]
		public HeuristicCover(Region reg, bool varlen, string dir) {
			_range = new Interval();
			_savePairs = new List<Int64Pair>();
			_dotrace = true;

			_reg = reg;
			_directory = dir;
			_varlen = varlen;


			if (_directory == null) {
				_mode = Mode.Normal;
			} else {
				_mode = Mode.File;
			}			
		}
		/// <summary>
		/// This is the main place of a triangle vs a Convex.
		/// </summary>
		/// <param name="index">one of eight octahedral triangles</param>
		private void testTrixel(int index) {

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

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

			hid = 8 + index;

            Trixel.TopVertices(index, out fa, out fb, out fc);
            v0 = new Cartesian(fa, false);
			v1 = new Cartesian(fb, false);
			v2 = new Cartesian(fc, false);

			depthFirstWLookahead(_minlevel, 0, hid, v0, v1, v2, 0);
		}
 
        
		private void saveTrixel(Int64 hid) {
			Int64Pair lohi;
			int depth = Trixel.LevelOfHid(hid);

			if (_varlen) {
				lohi.lo = hid;
				lohi.hi = -1L;	// < 0 indicates that list is variable length trixels HIDs
			} else {
				lohi = Trixel.Extend(hid, 20);
			}
			if (_dotrace && _mode == Mode.File) {
				_savedfile.dump("{0,20}", hid);
				_extendedfile.dump("{0,20}  {1,20}", lohi.lo, lohi.hi);
			}

            /// if buildist is null, we don't build a list (just exercise the method)

			if (_savePairs != null) {
				_savePairs.Add(lohi);
			} else {
				throw new Exception("EMPTY BUILDLIST");
			}
			return;
		}
		public static void getAxisandAngle(Cartesian P, Cartesian Q, Cartesian axis, out double angle) {
			Cartesian Pc = new Cartesian(P, false);
			Cartesian Qc = new Cartesian(Q, false);

			double n1, n2;
			n1 = Pc.Norm;
			n2 = Qc.Norm;

			angle = Math.Acos(Pc.Dot(Qc));

			Cartesian Tc;
			Tc = Pc.Cross(Qc, true);

			axis.X = Tc.X;
			axis.Y = Tc.Y;
			axis.Z = Tc.Z;
		}
		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;
		}
        
		public 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 == 3085L){ // 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
        };
		/// <summary>
		/// Test a triangle's subtriangle whether they are partial.
		//  if level is nonzero, recurse, unless some conditions are met
		/// </summary>
		/// <param name="hid">a</param>
		/// <param name="v0">a</param>
		/// <param name="v1">a</param>
		/// <param name="v2">a</param>
		/// <param name="PPrev">a</param>
		/// <returns>nothing</returns>
		/// \callgraph
		private void depthFirstWLookahead(int levelsForced, int level, Int64 hid,
			Cartesian v0, Cartesian v1, Cartesian v2, int PPrev) {
			Int64[] ids = new Int64[4];
			Int64 id0;
			Savecode why;
			Markup[] m = new Markup[4];
			int P, F;// count number of partials and fulls

			Cartesian w0 = v1.Add(v2);			w0.Normalize();
			Cartesian w1 = v2.Add(v0);			w1.Normalize();
			Cartesian w2 = v0.Add(v1);			w2.Normalize();

			ids[0] = id0 = hid << 2;
			ids[1] = id0 + 1;
			ids[2] = id0 + 2;
			ids[3] = id0 + 3;

			m[0] = testNode(v0, w2, w1, ids[0]);
			m[1] = testNode(v1, w0, w2, ids[1]);
			m[2] = testNode(v2, w1, w0, ids[2]);
			m[3] = testNode(w0, w1, w2, ids[3]);

			F = 0;
			if (m[0] == Markup.Full) F++;
			if (m[1] == Markup.Full) F++;
			if (m[2] == Markup.Full) F++;
			if (m[3] == Markup.Full) F++;

			P = 0;
			if (m[0] == Markup.Partial) P++;
			if (m[1] == Markup.Partial) P++;
			if (m[2] == Markup.Partial) P++;
			if (m[3] == Markup.Partial) P++;

			//
			// Several interesting cases for saving this (the parent) trixel.
			// Case P==4, all four children are partials, so pretend parent is full, we save and return
			// Case P==3, and F==1, most of the parent is in, so pretend that parent is full again
			// Case P==2 or 3, but the previous testPartial had three partials, so parent was in an arc
			// as opposed to previous partials being fewer, so parent was in a tiny corner...
			why = Savecode.Nosave;
			//
			// The Big Idea: several heuristics cut off further subdivisions
			// but this may cause huge overshoots!
			// So, I added an option to tweak the heuristics
			// with a flagarray to ignore certain heuristic on demand
			//
			// code  abbrev         summary
			// ----------------------------------------------------
			// p4    peq4	        all four children are partial
			// f2    fge2           more than two children are full
			// p3    peq3feq1       3 are partial one is full
			// p1    pgt1ppreveq3   more than 1 partial, parent had three partials
			// a good halfspace to test this on is
			// (1, 0.8, 0.5, 0.9998) normalized, of course.
			// 

			//  The rule P>1&Pprev=3 may get triggered to early causing huge
			//  overshoots. So, we force it 
			//  maxlevel - level to go
			

			if (level >= _maxlevel) why = Savecode.Reachedlevel;
			else if (P == 4) why = Savecode.Peq4;
			else if (F >= 2) why = Savecode.Fge2;
			else if ((P == 3) && (F == 1)) why = Savecode.Peq3Feq1;
			else if ((P > 1 && PPrev == 3)) why = Savecode.Pgt1Ppreveq3;
			

			if (levelsForced > 0) 
				why = Savecode.Nosave; 

			if (why != Savecode.Nosave) {
				if (_dotrace && _mode == Mode.File) {
					this._tracefile.dump("Saving {0} P why {1} (level {2}) ",
						hid, why, level);
					
				}
				saveTrixel(hid);
				return;
			} else {
				// look at each child, see if some need saving;
				if (F == 4 && (levelsForced <= 0))
					throw new ApplicationException("depthFirst inconsistency 1");
				for (int i = 0; i < 4; i++) {
					if (m[i] == Markup.Full) {
						//
						// <LOG>"saving <hid> code Child is full"</LOG>
						if (_dotrace && _mode == Mode.File) {
							this._tracefile.dump("Saving {0} F why {1} (level {2})",
								ids[i], Savecode.Full, level);
						}
						saveTrixel(ids[i]);
					}
				}
				// look at the four kids again, for partials
				level++; // <--- moved here from above, see legacy comment
				levelsForced--;
				if (m[0] == Markup.Partial) depthFirstWLookahead(levelsForced, level, ids[0], v0, w2, w1, P);
				if (m[1] == Markup.Partial) depthFirstWLookahead(levelsForced, level, ids[1], v1, w0, w2, P);
				if (m[2] == Markup.Partial) depthFirstWLookahead(levelsForced, level, ids[2], v2, w1, w0, P);
				if (m[3] == Markup.Partial) depthFirstWLookahead(levelsForced, level, ids[3], w0, w1, w2, P);
			}
		}
		private void checkplist(List<Int64Pair> plist) {
			for (int i = 1; i < plist.Count; i++) {
				if (plist[i - 1].hi > plist[i].lo) {
					StringBuilder sb = new StringBuilder();
					sb.AppendFormat("checkplist at line {0}: {1} - {2}", i, plist[i - 1].lo,
						plist[i - 1].hi);
					sb.AppendFormat("  after {0} - {1}", plist[i].lo, plist[i].hi);
					throw new ApplicationException(sb.ToString());
				}
			}
		}
        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();
        }
        private void IntersectCleanup(bool rangeMergeFlag) {
            if(_dotrace && _mode == Mode.File) {
                _tracefile.close();
                _savedfile.close();
                _extendedfile.close();
                _halfspacesfile.close();
            }
            //if(!_varlen) {
            //    _savePairs.Sort(Interval.ComparisonRanges);
            //    Interval.compact(_savePairs);
            //    Interval.IsWellFormed(_savePairs);
            //    if(rangeMergeFlag) {
            //        _range.merge(_savePairs);
            //    }
            //}
        }
       
		private List<Int64Pair> Rintersect() {

			_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("ConvextTrixel: Undefined debug directory (x4intersct)"));
				}
				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();

			for (int i = 0; i < 8; i++) {
				testTrixel(i);
			}

			if (_dotrace && _mode == Mode.File) {
				_tracefile.close();
				_savedfile.close();
				_extendedfile.close();
				_halfspacesfile.close();
			}

			if (!_varlen) {
				_savePairs.Sort(Interval.ComparisonRanges);
				Interval.compact(_savePairs);
				Interval.IsWellFormed(_savePairs);
				_range.merge(_savePairs);

				//_range.check();
				//_range.compact();
				//_range.check();

				return _range.pairList;
			} else {
				return _savePairs;
			}
		}

		// -------------------------------------------------------------------------------------- END
		/// <summary>
		/// Check if a hid is in the give list of ranges
		/// 
		/// Used for diagnostic purposes
		/// </summary>
		/// <param name="hid"></param>
		/// <param name="coverlist"></param>
		/// <returns></returns>
		public 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;
		}
	}
}
