点到折线最短距离所在点距离折线起点的累积距离

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.Geometry;
using RGeos.Geometry

namespace RGeos.Geometry
{
public class CulmulateDistance
{
/// <summary>
/// 点到折线最短距离处距离折线起点的累积距离
/// </summary>
/// <param name="P">任意一点</param>
/// <param name="polyline">折线</param>
/// <returns></returns>
public static double CulmulateDist_Point_to_Polyline(IPoint P, IPolyline polyline)
{
ISegmentCollection segs = polyline as ISegmentCollection;
double min = double.MaxValue;
int segIndex = -1;//最短距离所在档的索引
for (int i = 0; i < segs.SegmentCount; i++)
{
//点到每条线段的最短距离
double dis = Dist_Point_to_Segment(P, segs.get_Segment(i));
if (dis < min)
{
min = dis;
segIndex = i;//取出最小的一个
}
}
double culmulateDis = 0;
for (int i = 0; i < segs.SegmentCount; i++)
{
if (segIndex != i)
{
culmulateDis += segs.get_Segment(i).Length;
}
else
{
ISegment current = segs.get_Segment(i);
Vector3d v = new Vector3d();
v.X = current.ToPoint.X - current.FromPoint.X;
v.Y = current.ToPoint.Y - current.FromPoint.Y;

Vector3d w = new Vector3d();
w.X = P.X - current.FromPoint.X;
w.Y = P.Y - current.FromPoint.Y;

double c1 = dot(w, v);//投影长度
if (c1 <= 0)//这种情况最短距离在该线段的起点处
{
break;
}
double c2 = dot(v, v);
if (c2 <= c1)//这种情况最短距离在该线段的终点处
{
culmulateDis += Math.Sqrt(c2);
break;
}

double b = c1 / c2;
culmulateDis += Math.Sqrt(c1);
IPoint Pb = new PointClass();
Pb.X = current.FromPoint.X + b * v.X;
Pb.Y = current.FromPoint.Y + b * v.Y;
break;
}
}
return culmulateDis;
}

public static IPoint GetCentrePoint(IPolygon polygon)
{
IArea pArea = polygon as IArea;
IPoint pt = new PointClass();
pt.X = pArea.Centroid.X;
pt.Y = pArea.Centroid.Y;
return pt;
}

public static double Dist_Point_to_Segment(IPoint P, ISegment S)
{
Vector3d v = new Vector3d();
v.X = S.ToPoint.X - S.FromPoint.X;
v.Y = S.ToPoint.Y - S.FromPoint.Y;

Vector3d w = new Vector3d();
w.X = P.X - S.FromPoint.X;
w.Y = P.Y - S.FromPoint.Y;

double c1 = dot(w, v);
if (c1 <= 0)
return d(P, S.FromPoint);

double c2 = dot(v, v);
if (c2 <= c1)
return d(P, S.ToPoint);

double b = c1 / c2;
IPoint Pb = new PointClass();
Pb.X = S.FromPoint.X + b * v.X;
Pb.Y = S.FromPoint.Y + b * v.Y;
return d(P, Pb);
}
/// <summary>
/// 向量的模
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public static double norm(Vector3d v)
{
return Math.Sqrt(dot2(v, v));  // norm = length of vector
}
/// <summary>
/// 2D数量积，点乘
/// </summary>
/// <param name="u"></param>
/// <param name="v"></param>
/// <returns></returns>
public static double dot2(Vector3d u, Vector3d v)
{
return ((u).X * (v).X + (u).Y * (v).Y);
}
public static double dot(Vector3d u, Vector3d v)
{
return ((u).X * (v).X + (u).Y * (v).Y + (u).Z * (v).Z);
}

public static double d(IPoint P, IPoint P1)
{
return Math.Sqrt((P1.X - P.X) * (P1.X - P.X) + (P1.Y - P.Y) * (P1.Y - P.Y));
}

/// <param name="x">增量X</param>
/// <param name="y">增量Y</param>
/// <returns>象限角</returns>
public static double GetQuadrantAngle(double x, double y)
{
double theta = Math.Atan(y / x);
if (x > 0 && y == 0) return 0;
if (x == 0 && y > 0) return Math.PI / 2;
if (x < 0 && y == 0) return Math.PI;
if (x == 0 && y < 0) return 3 * Math.PI / 2;

if (x > 0 && y > 0) return theta;
if (x > 0 && y < 0) return Math.PI * 2 + theta;
if (x < 0 && y > 0) return theta + Math.PI;
if (x < 0 && y < 0) return theta + Math.PI;
return theta;
}
}
}