/// <summary>
        /// 
        /// </summary>
        /// <param name="polygon"></param>
        /// <param name="writer"></param>
        public void Write(IPolygon polygon, BinaryWriter writer)
        {
            writer.Write((int) ShapeGeometryTypes.Polygon);

            // Write BoundingBox            
            WriteBoundingBox(polygon, writer);

            // Write NumParts and NumPoints            
            writer.Write((int) (polygon.NumInteriorRings + 1));
            writer.Write((int)  polygon.NumPoints);

            // Write IndexParts
            int count = 0;
            writer.Write((int) count);
            if (polygon.NumInteriorRings != 0)
            {
                // Write external shell index
                count += polygon.ExteriorRing.NumPoints;
                writer.Write((int) count);
                for (int i = 1; i < polygon.NumInteriorRings; i++)
                {
                    // Write internal holes index
                    count += polygon.GetInteriorRingN(i - 1).NumPoints;
                    writer.Write((int) count);
                }
            }

            // Write Coordinates
            for (int i = 0; i < polygon.NumPoints; i++)
                Write(polygon.Coordinates[i], writer);
        }

 /**
  * Finds a best-fit plane for the polygon, 
  * by sampling a few points from the exterior ring.
  * <p>
  * The algorithm used is Newell's algorithm:
  * - a base point for the plane is determined from the average of all vertices
  * - the normal vector is determined by
  *   computing the area of the projections on each of the axis planes
  * 
  * @param poly the polygon to determine the plane for
  * @return the best-fit plane
  */
 private static Plane3D FindBestFitPlane(IPolygon poly)
 {
     var seq = poly.ExteriorRing.CoordinateSequence;
     var basePt = AveragePoint(seq);
     var normal = AverageNormal(seq);
     return new Plane3D(normal, basePt);
 }

        public void HandleNewInput(IPolygon geometry)
        {
            // Reset labels
            lbNumVert2.Text = "-";
            lbNumTri2.Text = "-";
            lbNumSeg2.Text = "-";

            lbNumVert.Text = geometry.Points.Count.ToString();
            lbNumSeg.Text = geometry.Segments.Count().ToString();
            lbNumTri.Text = "0";

            // Statistics labels
            lbAreaMin.Text = "-";
            lbAreaMax.Text = "-";
            lbEdgeMin.Text = "-";
            lbEdgeMax.Text = "-";
            lbAngleMin.Text = "-";
            lbAngleMax.Text = "-";

            // Quality labels
            lbQualAlphaMin.Text = "-";
            lbQualAlphaAve.Text = "-";
            lbQualAspectMin.Text = "-";
            lbQualAspectAve.Text = "-";

            angleHistogram1.SetData(null, null);
        }

 public void HandleMeshImport(IPolygon geometry, Mesh mesh)
 {
     // Previous mesh stats
     lbNumVert2.Text = "-";
     lbNumTri2.Text = "-";
     lbNumSeg2.Text = "-";
 }

 protected override void OnRenderInternal(Map map, IPolygon polygon, IGraphics g)
 {
     IPoint pt = polygon.Centroid;
     Point renderingOrigin = Point.Truncate(Transform.WorldtoMap(pt.Coordinate, map));
     g.RenderingOrigin = new PointStruct(renderingOrigin.X, renderingOrigin.Y);
     base.OnRenderInternal(map, polygon, g);
 }

        // 0 0 0 0
        // 0 1 0 0
        // 0 1 2 0
        // 0 1 2 3 0 ...
        private void AppendCurvePoint(IPolygon polygon, ICoordinate worldPos)
        {
            List<ICoordinate> vertices = new List<ICoordinate>();

            ILineString linearRing = polygon.ExteriorRing;
            for (int i = 0; i < linearRing.Coordinates.Length; i++)
            {
                if (linearRing.Coordinates.Length <= 4)
                {
                    if (1 == i)
                    {
                        if (linearRing.Coordinates[0].Equals2D(linearRing.Coordinates[1]))
                        {
                            // 0 0 ? 0 -> 0 1 ? 0 
                            vertices.Add(worldPos);
                        }
                        else
                        {
                            // 0 1 ? 0 -> 0 1 ? 0
                            vertices.Add(linearRing.Coordinates[i]);
                        }
                    }
                    else if (2 == i)
                    {
                        if (linearRing.Coordinates[1].Equals2D(linearRing.Coordinates[2]))
                        {
                            // 0 0 0 0 -> 0 1 1 0
                            vertices.Add(worldPos);
                        }
                        else
                        {
                            // 0 1 2 0 -> 0 1 2 3 0
                            vertices.Add(linearRing.Coordinates[i]);
                            vertices.Add(worldPos);
                        }
                    }
                    else
                    {
                        vertices.Add(linearRing.Coordinates[i]);
                    }
                }
                else
                {
                    if (i == (linearRing.Coordinates.Length - 1))
                    {
                        // insert before last point to keep ring closed
                        vertices.Add(worldPos);
                    }
                    vertices.Add(linearRing.Coordinates[i]);
                }
            }
            int index = FeatureProvider.GetFeatureCount() - 1;
            ILinearRing newLinearRing = GeometryFactory.CreateLinearRing(vertices.ToArray());
            IPolygon newPolygon = GeometryFactory.CreatePolygon(newLinearRing, null);
            //##layerEditor.UpdateCurvePointInserted(index, newPolygon, vertices.Count - 1);
            //((FeatureProvider)layerEditor.VectorLayer.DataSource).UpdateGeometry(index, newPolygon);
            ((Feature)FeatureProvider.Features[index]).Geometry = newPolygon;
            Layer.RenderRequired = true;
            // do not remove see newline MapControl.SelectTool.Select((VectorLayer)Layer, newPolygon, -1);
        }

		/// <summary>
		/// Create a new Station object typing all the property values
		/// </summary>
		/// <param name="lat"></param>
		/// <param name="lon"></param>
		/// <param name="alt"></param>
		/// <param name="compass"></param>
		public Station (string stationID, IPolygon grid, double lat, double lon, double alt, double compass)
		{
			_lat=lat;
			_lon=lon;
			_alt=alt;
			_compass=compass;
			_stationID=stationID;
		}

 public static void ComputeDistance(IPolygon poly, Coordinate pt, PointPairDistance ptDist)
 {
     ComputeDistance(poly.ExteriorRing, pt, ptDist);
     for (var i = 0; i < poly.NumInteriorRings; i++)
     {
         ComputeDistance(poly.GetInteriorRingN(i), pt, ptDist);
     }
 }

 /// <summary>
 /// 
 /// </summary>
 /// <param name="p"></param>
 /// <returns></returns>
 private bool HasRepeatedPoint(IPolygon p)
 {
     if (HasRepeatedPoint(p.ExteriorRing.Coordinates))
         return true;
     for (int i = 0; i < p.NumInteriorRings; i++)
         if (HasRepeatedPoint(p.GetInteriorRingN(i).Coordinates)) 
             return true;            
     return false;
 }

 private void Benchmark(IPolygon poly)
 {
     var start = DateTime.Now;
     var valid = poly.IsValid;
     var end = DateTime.Now;
     var diff = end - start;
     var td = diff.TotalSeconds;
     LogLine("{0} \t{1} \t{2} \t{3}", poly.NumInteriorRings, poly.NumPoints, valid, td);
 }

        private static bool IsAllCW(IPolygon poly)
        {
            IList<Coordinate> shell = poly.Shell.Coordinates;
            for (int i = 0, c = shell.Count - 3; i < c; i++)
            {
                if (CgAlgorithms.ComputeOrientation(shell[i], shell[i + 1], shell[i + 2]) == 1)
                    return false;
            }

            return true;
        }

        public Triangulatable(IPolygon polygon)
        {
            List<Vector2D> vPoints = new List<Vector2D>(polygon.VerticesArray);
             List<TriangulationPoint> tPoints = vPoints.ConvertAll(p => (TriangulationPoint) (p));

             if (tPoints.Count < 3) throw new ArgumentException("List has fewer than 3 points", "points");

             // Lets do one sanity check that first and last point haven't got same position
             // It's something that often happens when importing polygon data from other formats
             if (tPoints[0].Equals(tPoints[tPoints.Count - 1])) tPoints.RemoveAt(tPoints.Count - 1);

             points.AddRange(tPoints);
        }

        /// <summary>
        /// Save a polygon geometry to disk.
        /// </summary>
        /// <param name="mesh">An instance of the <see cref="IPolygon" /> class.</param>
        /// <param name="filename">The path of the file to save.</param>
        public static void Write(IPolygon polygon, string filename)
        {
            foreach (IPolygonFormat format in formats)
            {
                if (format != null && format.IsSupported(filename))
                {
                    format.Write(polygon, filename);
                    return;
                }
            }

            throw new Exception("File format not supported.");
        }

 public static Boolean ContainsPointInPolygon(Coordinate p, IPolygon poly)
 {
     if (poly.IsEmpty) return false;
     ILinearRing shell = (ILinearRing)poly.ExteriorRing;
     if (!IsPointInRing(p, shell)) return false;
     // now test if the point lies in or on the holes
     for (int i = 0; i < poly.NumInteriorRings; i++)
     {
         ILinearRing hole = (ILinearRing)poly.GetInteriorRingN(i);
         if (IsPointInRing(p, hole)) return false;
     }
     return true;
 }

        /// <summary>
        /// Insert segments into the mesh.
        /// </summary>
        /// <param name="input">The polygon.</param>
        /// <param name="options">Constraint options.</param>
        public void Apply(IPolygon input, ConstraintOptions options)
        {
            behavior.Poly = input.Segments.Count > 0;

            // Copy constraint options
            if (options != null)
            {
                behavior.ConformingDelaunay = options.ConformingDelaunay;
                behavior.Convex = options.Convex;
                behavior.NoBisect = options.SegmentSplitting;

                if (behavior.ConformingDelaunay)
                {
                    behavior.Quality = true;
                }
            }

            //if (input.EdgeMarkers != null)
            //{
            //    behavior.UseBoundaryMarkers = true;
            //}

            behavior.useRegions = input.Regions.Count > 0;

            // Ensure that no vertex can be mistaken for a triangular bounding
            // box vertex in insertvertex().
            mesh.infvertex1 = null;
            mesh.infvertex2 = null;
            mesh.infvertex3 = null;

            if (behavior.useSegments)
            {
                // Segments will be introduced next.
                mesh.checksegments = true;

                // Insert PSLG segments and/or convex hull segments.
                FormSkeleton(input);
            }

            if (behavior.Poly && (mesh.triangles.Count > 0))
            {
                // Copy holes and regions
                mesh.holes.AddRange(input.Holes);
                mesh.regions.AddRange(input.Regions);

                // Carve out holes and concavities.
                CarveHoles();
            }
        }

        private void AddPolygon(IPolygon polygon, GraphicsPath graphicsPath)
        {
            if (polygon == null)
                throw new ArgumentNullException("polygon");

            ILineString exterior = polygon.ExteriorRing;
            IEnumerable<PointF> coords = this.GetCoords(exterior);
            graphicsPath.AddPolygon(coords.ToArray());

            foreach (ILineString ring in polygon.InteriorRings)
            {
                coords = this.GetCoords(ring);
                graphicsPath.AddPolygon(coords.ToArray());
            }
        }

        public void Render(LittleSharpRenderEngine engine, Graphics graphics, IPolygon polygon, IAreaStyle style)
        {
            if (polygon == null || style == null) return;

            System.Drawing.Drawing2D.GraphicsPath gp = new System.Drawing.Drawing2D.GraphicsPath();

            gp.AddPolygon(RenderUtil.CoordToPoint(polygon.Shell.Coordinates));
            foreach (ILinearRing l in polygon.Holes)
                gp.AddPolygon(RenderUtil.CoordToPoint(l.Coordinates));
            gp.CloseFigure();

            if (style.Fill != null) RenderUtil.RenderFill(engine, graphics, gp, style.Fill);

            if (style.Outline != null) RenderUtil.RenderOutline(engine, graphics, gp, style.Outline);
        }

        public static void AssertPolygonsEqual(IPolygon poly1, IPolygon poly2)
        {
            Assert.IsNotNull(poly1);
            Assert.IsNotNull(poly2);

            ILineString line1 = poly1.Shell;
            ILineString line2 = poly2.Shell;

            Assert.AreEqual(line1.Coordinates.Length, line2.Coordinates.Length, "Number of coordinates between polygons doesn't match");

            for (int i = 0; i < line2.Coordinates.Length; i++)
            {
                AssertCoordinatesEqual(line2.Coordinates[i], line1.Coordinates[i]);
            }
        }

        public static string ConvertPolygonToJson(IPolygon poly)
        {
            // Note: According to the specification for GeoJSON, the CRS tags are deprecated and no
            // longer used. See: https://tools.ietf.org/html/draft-butler-geojson-06#page-9
            // GeoJSON is supposed to be defaulted and assumed to always be projected
            // as CRS 84 / EPSG:4326

            // with the above in mind, should we reproject the polygon?

            // Project polygon to GeoJSON default EPSG:4326
            //ISpatialReferenceFactory2 sre = (ISpatialReferenceFactory2)new SpatialReferenceEnvironment();
            //ISpatialReference sr = sre.CreateSpatialReference(4326);
            //poly.Project(sr);
            //poly.SpatialReference = sr;

            JsonPolygon jsonPoly = new JsonPolygon();

            IPointCollection pc = (IPointCollection)poly;
            IEnvelope bbox = poly.Envelope;

            List<List<double>> allPoints = new List<List<double>>();

            for (int i = 0; i < pc.PointCount; i++)
            {
                double x = pc.get_Point(i).X;
                double y = pc.get_Point(i).Y;

                List<double> point = new List<double>();
                point.Add(x);
                point.Add(y);

                allPoints.Add(point);
            }

            // the GeoJSON spec has the coordinate array a three-array stagger (points in arrays, holes in array around that, multipoly in another array)
            // we don't deal with multipoly or holes in SUITT, but we need to make sure the format matches, hence the wierd List
            jsonPoly.coordinates.Add(allPoints);

            //define bbox: minx, miny, maxx, maxy
            jsonPoly.bbox.Add(bbox.XMin);
            jsonPoly.bbox.Add(bbox.YMin);
            jsonPoly.bbox.Add(bbox.XMax);
            jsonPoly.bbox.Add(bbox.YMax);

            // set the projection
            //jsonPoly.crs = new GeoJsonCRS("EPSG:" + poly.SpatialReference.FactoryCode.ToString());
            return JsonConvert.SerializeObject(jsonPoly);
        }

        public static IMapPoint Center(IPolygon polygon)
        {
            IEnumerable<IMapPoint> pc = polygon.GetPoints();
            double x = 0;
            double y = 0;

            foreach (IMapPoint p in pc)
            {
                x += p.X;
                y += p.Y;
            }
            x /= pc.Count();
            y /= pc.Count();
            IMapPoint center = Runtime.geometryEngine.newMapPoint(x, y);
            return center;
        }