离散点最大边界算法

还记得当初找工作时的第一期望就是能找一个和算法相关的工作,可是事与愿违,被拉进来一家gis方面的软件公司,既来之则安之,工作了半年多,终于也遇见了一个算法,最优路径,当初在学校是曾经稍微写过一点,然而postgis又提供了这样一个函数,直接拿来用就好了,自己也就没太走心.如今,又有一个新的需求,爆管分析,就是有一个管点出现问题了,我能用postgis函数找到它周围所有相连的管点,再从这些管点中找出它们的最大范围.简而言之,就是求离散点的最大边界.

定义离散点类

public class point
{
    public Nullable<decimal> x { get; set; }
    public Nullable<decimal> y { get; set; }
    /**
     * 边界查找算法中 是否被找到
     */
    public bool founded = false;
    public point(Nullable<decimal> x, Nullable<decimal> y)
    {
        this.x = x;
        this.y = y;
    }
    public bool equals(point other)
    {
        if (this.x != other.x)
            return false;
        if (this.y != other.y)
            return false;
        return true;
    }
}

离散点工具类

public class Utils
{
/**
 * <p>
 * <b>求两点之间的长度</b>
 * <p>
 * <pre>
 * 求两点之间的长度
 * </pre>
 * @param   ws  西南点
 * @param   en  东北点
 * @return  两点之间的长度
 */
public static double calLineLen(point ws, point en) {
    if (null == ws || null == en) {
        return .0;
    }
    if (ws.equals(en)) {
        return .0;
    }
    double a = Math.Abs(Convert.ToDouble(ws.x) - Convert.ToDouble(en.x)); // 直角三角形的直边a
    double b = Math.Abs(Convert.ToDouble(ws.y) - Convert.ToDouble(en.y)); // 直角三角形的直边b
    double min = Math.Min(a, b); // 短直边
    double max = Math.Max(a, b); // 长直边
    /**
     * 为防止计算平方时float溢出，做如下转换
     * √(min²+max²) = √((min/max)²+1) * abs(max)
     */
    double inner = min / max;
    return Math.Sqrt(inner * inner + 1.0) * max;
}

/**
 * <p>
 * <b>计算向量角</b>
 * <p>
 * <pre>
 * 计算两点组成的向量与x轴正方向的向量角
 * </pre>
 * @param   s   向量起点
 * @param   d   向量终点
 * @return  向量角
 */
public static double angleOf(point s, point d) {
    double dist = calLineLen(s, d);
    if (dist <= 0) {
        return -1;
    }
    double x = Convert.ToDouble(d.x - s.x); // 直角三角形的直边a
    double y = Convert.ToDouble(d.y - s.y); // 直角三角形的直边b
    if (y >= 0) { /* 1 2 象限 */
        return Math.Acos(x / dist);
    } else { /* 3 4 象限 */
        return Math.Acos(-x / dist) + Math.PI;
    }
}
/**
 * <p>
 * <b>修正角度</b>
 * <p>
 * <pre>
 * 修正角度到 [0, 2PI]
 * </pre>
 *
 * @author  ManerFan 2015年4月17日
 * @param   angle   原始角度
 * @return  修正后的角度
 */
public static double reviseAngle(double angle) {
    while (angle < 0) {
        angle += 2 * Math.PI;
    }
    while (angle >= 2 * Math.PI) {
        angle -= 2 * Math.PI;
    }
    return angle;
}
}

寻找起始点坐标

/** 查找起始点（保证y最大的情况下、尽量使x最小的点） */
 private static point findStartpoint(List<point> ps)
    {
        if (null == ps || ps.Count() == 0)
        {
            return null;
        }
        point p = ps[0];
        for (int i = 0; i < ps.Count(); i++)
        {
            if (ps[i].y > p.y || (ps[i].y == p.y && ps[i].x < p.x))
            { /* 找到最靠上靠左的点 */
                p = ps[i];
            }
        }
        return p;
    }

最大边界算法

public List<point> findSmallestPolygon(List<point> ps)
   {
       JsonResult result = new JsonResult();
       //ps = test();  //测试数据
       if (null == ps || ps.Count()==0)
       {
           //ps = test();  //测试数据
           return null;
       }
       point corner = findStartpoint(ps);
       if (null == corner)
       {
           return null;
       }
       double minAngleDif, oldAngle = 2 * Math.PI;
       List<point> bound = new List<point>();
       do
       {
           minAngleDif = 2 * Math.PI;
           bound.Add(corner);
           point nextpoint = corner;
           double nextAngle = oldAngle;
           foreach (point p in ps)
           {
               if (p.founded)
               { // 已被加入边界链表的点
                   continue;
               }
               if (p.equals(corner))
               { // 重合点
                   /*if (!p.equals(bound.getFirst())) {
                       p.founded = true;
                   }*/
                   continue;
               }
               double currAngle = Utils.angleOf(corner, p); /* 当前向量与x轴正方向的夹角 */
               double angleDif = Utils.reviseAngle(oldAngle - currAngle); /* 两条向量之间的夹角（顺时针旋转的夹角） */
               if (angleDif < minAngleDif)
               {
                   minAngleDif = angleDif;
                   nextpoint = p;
                   nextAngle = currAngle;
               }
           }
           oldAngle = nextAngle;
           corner = nextpoint;
           corner.founded = true;
       } while (!corner.equals(bound[0])); /* 判断边界是否闭合 */
       //result.Data = bound;
       //result.JsonRequestBehavior = JsonRequestBehavior.AllowGet;
       //return result;
       return bound;
   }

测试数据

public static List<point> test() {
        List<point> ps = new List<point>();
        ps.Add(new point(1, 1));
        ps.Add(new point(1, 2));
        ps.Add(new point(2, 3));
        ps.Add(new point(2, 4));
        ps.Add(new point(3, 5));
        ps.Add(new point(4, 4));
        ps.Add(new point(5, 4));
        ps.Add(new point(5, 3));
        ps.Add(new point(4, 2));
        ps.Add(new point(4, 1)); 
        ps.Add(new point(3, 1));
        ps.Add(new point(2, 2));
        // 以上为所求结果
        ps.Add(new point(3, 2));
        ps.Add(new point(3, 3));
        ps.Add(new point(4, 3));
        ps.Add(new point(3, 4));
        return ps;
    }