context: using System;
  using GeoAPI.Geometries;
  using GisSharpBlog.NetTopologySuite.Geometries;
  
  namespace GisSharpBlog.NetTopologySuite.Algorithm
  {
      /// <summary>
      /// Computes the minimum diameter of a <c>Geometry</c>.
      /// The minimum diameter is defined to be the
      /// width of the smallest band that contains the point,
      /// where a band is a strip of the plane defined
      /// by two parallel lines.
      /// This can be thought of as the smallest hole that the point can be
      /// moved through, with a single rotation.
      /// The first step in the algorithm is computing the convex hull of the Geometry.
      /// If the input Geometry is known to be convex, a hint can be supplied to
      /// avoid this computation.
      /// </summary>
      public class MinimumDiameter
      {
          private readonly IGeometry inputGeom;
          private readonly bool isConvex;
  
          private LineSegment minBaseSeg = new LineSegment();
          private ICoordinate minWidthPt = null;
          private int minPtIndex;
          private double minWidth = 0.0;
  
          /// <summary> 
          /// Compute a minimum diameter for a giver <c>Geometry</c>.
          /// </summary>
          /// <param name="inputGeom">a Geometry.</param>
          public MinimumDiameter(IGeometry inputGeom) 
              : this(inputGeom, false) { }
  
          /// <summary> 
          /// Compute a minimum diameter for a giver <c>Geometry</c>,
          /// with a hint if
          /// the Geometry is convex
          /// (e.g. a convex Polygon or LinearRing,
          /// or a two-point LineString, or a Point).
          /// </summary>
          /// <param name="inputGeom">a Geometry which is convex.</param>
          /// <param name="isConvex"><c>true</c> if the input point is convex.</param>
          public MinimumDiameter(IGeometry inputGeom, bool isConvex)
          {
              this.inputGeom = inputGeom;
              this.isConvex = isConvex;
          }
  
          /// <summary> 
          /// Gets the length of the minimum diameter of the input Geometry.
          /// </summary>
          /// <returns>The length of the minimum diameter.</returns>
          public double Length
          {
              get
              {
                  ComputeMinimumDiameter();
                  return minWidth;
              }
          }
  
          /// <summary>
          /// Gets the <c>Coordinate</c> forming one end of the minimum diameter.
          /// </summary>
          /// <returns>A coordinate forming one end of the minimum diameter.</returns>
          public ICoordinate WidthCoordinate
          {
              get
              {
                  ComputeMinimumDiameter();
                  return minWidthPt;
              }
          }
  
          /// <summary>
          /// Gets the segment forming the base of the minimum diameter.
          /// </summary>
          /// <returns>The segment forming the base of the minimum diameter.</returns>
          public ILineString SupportingSegment
          {
              get
              {
                  ComputeMinimumDiameter();
                  return inputGeom.Factory.CreateLineString(new ICoordinate[] { minBaseSeg.P0, minBaseSeg.P1 });
              }
          }
  
          /// <summary>
          /// Gets a <c>LineString</c> which is a minimum diameter.
          /// </summary>
          /// <returns>A <c>LineString</c> which is a minimum diameter.</returns>
          public ILineString Diameter
          {
              get
              {
                  ComputeMinimumDiameter();
  
                  // return empty linearRing if no minimum width calculated
                  if (minWidthPt == null)
                  {
                      ICoordinate[] nullCoords = null;
                      return inputGeom.Factory.CreateLineString(nullCoords);
                  }
  
                  ICoordinate basePt = minBaseSeg.Project(minWidthPt);
                  return inputGeom.Factory.CreateLineString(new ICoordinate[] { basePt, minWidthPt });
              }
          }
  
          /// <summary>
          /// 
          /// </summary>
          private void ComputeMinimumDiameter()
          {
              // check if computation is cached
              if (minWidthPt != null)
                  return;
  
              if (isConvex) ComputeWidthConvex(inputGeom);
              else
              {
                  IGeometry convexGeom = (new ConvexHull(inputGeom)).GetConvexHull();
                  ComputeWidthConvex(convexGeom);
              }
          }
  
          /// <summary>
          /// 
          /// </summary>
          /// <param name="geom"></param>
          private void ComputeWidthConvex(IGeometry geom)
          {
              ICoordinate[] pts = null;
              if (geom is IPolygon)
                   pts = ((IPolygon) geom).ExteriorRing.Coordinates;
              else pts = geom.Coordinates;
  
              // special cases for lines or points or degenerate rings
              if (pts.Length == 0) 
              {
                  minWidth = 0.0;
                  minWidthPt = null;
                  minBaseSeg = null;
              }
              else if (pts.Length == 1) 
              {
                  minWidth = 0.0;
                  minWidthPt = pts[0];
                  minBaseSeg.P0 = pts[0];
                  minBaseSeg.P1 = pts[0];
              }
              else if (pts.Length == 2 || pts.Length == 3) 
              {
                  minWidth = 0.0;
                  minWidthPt = pts[0];
                  minBaseSeg.P0 = pts[0];
                  minBaseSeg.P1 = pts[1];
              }
              else ComputeConvexRingMinDiameter(pts);
          }
  
          /// <summary> 
          /// Compute the width information for a ring of <c>Coordinate</c>s.
          /// Leaves the width information in the instance variables.
          /// </summary>
          /// <param name="pts"></param>
          private void ComputeConvexRingMinDiameter(ICoordinate[] pts)
          {
              // for each segment in the ring
              minWidth = Double.MaxValue;
              int currMaxIndex = 1;
  
              LineSegment seg = new LineSegment();
              // compute the max distance for all segments in the ring, and pick the minimum
              for (int i = 0; i < pts.Length - 1; i++) 
              {
                  seg.P0 = pts[i];
                  seg.P1 = pts[i + 1];
                  currMaxIndex = FindMaxPerpDistance(pts, seg, currMaxIndex);
              }
          }
  
          /// <summary>
          /// 
          /// </summary>
          /// <param name="pts"></param>
          /// <param name="seg"></param>
          /// <param name="startIndex"></param>
          /// <returns></returns>
          private int FindMaxPerpDistance(ICoordinate[] pts, LineSegment seg, int startIndex)
          {
              double maxPerpDistance = seg.DistancePerpendicular(pts[startIndex]);
              double nextPerpDistance = maxPerpDistance;
              int maxIndex = startIndex;
              int nextIndex = maxIndex;
              while (nextPerpDistance >= maxPerpDistance) 
              {
                  maxPerpDistance = nextPerpDistance;
                  maxIndex = nextIndex;
  
                  nextIndex = NextIndex(pts, maxIndex);
  
gt:                 nextPerpDistance = seg.DistancePerpendicular(pts[nextIndex]);