// The following handle manipulation primitives are all described by Guibas
        // and Stolfi. However, Guibas and Stolfi use an edge-based data structure,
        // whereas I use a triangle-based data structure.

        /// <summary>
        /// Find the abutting triangle; same edge. [sym(abc) -> ba*]
        /// </summary>
        /// <remarks>
        /// Note that the edge direction is necessarily reversed, because the handle specified 
        /// by an oriented triangle is directed counterclockwise around the triangle.
        /// </remarks>
        public void Sym(ref Otri o2)
        {
            //o2 = tri.triangles[orient];
            // decode(ptr, otri2);

            o2.triangle = triangle.neighbors[orient].triangle;
            o2.orient = triangle.neighbors[orient].orient;
        }

        /// <summary>
        /// Find a new location for a Steiner point.
        /// </summary>
        /// <param name="torg"></param>
        /// <param name="tdest"></param>
        /// <param name="tapex"></param>
        /// <param name="xi"></param>
        /// <param name="eta"></param>
        /// <param name="offcenter"></param>
        /// <param name="badotri"></param>
        /// <returns></returns>
        public Point FindLocation(Vertex torg, Vertex tdest, Vertex tapex,
            ref double xi, ref double eta, bool offcenter, Otri badotri)
        {
            // Based on using -U switch, call the corresponding function
            if (behavior.MaxAngle == 0.0)
            {
                return FindNewLocationWithoutMaxAngle(torg, tdest, tapex, ref xi, ref eta, true, badotri);
            }

            // With max angle
            return FindNewLocation(torg, tdest, tapex, ref xi, ref eta, true, badotri);
        }

        /// <summary>
        /// Finds the star of a given point.
        /// </summary>
        /// <param name="badotri"></param>
        /// <param name="p"></param>
        /// <param name="q"></param>
        /// <param name="r"></param>
        /// <param name="whichPoint"></param>
        /// <param name="points">List of points on the star of the given point.</param>
        /// <returns>Number of points on the star of the given point.</returns>
        private int GetStarPoints(Otri badotri, Vertex p, Vertex q, Vertex r,
                    int whichPoint, ref double[] points)
        {

            Otri neighotri = default(Otri);  // for return value of the function
            Otri tempotri;   // for temporary usage
            double first_x = 0, first_y = 0;	  // keeps the first point to be considered
            double second_x = 0, second_y = 0;  // for determining the edge we will begin
            double third_x = 0, third_y = 0;	  // termination
            double[] returnPoint = new double[2];	  // for keeping the returned point	
            int numvertices = 0;	  // for keeping number of surrounding vertices

            // first determine which point to be used to find its neighbor triangles
            switch (whichPoint)
            {
                case 1:
                    first_x = p.x;	// point at the center
                    first_y = p.y;
                    second_x = r.x; // second vertex of first edge to consider
                    second_y = r.y;
                    third_x = q.x;  // for terminating the search
                    third_y = q.y;
                    break;
                case 2:
                    first_x = q.x;  // point at the center
                    first_y = q.y;
                    second_x = p.x; // second vertex of first edge to consider
                    second_y = p.y;
                    third_x = r.x;	// for terminating the search
                    third_y = r.y;
                    break;
                case 3:
                    first_x = r.x;	// point at the center
                    first_y = r.y;
                    second_x = q.x; // second vertex of first edge to consider
                    second_y = q.y;
                    third_x = p.x;	// for terminating the search
                    third_y = p.y;
                    break;
            }
            tempotri = badotri;
            // add first point as the end of first edge
            points[numvertices] = second_x;
            numvertices++;
            points[numvertices] = second_y;
            numvertices++;
            // assign as dummy value
            returnPoint[0] = second_x; returnPoint[1] = second_y;
            // until we reach the third point of the beginning triangle	
            do
            {
                // find the neighbor's third point where it is incident to given edge
                if (!GetNeighborsVertex(tempotri, first_x, first_y, second_x, second_y, ref returnPoint, ref neighotri))
                {
                    // go to next triangle
                    tempotri = neighotri;
                    // now the second point is the neighbor's third vertex			
                    second_x = returnPoint[0];
                    second_y = returnPoint[1];
                    // add a new point to the list of surrounding points
                    points[numvertices] = returnPoint[0];
                    numvertices++;
                    points[numvertices] = returnPoint[1];
                    numvertices++;
                }
                else
                {
                    numvertices = 0;
                    break;
                }

            } while (!((Math.Abs(returnPoint[0] - third_x) <= EPS) &&
                     (Math.Abs(returnPoint[1] - third_y) <= EPS)));
            return numvertices / 2;

        }

        /// <summary>
        /// Removes ghost triangles.
        /// </summary>
        /// <param name="startghost"></param>
        /// <returns>Number of vertices on the hull.</returns>
        int RemoveGhosts(ref Otri startghost)
        {
            Otri searchedge = default(Otri);
            Otri dissolveedge = default(Otri);
            Otri deadtriangle = default(Otri);
            Vertex markorg;
            int hullsize;

            bool noPoly = !mesh.behavior.Poly;

            // Find an edge on the convex hull to start point location from.
            startghost.Lprev(ref searchedge);
            searchedge.SymSelf();
            Mesh.dummytri.neighbors[0] = searchedge;
            // Remove the bounding box and count the convex hull edges.
            startghost.Copy(ref dissolveedge);
            hullsize = 0;
            do
            {
                hullsize++;
                dissolveedge.Lnext(ref deadtriangle);
                dissolveedge.LprevSelf();
                dissolveedge.SymSelf();

                // If no PSLG is involved, set the boundary markers of all the vertices
                // on the convex hull.  If a PSLG is used, this step is done later.
                if (noPoly)
                {
                    // Watch out for the case where all the input vertices are collinear.
                    if (dissolveedge.triangle != Mesh.dummytri)
                    {
                        markorg = dissolveedge.Org();
                        if (markorg.mark == 0)
                        {
                            markorg.mark = 1;
                        }
                    }
                }
                // Remove a bounding triangle from a convex hull triangle.
                dissolveedge.Dissolve();
                // Find the next bounding triangle.
                deadtriangle.Sym(ref dissolveedge);

                // Delete the bounding triangle.
                mesh.TriangleDealloc(deadtriangle.triangle);
            } while (!dissolveedge.Equal(startghost));

            return hullsize;
        }

        void Check4DeadEvent(ref Otri checktri, SweepEvent[] eventheap, ref int heapsize)
        {
            SweepEvent deadevent;
            SweepEventVertex eventvertex;
            int eventnum = -1;

            eventvertex = checktri.Org() as SweepEventVertex;
            if (eventvertex != null)
            {
                deadevent = eventvertex.evt;
                eventnum = deadevent.heapposition;

                HeapDelete(eventheap, heapsize, eventnum);
                heapsize--;
                checktri.SetOrg(null);
            }
        }

        SplayNode CircleTopInsert(SplayNode splayroot, Otri newkey,
                                  Vertex pa, Vertex pb, Vertex pc, double topy)
        {
            double ccwabc;
            double xac, yac, xbc, ybc;
            double aclen2, bclen2;
            Point searchpoint = new Point(); // TODO: mesh.nextras
            Otri dummytri = default(Otri);

            ccwabc = Primitives.CounterClockwise(pa, pb, pc);
            xac = pa.x - pc.x;
            yac = pa.y - pc.y;
            xbc = pb.x - pc.x;
            ybc = pb.y - pc.y;
            aclen2 = xac * xac + yac * yac;
            bclen2 = xbc * xbc + ybc * ybc;
            searchpoint.x = pc.x - (yac * bclen2 - ybc * aclen2) / (2.0 * ccwabc);
            searchpoint.y = topy;
            return SplayInsert(Splay(splayroot, searchpoint, ref dummytri), newkey, searchpoint);
        }

        SplayNode Splay(SplayNode splaytree, Point searchpoint, ref Otri searchtri)
        {
            SplayNode child, grandchild;
            SplayNode lefttree, righttree;
            SplayNode leftright;
            Vertex checkvertex;
            bool rightofroot, rightofchild;

            if (splaytree == null)
            {
                return null;
            }
            checkvertex = splaytree.keyedge.Dest();
            if (checkvertex == splaytree.keydest)
            {
                rightofroot = RightOfHyperbola(ref splaytree.keyedge, searchpoint);
                if (rightofroot)
                {
                    splaytree.keyedge.Copy(ref searchtri);
                    child = splaytree.rchild;
                }
                else
                {
                    child = splaytree.lchild;
                }
                if (child == null)
                {
                    return splaytree;
                }
                checkvertex = child.keyedge.Dest();
                if (checkvertex != child.keydest)
                {
                    child = Splay(child, searchpoint, ref searchtri);
                    if (child == null)
                    {
                        if (rightofroot)
                        {
                            splaytree.rchild = null;
                        }
                        else
                        {
                            splaytree.lchild = null;
                        }
                        return splaytree;
                    }
                }
                rightofchild = RightOfHyperbola(ref child.keyedge, searchpoint);
                if (rightofchild)
                {
                    child.keyedge.Copy(ref searchtri);
                    grandchild = Splay(child.rchild, searchpoint, ref searchtri);
                    child.rchild = grandchild;
                }
                else
                {
                    grandchild = Splay(child.lchild, searchpoint, ref searchtri);
                    child.lchild = grandchild;
                }
                if (grandchild == null)
                {
                    if (rightofroot)
                    {
                        splaytree.rchild = child.lchild;
                        child.lchild = splaytree;
                    }
                    else
                    {
                        splaytree.lchild = child.rchild;
                        child.rchild = splaytree;
                    }
                    return child;
                }
                if (rightofchild)
                {
                    if (rightofroot)
                    {
                        splaytree.rchild = child.lchild;
                        child.lchild = splaytree;
                    }
                    else
                    {
                        splaytree.lchild = grandchild.rchild;
                        grandchild.rchild = splaytree;
                    }
                    child.rchild = grandchild.lchild;
                    grandchild.lchild = child;
                }
                else
                {
                    if (rightofroot)
                    {
                        splaytree.rchild = grandchild.lchild;
                        grandchild.lchild = splaytree;
                    }
                    else
                    {
                        splaytree.lchild = child.rchild;
                        child.rchild = splaytree;
                    }
                    child.lchild = grandchild.rchild;
//.........这里部分代码省略.........

        /// <summary>
        /// Insert a vertex into a Delaunay triangulation, performing flips as necessary 
        /// to maintain the Delaunay property.
        /// </summary>
        /// <param name="newvertex">The point to be inserted.</param>
        /// <param name="searchtri">The triangle to start the search.</param>
        /// <param name="splitseg">Segment to split.</param>
        /// <param name="segmentflaws">Check for creation of encroached subsegments.</param>
        /// <param name="triflaws">Check for creation of bad quality triangles.</param>
        /// <returns>If a duplicate vertex or violated segment does not prevent the 
        /// vertex from being inserted, the return value will be ENCROACHINGVERTEX if 
        /// the vertex encroaches upon a subsegment (and checking is enabled), or
        /// SUCCESSFULVERTEX otherwise. In either case, 'searchtri' is set to a handle
        /// whose origin is the newly inserted vertex.</returns>
        /// <remarks>
        /// The point 'newvertex' is located. If 'searchtri.triangle' is not NULL,
        /// the search for the containing triangle begins from 'searchtri'.  If
        /// 'searchtri.triangle' is NULL, a full point location procedure is called.
        /// If 'insertvertex' is found inside a triangle, the triangle is split into
        /// three; if 'insertvertex' lies on an edge, the edge is split in two,
        /// thereby splitting the two adjacent triangles into four. Edge flips are
        /// used to restore the Delaunay property. If 'insertvertex' lies on an
        /// existing vertex, no action is taken, and the value DUPLICATEVERTEX is
        /// returned. On return, 'searchtri' is set to a handle whose origin is the
        /// existing vertex.
        ///
        /// InsertVertex() does not use flip() for reasons of speed; some
        /// information can be reused from edge flip to edge flip, like the
        /// locations of subsegments.
        /// 
        /// Param 'splitseg': Normally, the parameter 'splitseg' is set to NULL, 
        /// implying that no subsegment should be split. In this case, if 'insertvertex' 
        /// is found to lie on a segment, no action is taken, and the value VIOLATINGVERTEX 
        /// is returned. On return, 'searchtri' is set to a handle whose primary edge is the 
        /// violated subsegment.
        /// If the calling routine wishes to split a subsegment by inserting a vertex in it, 
        /// the parameter 'splitseg' should be that subsegment. In this case, 'searchtri' 
        /// MUST be the triangle handle reached by pivoting from that subsegment; no point 
        /// location is done.
        /// 
        /// Param 'segmentflaws': Flags that indicate whether or not there should
        /// be checks for the creation of encroached subsegments. If a newly inserted 
        /// vertex encroaches upon subsegments, these subsegments are added to the list 
        /// of subsegments to be split if 'segmentflaws' is set.
        /// 
        /// Param 'triflaws': Flags that indicate whether or not there should be
        /// checks for the creation of bad quality triangles. If bad triangles are 
        /// created, these are added to the queue if 'triflaws' is set.
        /// </remarks>
        internal InsertVertexResult InsertVertex(Vertex newvertex, ref Otri searchtri,
            ref Osub splitseg, bool segmentflaws, bool triflaws)
        {
            Otri horiz = default(Otri);
            Otri top = default(Otri);
            Otri botleft = default(Otri), botright = default(Otri);
            Otri topleft = default(Otri), topright = default(Otri);
            Otri newbotleft = default(Otri), newbotright = default(Otri);
            Otri newtopright = default(Otri);
            Otri botlcasing = default(Otri), botrcasing = default(Otri);
            Otri toplcasing = default(Otri), toprcasing = default(Otri);
            Otri testtri = default(Otri);
            Osub botlsubseg = default(Osub), botrsubseg = default(Osub);
            Osub toplsubseg = default(Osub), toprsubseg = default(Osub);
            Osub brokensubseg = default(Osub);
            Osub checksubseg = default(Osub);
            Osub rightsubseg = default(Osub);
            Osub newsubseg = default(Osub);
            BadSubseg encroached;
            //FlipStacker newflip;
            Vertex first;
            Vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
            Vertex segmentorg, segmentdest;
            int region;
            double area;
            InsertVertexResult success;
            LocateResult intersect;
            bool doflip;
            bool mirrorflag;
            bool enq;

            if (splitseg.seg == null)
            {
                // Find the location of the vertex to be inserted.  Check if a good
                // starting triangle has already been provided by the caller.
                if (searchtri.triangle == dummytri)
                {
                    // Find a boundary triangle.
                    horiz.triangle = dummytri;
                    horiz.orient = 0;
                    horiz.SymSelf();
                    // Search for a triangle containing 'newvertex'.
                    intersect = locator.Locate(newvertex, ref horiz);
                }
                else
                {
                    // Start searching from the triangle provided by the caller.
                    searchtri.Copy(ref horiz);
                    intersect = locator.PreciseLocate(newvertex, ref horiz, true);
                }
            }
//.........这里部分代码省略.........

        /// <summary>
        /// Transform two triangles to two different triangles by flipping an edge 
        /// counterclockwise within a quadrilateral.
        /// </summary>
        /// <param name="flipedge">Handle to the edge that will be flipped.</param>
        /// <remarks>Imagine the original triangles, abc and bad, oriented so that the
        /// shared edge ab lies in a horizontal plane, with the vertex b on the left
        /// and the vertex a on the right. The vertex c lies below the edge, and
        /// the vertex d lies above the edge. The 'flipedge' handle holds the edge
        /// ab of triangle abc, and is directed left, from vertex a to vertex b.
        ///
        /// The triangles abc and bad are deleted and replaced by the triangles cdb
        /// and dca.  The triangles that represent abc and bad are NOT deallocated;
        /// they are reused for dca and cdb, respectively.  Hence, any handles that
        /// may have held the original triangles are still valid, although not
        /// directed as they were before.
        ///
        /// Upon completion of this routine, the 'flipedge' handle holds the edge
        /// dc of triangle dca, and is directed down, from vertex d to vertex c.
        /// (Hence, the two triangles have rotated counterclockwise.)
        ///
        /// WARNING:  This transformation is geometrically valid only if the
        /// quadrilateral adbc is convex.  Furthermore, this transformation is
        /// valid only if there is not a subsegment between the triangles abc and
        /// bad.  This routine does not check either of these preconditions, and
        /// it is the responsibility of the calling routine to ensure that they are
        /// met.  If they are not, the streets shall be filled with wailing and
        /// gnashing of teeth.
        /// 
        /// Terminology
        ///
        /// A "local transformation" replaces a small set of triangles with another
        /// set of triangles.  This may or may not involve inserting or deleting a
        /// vertex.
        ///
        /// The term "casing" is used to describe the set of triangles that are
        /// attached to the triangles being transformed, but are not transformed
        /// themselves.  Think of the casing as a fixed hollow structure inside
        /// which all the action happens.  A "casing" is only defined relative to
        /// a single transformation; each occurrence of a transformation will
        /// involve a different casing.
        /// </remarks>
        internal void Flip(ref Otri flipedge)
        {
            Otri botleft = default(Otri), botright = default(Otri);
            Otri topleft = default(Otri), topright = default(Otri);
            Otri top = default(Otri);
            Otri botlcasing = default(Otri), botrcasing = default(Otri);
            Otri toplcasing = default(Otri), toprcasing = default(Otri);
            Osub botlsubseg = default(Osub), botrsubseg = default(Osub);
            Osub toplsubseg = default(Osub), toprsubseg = default(Osub);
            Vertex leftvertex, rightvertex, botvertex;
            Vertex farvertex;

            // Identify the vertices of the quadrilateral.
            rightvertex = flipedge.Org();
            leftvertex = flipedge.Dest();
            botvertex = flipedge.Apex();
            flipedge.Sym(ref top);

            // SELF CHECK

            //if (top.triangle == dummytri)
            //{
            //    logger.Error("Attempt to flip on boundary.", "Mesh.Flip()");
            //    flipedge.LnextSelf();
            //    return;
            //}

            //if (checksegments)
            //{
            //    flipedge.SegPivot(ref toplsubseg);
            //    if (toplsubseg.ss != dummysub)
            //    {
            //        logger.Error("Attempt to flip a segment.", "Mesh.Flip()");
            //        flipedge.LnextSelf();
            //        return;
            //    }
            //}

            farvertex = top.Apex();

            // Identify the casing of the quadrilateral.
            top.Lprev(ref topleft);
            topleft.Sym(ref toplcasing);
            top.Lnext(ref topright);
            topright.Sym(ref toprcasing);
            flipedge.Lnext(ref botleft);
            botleft.Sym(ref botlcasing);
            flipedge.Lprev(ref botright);
            botright.Sym(ref botrcasing);
            // Rotate the quadrilateral one-quarter turn counterclockwise.
            topleft.Bond(ref botlcasing);
            botleft.Bond(ref botrcasing);
            botright.Bond(ref toprcasing);
            topright.Bond(ref toplcasing);

            if (checksegments)
            {
                // Check for subsegments and rebond them to the quadrilateral.
//.........这里部分代码省略.........

 public void Update(ref Otri otri)
 {
     otri.Copy(ref recenttri);
 }

        /// <summary>
        /// Find a triangle or edge containing a given point.
        /// </summary>
        /// <param name="searchpoint">The point to locate.</param>
        /// <param name="searchtri">The triangle to start the search at.</param>
        /// <returns>Location information.</returns>
        /// <remarks>
        /// Searching begins from one of:  the input 'searchtri', a recently
        /// encountered triangle 'recenttri', or from a triangle chosen from a
        /// random sample. The choice is made by determining which triangle's
        /// origin is closest to the point we are searching for. Normally,
        /// 'searchtri' should be a handle on the convex hull of the triangulation.
        ///
        /// Details on the random sampling method can be found in the Mucke, Saias,
        /// and Zhu paper cited in the header of this code.
        ///
        /// On completion, 'searchtri' is a triangle that contains 'searchpoint'.
        ///
        /// Returns ONVERTEX if the point lies on an existing vertex. 'searchtri'
        /// is a handle whose origin is the existing vertex.
        ///
        /// Returns ONEDGE if the point lies on a mesh edge. 'searchtri' is a
        /// handle whose primary edge is the edge on which the point lies.
        ///
        /// Returns INTRIANGLE if the point lies strictly within a triangle.
        /// 'searchtri' is a handle on the triangle that contains the point.
        ///
        /// Returns OUTSIDE if the point lies outside the mesh. 'searchtri' is a
        /// handle whose primary edge the point is to the right of.  This might
        /// occur when the circumcenter of a triangle falls just slightly outside
        /// the mesh due to floating-point roundoff error. It also occurs when
        /// seeking a hole or region point that a foolish user has placed outside
        /// the mesh.
        ///
        /// WARNING:  This routine is designed for convex triangulations, and will
        /// not generally work after the holes and concavities have been carved.
        /// </remarks>
        public LocateResult Locate(Point searchpoint, ref Otri searchtri)
        {
            Otri sampletri = default(Otri);
            Vertex torg, tdest;
            float searchdist, dist;
            float ahead;

            // Record the distance from the suggested starting triangle to the
            // point we seek.
            torg = searchtri.Org();
            searchdist = (searchpoint.X - torg.x) * (searchpoint.X - torg.x) +
                         (searchpoint.Y - torg.y) * (searchpoint.Y - torg.y);

            // If a recently encountered triangle has been recorded and has not been
            // deallocated, test it as a good starting point.
            if (recenttri.triangle != null)
            {
                if (!Otri.IsDead(recenttri.triangle))
                {
                    torg = recenttri.Org();
                    if ((torg.x == searchpoint.X) && (torg.y == searchpoint.Y))
                    {
                        recenttri.Copy(ref searchtri);
                        return LocateResult.OnVertex;
                    }
                    dist = (searchpoint.X - torg.x) * (searchpoint.X - torg.x) +
                           (searchpoint.Y - torg.y) * (searchpoint.Y - torg.y);
                    if (dist < searchdist)
                    {
                        recenttri.Copy(ref searchtri);
                        searchdist = dist;
                    }
                }
            }

            // TODO: Improve sampling.
            sampler.Update(mesh);
            int[] samples = sampler.GetSamples(mesh);

            foreach (var key in samples)
            {
                sampletri.triangle = mesh.triangles[key];
                if (!Otri.IsDead(sampletri.triangle))
                {
                    torg = sampletri.Org();
                    dist = (searchpoint.X - torg.x) * (searchpoint.X - torg.x) +
                           (searchpoint.Y - torg.y) * (searchpoint.Y - torg.y);
                    if (dist < searchdist)
                    {
                        sampletri.Copy(ref searchtri);
                        searchdist = dist;
                    }
                }
            }

            // Where are we?
            torg = searchtri.Org();
            tdest = searchtri.Dest();
            // Check the starting triangle's vertices.
            if ((torg.x == searchpoint.X) && (torg.y == searchpoint.Y))
            {
                return LocateResult.OnVertex;
            }
//.........这里部分代码省略.........

        /// <summary>
        /// Find a triangle or edge containing a given point.
        /// </summary>
        /// <param name="searchpoint">The point to locate.</param>
        /// <param name="searchtri">The triangle to start the search at.</param>
        /// <param name="stopatsubsegment"> If 'stopatsubsegment' is set, the search 
        /// will stop if it tries to walk through a subsegment, and will return OUTSIDE.</param>
        /// <returns>Location information.</returns>
        /// <remarks>
        /// Begins its search from 'searchtri'. It is important that 'searchtri'
        /// be a handle with the property that 'searchpoint' is strictly to the left
        /// of the edge denoted by 'searchtri', or is collinear with that edge and
        /// does not intersect that edge. (In particular, 'searchpoint' should not
        /// be the origin or destination of that edge.)
        ///
        /// These conditions are imposed because preciselocate() is normally used in
        /// one of two situations:
        ///
        /// (1)  To try to find the location to insert a new point.  Normally, we
        ///      know an edge that the point is strictly to the left of. In the
        ///      incremental Delaunay algorithm, that edge is a bounding box edge.
        ///      In Ruppert's Delaunay refinement algorithm for quality meshing,
        ///      that edge is the shortest edge of the triangle whose circumcenter
        ///      is being inserted.
        ///
        /// (2)  To try to find an existing point.  In this case, any edge on the
        ///      convex hull is a good starting edge. You must screen out the
        ///      possibility that the vertex sought is an endpoint of the starting
        ///      edge before you call preciselocate().
        ///
        /// On completion, 'searchtri' is a triangle that contains 'searchpoint'.
        ///
        /// This implementation differs from that given by Guibas and Stolfi.  It
        /// walks from triangle to triangle, crossing an edge only if 'searchpoint'
        /// is on the other side of the line containing that edge. After entering
        /// a triangle, there are two edges by which one can leave that triangle.
        /// If both edges are valid ('searchpoint' is on the other side of both
        /// edges), one of the two is chosen by drawing a line perpendicular to
        /// the entry edge (whose endpoints are 'forg' and 'fdest') passing through
        /// 'fapex'. Depending on which side of this perpendicular 'searchpoint'
        /// falls on, an exit edge is chosen.
        ///
        /// This implementation is empirically faster than the Guibas and Stolfi
        /// point location routine (which I originally used), which tends to spiral
        /// in toward its target.
        ///
        /// Returns ONVERTEX if the point lies on an existing vertex. 'searchtri'
        /// is a handle whose origin is the existing vertex.
        ///
        /// Returns ONEDGE if the point lies on a mesh edge. 'searchtri' is a
        /// handle whose primary edge is the edge on which the point lies.
        ///
        /// Returns INTRIANGLE if the point lies strictly within a triangle.
        /// 'searchtri' is a handle on the triangle that contains the point.
        ///
        /// Returns OUTSIDE if the point lies outside the mesh. 'searchtri' is a
        /// handle whose primary edge the point is to the right of.  This might
        /// occur when the circumcenter of a triangle falls just slightly outside
        /// the mesh due to floating-point roundoff error. It also occurs when
        /// seeking a hole or region point that a foolish user has placed outside
        /// the mesh.
        ///
        /// WARNING:  This routine is designed for convex triangulations, and will
        /// not generally work after the holes and concavities have been carved.
        /// However, it can still be used to find the circumcenter of a triangle, as
        /// long as the search is begun from the triangle in question.</remarks>
        public LocateResult PreciseLocate(Point searchpoint, ref Otri searchtri,
                                        bool stopatsubsegment)
        {
            Otri backtracktri = default(Otri);
            Osub checkedge = default(Osub);
            Vertex forg, fdest, fapex;
            float orgorient, destorient;
            bool moveleft;

            // Where are we?
            forg = searchtri.Org();
            fdest = searchtri.Dest();
            fapex = searchtri.Apex();
            while (true)
            {
                // Check whether the apex is the point we seek.
                if ((fapex.x == searchpoint.X) && (fapex.y == searchpoint.Y))
                {
                    searchtri.LprevSelf();
                    return LocateResult.OnVertex;
                }
                // Does the point lie on the other side of the line defined by the
                // triangle edge opposite the triangle's destination?
                destorient = Primitives.CounterClockwise(forg, fapex, searchpoint);
                // Does the point lie on the other side of the line defined by the
                // triangle edge opposite the triangle's origin?
                orgorient = Primitives.CounterClockwise(fapex, fdest, searchpoint);
                if (destorient > 0.0)
                {
                    if (orgorient > 0.0)
                    {
                        // Move left if the inner product of (fapex - searchpoint) and
                        // (fdest - forg) is positive.  This is equivalent to drawing
                        // a line perpendicular to the line (forg, fdest) and passing
//.........这里部分代码省略.........

        /// <summary>
        /// Find a new location for a Steiner point.
        /// </summary>
        /// <param name="torg"></param>
        /// <param name="tdest"></param>
        /// <param name="tapex"></param>
        /// <param name="circumcenter"></param>
        /// <param name="xi"></param>
        /// <param name="eta"></param>
        /// <param name="offcenter"></param>
        /// <param name="badotri"></param>
        private Point FindNewLocation(Vertex torg, Vertex tdest, Vertex tapex,
            ref double xi, ref double eta, bool offcenter, Otri badotri)
        {
            double offconstant = behavior.offconstant;

            // for calculating the distances of the edges
            double xdo, ydo, xao, yao, xda, yda;
            double dodist, aodist, dadist;
            // for exact calculation
            double denominator;
            double dx, dy, dxoff, dyoff;

            ////////////////////////////// HALE'S VARIABLES //////////////////////////////
            // keeps the difference of coordinates edge 
            double xShortestEdge = 0, yShortestEdge = 0 /*, xMiddleEdge, yMiddleEdge, xLongestEdge, yLongestEdge*/;

            // keeps the square of edge lengths
            double shortestEdgeDist = 0, middleEdgeDist = 0, longestEdgeDist = 0;

            // keeps the vertices according to the angle incident to that vertex in a triangle
            Point smallestAngleCorner, middleAngleCorner, largestAngleCorner;

            // keeps the type of orientation if the triangle
            int orientation = 0;
            // keeps the coordinates of circumcenter of itself and neighbor triangle circumcenter	
            Point myCircumcenter, neighborCircumcenter;

            // keeps if bad triangle is almost good or not
            int almostGood = 0;
            // keeps the cosine of the largest angle
            double cosMaxAngle;
            bool isObtuse; // 1: obtuse 0: nonobtuse
            // keeps the radius of petal
            double petalRadius;
            // for calculating petal center
            double xPetalCtr_1, yPetalCtr_1, xPetalCtr_2, yPetalCtr_2, xPetalCtr, yPetalCtr, xMidOfShortestEdge, yMidOfShortestEdge;
            double dxcenter1, dycenter1, dxcenter2, dycenter2;
            // for finding neighbor
            Otri neighborotri = default(Otri);
            double[] thirdPoint = new double[2];
            //int neighborNotFound = -1;
            // for keeping the vertices of the neighbor triangle
            Vertex neighborvertex_1;
            Vertex neighborvertex_2;
            Vertex neighborvertex_3;
            // dummy variables 
            double xi_tmp = 0, eta_tmp = 0;
            //vertex thirdVertex;
            // for petal intersection
            double vector_x, vector_y, xMidOfLongestEdge, yMidOfLongestEdge, inter_x, inter_y;
            double[] p = new double[5], voronoiOrInter = new double[4];
            bool isCorrect;

            // for vector calculations in perturbation
            double ax, ay, d;
            double pertConst = 0.06; // perturbation constant

            double lengthConst = 1; // used at comparing circumcenter's distance to proposed point's distance
            double justAcute = 1; // used for making the program working for one direction only
            // for smoothing
            int relocated = 0;// used to differentiate between calling the deletevertex and just proposing a steiner point
            double[] newloc = new double[2];   // new location suggested by smoothing
            double origin_x = 0, origin_y = 0; // for keeping torg safe
            Otri delotri; // keeping the original orientation for relocation process
            // keeps the first and second direction suggested points
            double dxFirstSuggestion, dyFirstSuggestion, dxSecondSuggestion, dySecondSuggestion;
            // second direction variables
            double xMidOfMiddleEdge, yMidOfMiddleEdge;

            double minangle;	// in order to make sure that the circumcircle of the bad triangle is greater than petal
            // for calculating the slab
            double linepnt1_x, linepnt1_y, linepnt2_x, linepnt2_y;	// two points of the line
            double line_inter_x = 0, line_inter_y = 0;
            double line_vector_x, line_vector_y;
            double[] line_p = new double[3]; // used for getting the return values of functions related to line intersection
            double[] line_result = new double[4];
            // intersection of slab and the petal
            double petal_slab_inter_x_first, petal_slab_inter_y_first, petal_slab_inter_x_second, petal_slab_inter_y_second, x_1, y_1, x_2, y_2;
            double petal_bisector_x, petal_bisector_y, dist;
            double alpha;
            bool neighborNotFound_first;
            bool neighborNotFound_second;
            ////////////////////////////// END OF HALE'S VARIABLES //////////////////////////////

            Statistic.CircumcenterCount++;

            // Compute the circumcenter of the triangle.
            xdo = tdest.x - torg.x;
            ydo = tdest.y - torg.y;
//.........这里部分代码省略.........

        /// <summary>
        /// Given the triangulation, and a vertex returns the minimum distance to the 
        /// vertices of the triangle where the given vertex located.
        /// </summary>
        /// <param name="newlocX"></param>
        /// <param name="newlocY"></param>
        /// <param name="searchtri"></param>
        /// <returns></returns>
        private double MinDistanceToNeighbor(double newlocX, double newlocY, ref Otri searchtri)
        {
            Otri horiz = default(Otri);	// for search operation
            LocateResult intersect = LocateResult.Outside;
            Vertex v1, v2, v3, torg, tdest;
            double d1, d2, d3, ahead;
            //triangle ptr;                         // Temporary variable used by sym().

            Point newvertex = new Point(newlocX, newlocY);

            // 	printf("newvertex %f,%f\n", newvertex[0], newvertex[1]);
            // Find the location of the vertex to be inserted.  Check if a good
            //   starting triangle has already been provided by the caller.	
            // Find a boundary triangle.
            //horiz.tri = m.dummytri;
            //horiz.orient = 0;
            //horiz.symself();
            // Search for a triangle containing 'newvertex'.
            // Start searching from the triangle provided by the caller.
            // Where are we?
            torg = searchtri.Org();
            tdest = searchtri.Dest();
            // Check the starting triangle's vertices.
            if ((torg.x == newvertex.x) && (torg.y == newvertex.y))
            {
                intersect = LocateResult.OnVertex;
                searchtri.Copy(ref horiz);

            }
            else if ((tdest.x == newvertex.x) && (tdest.y == newvertex.y))
            {
                searchtri.LnextSelf();
                intersect = LocateResult.OnVertex;
                searchtri.Copy(ref horiz);
            }
            else
            {
                // Orient 'searchtri' to fit the preconditions of calling preciselocate().
                ahead = Primitives.CounterClockwise(torg, tdest, newvertex);
                if (ahead < 0.0)
                {
                    // Turn around so that 'searchpoint' is to the left of the
                    //   edge specified by 'searchtri'.
                    searchtri.SymSelf();
                    searchtri.Copy(ref horiz);
                    intersect = mesh.locator.PreciseLocate(newvertex, ref horiz, false);
                }
                else if (ahead == 0.0)
                {
                    // Check if 'searchpoint' is between 'torg' and 'tdest'.
                    if (((torg.x < newvertex.x) == (newvertex.x < tdest.x)) &&
                        ((torg.y < newvertex.y) == (newvertex.y < tdest.y)))
                    {
                        intersect = LocateResult.OnEdge;
                        searchtri.Copy(ref horiz);

                    }
                }
                else
                {
                    searchtri.Copy(ref horiz);
                    intersect = mesh.locator.PreciseLocate(newvertex, ref horiz, false);
                }
            }
            if (intersect == LocateResult.OnVertex || intersect == LocateResult.Outside)
            {
                // set distance to 0
                //m.VertexDealloc(newvertex);
                return 0.0;
            }
            else
            { // intersect == ONEDGE || intersect == INTRIANGLE
                // find the triangle vertices
                v1 = horiz.Org();
                v2 = horiz.Dest();
                v3 = horiz.Apex();
                d1 = (v1.x - newvertex.x) * (v1.x - newvertex.x) + (v1.y - newvertex.y) * (v1.y - newvertex.y);
                d2 = (v2.x - newvertex.x) * (v2.x - newvertex.x) + (v2.y - newvertex.y) * (v2.y - newvertex.y);
                d3 = (v3.x - newvertex.x) * (v3.x - newvertex.x) + (v3.y - newvertex.y) * (v3.y - newvertex.y);
                //m.VertexDealloc(newvertex);
                // find minimum of the distance
                if (d1 <= d2 && d1 <= d3)
                {
                    return d1;
                }
                else if (d2 <= d3)
                {
                    return d2;
                }
                else
                {
                    return d3;
//.........这里部分代码省略.........