求任意多边形内部水平方向似最大矩形算法实现

背景说明

前段时间有个求点是否在多边形内部的需求，折腾了不少时间，现截取其中的的重点部分——求任意多边形内部水平方向似最大矩形——来搞篇博客。

求点是否在多边形内部这个算法很容易搞，一搜一大把，但数据量大的时候，算法就必须进行优化。一个显然的优化点就是求最大内接矩形，毕竟判断点是否在矩形内，最多只需要执行四个判断语句，执行速度非常快；而要判断多边形，则需要与每条边比较，相对于矩形会慢很多，特别是在做GIS数据的时候，基本全是复杂多边形。

原算法在这：https://www.cnblogs.com/naaoveGIS/p/7218634.html
虽然原始文章写了算法实现这个章节，但实际上是个空壳，真正实现还得靠自己，以下进行详细说明。

算法流程

既然要实现，这里先罗列一下原算法流程。

• 第一步：获取任意多边形的四角坐标，通过四角坐标构造矩形，将该矩形划分成M*N个规则格网；
• 第二步：遍历所有格网，判断每个格网和多边形的包含关系。格网在多边形中，则标记为1，否则为0；
• 第三步：计算由0和1组成的矩形中，由1组成的最大矩形；
• 第四步：求得所得最大矩形代表的四角坐标，构造成真实地理矩形。

原作者提了两个注意点，一是第三步比较难，二是第一步的矩形划分粒度根据需求精度来。

第三步确实挺难，但是也基本是一搜一大把，这里贴上leetcode 85.最大矩形，这题是求最大矩形面积，做一点小小的改动即可得到最大矩形的四角坐标。然后第二点的话，需要根据实际项目需求以及硬件配置来选择精度。

算法实现

注：本文使用groovy实现，这玩意与java无缝对接，同时又可以当脚本语言用，好使。

首先定义几个实体类，分别是点、矩形、多边形：

class Point {
    /**
     * 经度.
     */
    def lon
    /**
     * 纬度.
     */
    def lat
}

class Rectangle {
    /**
     * 最小经度.
     */
    def minLon
    /**
     * 最大经度.
     */
    def maxLon
    /**
     * 最小纬度.
     */
    def minLat
    /**
     * 最大纬度.
     */
    def maxLat
}

class MultiPolygon {
    /**
     * 多边形最大内接矩形.
     */
    List<Rectangle> largestEnclosingRectangle
    /**
     * 多边形数据点.
     */
    List<Point> points
}

多边形最大内接矩形搁成List是因为后续可能需要多个内接矩形，提高效率。

第一步：划分矩形

划分矩形是先求多边形的最小外接矩形，然后对这个最小外接矩形进行划分，求最小外接矩形可太简单了，就是遍历多边形的每个点，求得最大/最小经纬度即可，代码如下：

    /**
     * 获取多边形区域的最小外接矩形.
     */
    static Rectangle computeMinimumEnclosingRectangle(List<Point> pts) {
        def rec = new Rectangle(
                minLon: Double.MAX_VALUE, maxLon: Double.MIN_VALUE,
                minLat: Double.MAX_VALUE, maxLat: Double.MIN_VALUE)
        pts.each { p ->
            rec.minLon = [p.lon, rec.minLon].min()
            rec.maxLon = [p.lon, rec.maxLon].max()
            rec.minLat = [p.lat, rec.minLat].min()
            rec.maxLat = [p.lat, rec.maxLat].max()
        }
        rec
    }

求得外接最小矩形之后，为了数据好看一点，先将其坐标按精度标准化。假设你准备划分的小矩形精度为小数点后两位，那么最好把这个外接最小矩形的坐标搞成小数点后两位。思路很简单，小的往更小去，大的往更大去即可，代码如下：

 //1.根据切分尺度标准化矩形
def minLon = floor(rec.minLon, scale), maxLon = ceil(rec.maxLon, scale),
    minLat = floor(rec.minLat, scale), maxLat = ceil(rec.maxLat, scale)

private static def ceil(d, scale) {
    def n = (1 / scale) as int
    (BigDecimal) Math.ceil(d * n) / n
}

private static def floor(d, scale) {
    def n = (1 / scale) as int
    (BigDecimal) Math.floor(d * n) / n
}

实测使用double的话，会出现奇奇怪怪的精度问题，比如25.1+0.01=25.110000000000003，25.2-0.01=25.189999999999998，实测了C、python、java、groovy，皆会出现这种问题，解决起来不难，根据需求精度加ceil和floor函数即可。但是还不如直接用BigDecimal，只要把内存管好就行。

在得到标准化的外接最小矩形后，直接按精度切它。注意，我这里切完之后，保存的并不是一个个矩形，而是按顺序来的一个个数据点，因为矩形除边缘外，上下左右都是连续的，存矩形的话，会导致后面的计算把一个数据点计算4次，很不划算。这里还是说一下数据点矩阵怎么转换得到矩形矩阵吧。

先直观一点看一个被切分的矩形：

假设存矩形数据的矩阵是rec[][]，存数据点的矩阵是dot[][]。

看上图左上角那个矩形①，它是rec[0][0]；再看最左上角那个数据点，它是dot[0][0]，再自己随便找几个点看看，会发现矩阵在rec中的坐标[i][j]等于其左上角的点在dot中的坐标，知道左上角就自然知道整个矩形的数据点了。

单独看一个数据点dot[i][j]，则它的四个方向的矩形（如果存在）坐标应该如下（假设中间那点为dot[i][j]）：

代码如下，基本就是一个等距划分。其实也可以不存储划分的数据点，直接等距操作，不过就是操作有点麻烦，可以省一点内存。

/**
 * 根据矩形和切分尺度获取切分矩阵的点阵.
 * 为了减少计算量，这里输出的并不是真的矩形，由于是等距连续切分，直接输出切分的数据点即可，以左上角数据点作为标的.
 */
private static List<List<Point>> computeRegionalSlices(Rectangle rec, scale) {
    //1.根据切分尺度标准化矩形
    def minLon = floor(rec.minLon, scale), maxLon = ceil(rec.maxLon, scale),
        minLat = floor(rec.minLat, scale), maxLat = ceil(rec.maxLat, scale)

    //2.切分(出于人类习惯，暂定从上往下按行切分，即按纬度从大到小，经度从小到大)
    List<List<Point>> matrix = []
    for (def posLat = maxLat; posLat >= minLat; posLat -= scale) {
        List<Point> row = []
        for (def posLon = minLon; posLon <= maxLon; posLon += scale) {
            row << new Point(lon: posLon, lat: posLat)
        }
        matrix << row
    }
    matrix
}

这里上几张实际切分图。
整体：

局部细节：

第二步：遍历网格，进行标记，在多边形内标1，不在标0

这步比较简单，根据上面的一个点对应的四个矩形，进行标记即可。代码如下：

/**
 * 根据切分矩形和多边形计算出标记矩阵.
 */
private static int[][] computeMarkMatrix(List<Point> pts, List<List<Point>> regionalSlices) {
    def m = regionalSlices.size(), n = regionalSlices[0].size(),
        rectangleMarks = new int[m - 1][n - 1]

    //先将矩形标记矩阵全部标记为1，再根据点阵计算，将不在多边形内的矩阵标记为0
    rectangleMarks.each { Arrays.fill(it, 1) }

    //遍历切分点阵的每一个点，得到标记矩阵
    def inRange = { num, min, max -> num >= min && num <= max }
    (0..<m).each { posM ->
        (0..<n).each { posN ->
            def p = regionalSlices[posM][posN]
            //处理不在多边形内部的点
            if (!isPolygonContainsPoint(pts, p)) {
                //该点不在多边形内部，处理该点对应的四个方向的矩阵，行范围是[0, m-2]，列范围时[0, n-2]
                //左上角[posM-1, posN-1]
                if (inRange(posM - 1, 0, m - 2) && inRange(posN - 1, 0, n - 2)) {
                    rectangleMarks[posM - 1][posN - 1] = 0
                }
                //右上角[posM-1, posN]
                if (inRange(posM - 1, 0, m - 2) && inRange(posN, 0, n - 2)) {
                    rectangleMarks[posM - 1][posN] = 0
                }
                //左下角[posM, posN1-1]
                if (inRange(posM, 0, m - 2) && inRange(posN - 1, 0, n - 2)) {
                    rectangleMarks[posM][posN - 1] = 0
                }
                //右下角[posM, posN]
                if (inRange(posM, 0, m - 2) && inRange(posN, 0, n - 2)) {
                    rectangleMarks[posM][posN] = 0
                }
            }
        }
    }
    rectangleMarks
}

/**
 * 返回一个点是否在一个多边形区域内(开区间).
 */
static Boolean isPolygonContainsPoint(List<Point> pts, Point p) {
    int nCross = 0
    pts.size().times { i ->
        //算法思路：取多边形任意一个边,做点point的水平延长线,求解与当前边的交点个数，这里只求解点右边的交点
        def p1 = pts[i], p2 = pts[(i + 1) % pts.size()]
        //p1p2是水平线段,要么没有交点,要么有无限个交点
        if (p1.lat == p2.lat) {
            return
        }
        //point在p1p2 底部或者顶部 --> 无交点
        if (p.lat < [p1.lat, p2.lat].min() || p.lat >= [p1.lat, p2.lat].max()) {
            return
        }
        //求解point点水平线与当前p1p2边的交点的X坐标
        double x = p1.lon + (p2.lon - p1.lon) * (p.lat - p1.lat) / (p2.lat - p1.lat)
        //只记录单边交点，这里只记录点右边的交点
        if (x > p.lon) {
            ++nCross
        }
    }
    //交点奇数，点在多边形内；反则反之
    nCross % 2 == 1
}

注意，这一步返回的是01矩阵矩阵，因为下一步计算最大矩形如果使用数据点进行计算，会比较麻烦。

第三步：计算最大矩形

这一步直接去看上面的LeetCode算法就好，这里只加了几行，记录最大矩阵的坐标：

/**
 * 根据标记矩阵求最大矩形，返回[最小行标 最大行标 最小列标 最大列标 最大面积]
 */
static def maximalRectangle(int[][] matrix) {
    def m = matrix.length, n = matrix[0].length
    int[][] left = new int[m][n]
    (0..<m).each { i ->
        (0..<n).each { j ->
            if (matrix[i][j] == 1) {
                left[i][j] = (!j ? 0 : left[i][j - 1]) + 1
            }
        }
    }

    def minC = -1, maxC = -1, minR = -1, maxR = -1, ret = 0
    (0..<n).each { j ->  // 对于每一列，使用基于柱状图的方法
        int[] up = new int[m]
        int[] down = new int[m]

        Deque<Integer> stack = new LinkedList<>()   //链栈
        for (int i = 0; i < m; i++) {
            while (!stack.isEmpty() && left[stack.peek()][j] >= left[i][j]) {
                stack.pop()
            }
            up[i] = stack.isEmpty() ? -1 : stack.peek()
            stack.push(i)
        }
        stack.clear()
        for (int i = m - 1; i >= 0; i--) {
            while (!stack.isEmpty() && left[stack.peek()][j] >= left[i][j]) {
                stack.pop()
            }
            down[i] = stack.isEmpty() ? m : stack.peek()
            stack.push(i)
        }
        for (int i = 0; i < m; i++) {
            int height = down[i] - up[i] - 1
            int area = height * left[i][j]
            //记录最大矩形的位置
            if (area > ret) {
                ret = area
                minC = up[i] + 1
                maxC = down[i] - 1
                minR = j - left[i][j] + 1
                maxR = j
            }
        }
    }
    return [minC, maxC, minR, maxR, ret]
}

由于groovy是可以返回多值的，所以上面的代码直接返回了该最大矩形在矩形矩阵中的最小/最大行/列标记，同时返回了该矩形的面积。

第四步：求得实际坐标

这步就非常容易了，自己根据上面的几个图稍微推导一下即可，就直接一行代码：

new Rectangle(minLon: regionalSlices[0][minR].lon, maxLon: regionalSlices[0][maxR + 1].lon,
        minLat: regionalSlices[maxC + 1][0].lat, maxLat: regionalSlices[minC][0].lat)

下面放几张效果图。

中间最大的那个是最大内接矩形，然后我算了算，这个图的话，这个最大矩形覆盖面积只有56%左右，所以可以按照面积需求，或者其他需求继续切下去。

完整代码放到下面，里面还有一个扩展多边形的算法，也是我这需求之一，这里就不进行详细阐述了。

package com.fiberhome.utils.position

import com.fiberhome.utils.position.base.Point
import com.fiberhome.utils.position.base.Rectangle
import com.fiberhome.utils.position.base.ScopeOfCity
import groovy.util.logging.Slf4j

import java.text.DecimalFormat

/**
 * created with IntelliJ IDEA 2020.3
 * author: hxw
 * date: 2021/8/12 9:46
 * version: 1.0
 * description:
 *      位置工具类.
 */
@Slf4j
class PositionUtils {

    /**
     * 多边形扩展算法，算法参考：https://blog.csdn.net/lweiyue/article/details/103033630
     *
     * @param pts 需要扩展的多边形，要求边数不少于3
     * @param expand 需要扩展的距离
     */
    static List<Point> polygonExpand(List<Point> pts, double expand) {
        def ans = []
        def norm = { x, y -> Math.sqrt(x * x + y * y) }
        def df = { x ->
            def df = new DecimalFormat('#.######')
            df.format(x) as double
        }
        pts.eachWithIndex { p, index ->
            def p1 = pts[index - 1], p2 = pts[(index + 1) % pts.size()]

            double v1x = p1.lon - p.lon
            double v1y = p1.lat - p.lat
            double n1 = norm(v1x, v1y)
            double vv1x = v1x / n1
            double vv1y = v1y / n1

            double v2x = p2.lon - p.lon
            double v2y = p2.lat - p.lat
            double n2 = norm(v2x, v2y)
            double vv2x = v2x / n2
            double vv2y = v2y / n2

            double judge = v1x * v2y - v2x * v1y
            double vectorLen = -expand / Math.sqrt((1 - (vv1x * vv2x + vv1y * vv2y)) / 2.0f)
            if (judge > 0)
                vectorLen *= -1
            double vx = vv1x + vv2x
            double vy = vv1y + vv2y
            vectorLen = vectorLen / norm(vx, vy)
            vx *= vectorLen
            vy *= vectorLen

            ans << new Point(lon: df.call(vx + p.lon), lat: df.call(vy + p.lat))
        }
        ans
    }

    /**
     * 返回一个点是否在一个多边形区域内(开区间).
     */
    static Boolean isPolygonContainsPoint(List<Point> pts, Point p) {
        int nCross = 0
        pts.size().times { i ->
            //算法思路：取多边形任意一个边,做点point的水平延长线,求解与当前边的交点个数，这里只求解点右边的交点
            def p1 = pts[i], p2 = pts[(i + 1) % pts.size()]
            //p1p2是水平线段,要么没有交点,要么有无限个交点
            if (p1.lat == p2.lat) {
                return
            }
            //point在p1p2 底部或者顶部 --> 无交点
            if (p.lat < [p1.lat, p2.lat].min() || p.lat >= [p1.lat, p2.lat].max()) {
                return
            }
            //求解point点水平线与当前p1p2边的交点的X坐标
            double x = p1.lon + (p2.lon - p1.lon) * (p.lat - p1.lat) / (p2.lat - p1.lat)
            //只记录单边交点，这里只记录点右边的交点
            if (x > p.lon) {
                ++nCross
            }
        }
        //交点奇数，点在多边形内；反则反之
        nCross % 2 == 1
    }

    /**
     * 求一个多边形区域的水平方向最大内接矩形，由于是经纬度数据，精确到小数点后两位，误差(只小不大)约一公里.
     * [算法参考]
     *      算法流程：https://www.cnblogs.com/naaoveGIS/p/7218634.html
     *      最大矩形：https://leetcode-cn.com/problems/maximal-rectangle/solution/zui-da-ju-xing-by-leetcode-solution-bjlu/
     *
     * 算法步骤：
     *      1.根据最小外接矩形进行区域切块，这里以0.01作为切分点，经纬度每隔0.01切分一次
     *      2.检查该区域是否在多边形内部，true标记1，false标记0，得到标记矩阵
     *      3.根据上述标记矩形根据leetcode算法获取最大内接矩阵，注意空间复杂度为O(mn)
     */
    static Rectangle computeLargestEnclosingRectangle(List<Point> pts) {
        //1.区域切块，不是真的切成矩形，而是切分成数据点，减少75%的计算量
        def scale = 0.01
        def minimumEnclosingRectangle = computeMinimumEnclosingRectangle(pts)
        def regionalSlices = computeRegionalSlices(minimumEnclosingRectangle, scale)

        //2.标记矩阵，这里将点阵经纬度转换为矩形标记矩阵，每个矩形以左上角作为标的，
        //  比如矩形marks[0][0]的左上角坐标为regionalSlices[0][0]，右下角坐标为regionalSlices[1][1]
        def marks = computeMarkMatrix(pts, regionalSlices)

        //3.计算最大内接矩阵，获取矩形
        def (minC, maxC, minR, maxR, area) = maximalRectangle(marks)
        new Rectangle(minLon: regionalSlices[0][minR].lon, maxLon: regionalSlices[0][maxR + 1].lon,
                minLat: regionalSlices[maxC + 1][0].lat, maxLat: regionalSlices[minC][0].lat)
    }

    /**
     * 递归求得最最大矩形.
     * 这里最难得到的是停机条件，粗略计算停机条件如下.
     *
     * 设外接矩形面积为a1，内接矩形面积为a4，多边形边数为n1，内接矩形数量为n2，判断一条边耗时t1，判断一个内接矩形耗时t2，
     * 由程序执行语句得到t1约是t2的三倍，则总耗时
     *          t = n2 * a4 / a1 + (n2 + 3 * n1)* (1 - a4 / a1)
     * 停机条件:
     *      本次迭代的t大于上次迭代的t，则停机，取上次数据为最终结果.
     */
    static List<Rectangle> computeLargestEnclosingRectangleRecursion(List<Point> pts) {
        //1.区域切块，不是真的切成矩形，而是切分成数据点，减少75%的计算量
        def scale = 0.01
        def minimumEnclosingRectangle = computeMinimumEnclosingRectangle(pts)
        def regionalSlices = computeRegionalSlices(minimumEnclosingRectangle, scale)

        //2.标记矩阵，这里将点阵经纬度转换为矩形标记矩阵，每个矩形以左上角作为标的，
        //  比如矩形marks[0][0]的左上角坐标为regionalSlices[0][0]，右下角坐标为regionalSlices[1][1]
        def marks = computeMarkMatrix(pts, regionalSlices)

        //3.迭代计算最优解
        def ret = []
        //计算a1、n1，这两个值不会随着迭代次数改变而改变
        def a1 = marks.length * marks[0].length, n1 = pts.size(),
            a4 = 0, n2 = 0

        def ot, nt = Double.MAX_VALUE
        do {
            ot = nt
            //计算a4、n2，这两个值每次迭代都会改变
            def (minC, maxC, minR, maxR, area) = maximalRectangle(marks)
            a4 += area
            ++n2
            nt = n2 * a4 / a1 + (n2 + 3.0 * n1) * (1.0 - a4 / a1)
            minC.upto(maxC) { i ->
                minR.upto(maxR) { j -> marks[i][j] = 0 }
            }
            ret << new Rectangle(minLon: regionalSlices[0][minR].lon, maxLon: regionalSlices[0][maxR + 1].lon,
                    minLat: regionalSlices[maxC + 1][0].lat, maxLat: regionalSlices[minC][0].lat)
        } while (nt < ot)
        ret.removeLast()
        ret
    }

    /**
     * 判断一个点是否在矩形内.
     */
    static boolean isRectangleContainsPoint(Rectangle rec, Point p) {
        return p.lon >= rec.minLon && p.lon <= rec.maxLon
                && p.lat >= rec.minLat && p.lat <= rec.maxLat
    }

    /**
     * 获取多边形区域的最小外接矩形.
     */
    static Rectangle computeMinimumEnclosingRectangle(List<Point> pts) {
        def rec = new Rectangle(
                minLon: Double.MAX_VALUE, maxLon: Double.MIN_VALUE,
                minLat: Double.MAX_VALUE, maxLat: Double.MIN_VALUE)
        pts.each { p ->
            rec.minLon = [p.lon, rec.minLon].min()
            rec.maxLon = [p.lon, rec.maxLon].max()
            rec.minLat = [p.lat, rec.minLat].min()
            rec.maxLat = [p.lat, rec.maxLat].max()
        }
        rec
    }

    /**
     * 根据中心点计算该省份扩展一公里所需要的公里数.
     */
    static def computeExpandDistance(Point p) {
        def lon = p.lon, lat = p.lat, scale = 0.0000001
        do {
            lat += scale
        } while (dist(p, new Point(lon: lon, lat: lat)) < 1000.0)
        return lat - p.lat
    }

    /**
     * 判断点是否在城市里.
     */
    static boolean isCityContainsPoint(ScopeOfCity soc, Point p) {
        //先判断外接矩形
        if (!isRectangleContainsPoint(soc.minimumEnclosingRectangle, p)) {
            return false
        }
        //判断内接矩形
        soc.regions.find { mp ->
            mp.largestEnclosingRectangle.find { rec -> isRectangleContainsPoint(rec, p) }
        }.asBoolean() ?:
                //判断多边形
                soc.regions.find { mp -> isPolygonContainsPoint(mp.points, p) }.asBoolean()

    }

    /*---------------------------------------------------------------------------------------------------------------*/

    /**
     * 根据矩形和切分尺度获取切分矩阵的点阵.
     * 为了减少计算量，这里输出的并不是真的矩形，由于是等距连续切分，直接输出切分的数据点即可，以左上角数据点作为标的.
     */
    private static List<List<Point>> computeRegionalSlices(Rectangle rec, scale) {
        //1.根据切分尺度标准化矩形
        def minLon = floor(rec.minLon, scale), maxLon = ceil(rec.maxLon, scale),
            minLat = floor(rec.minLat, scale), maxLat = ceil(rec.maxLat, scale)

        //2.切分(出于人类习惯，暂定从上往下按行切分，即按纬度从大到小，经度从小到大)
        List<List<Point>> matrix = []
        for (def posLat = maxLat; posLat >= minLat; posLat -= scale) {
            List<Point> row = []
            for (def posLon = minLon; posLon <= maxLon; posLon += scale) {
                row << new Point(lon: posLon, lat: posLat)
            }
            matrix << row
        }
        matrix
    }

    /**
     * 根据切分矩形和多边形计算出标记矩阵.
     */
    private static int[][] computeMarkMatrix(List<Point> pts, List<List<Point>> regionalSlices) {
        def m = regionalSlices.size(), n = regionalSlices[0].size(),
            rectangleMarks = new int[m - 1][n - 1]

        //先将矩形标记矩阵全部标记为1，再根据点阵计算，将不在多边形内的矩阵标记为0
        rectangleMarks.each { Arrays.fill(it, 1) }

        //遍历切分点阵的每一个点，得到标记矩阵
        def inRange = { num, min, max -> num >= min && num <= max }
        (0..<m).each { posM ->
            (0..<n).each { posN ->
                def p = regionalSlices[posM][posN]
                //处理不在多边形内部的点
                if (!isPolygonContainsPoint(pts, p)) {
                    //该点不在多边形内部，处理该点对应的四个方向的矩阵，行范围是[0, m-2]，列范围时[0, n-2]
                    //左上角[posM-1, posN-1]
                    if (inRange(posM - 1, 0, m - 2) && inRange(posN - 1, 0, n - 2)) {
                        rectangleMarks[posM - 1][posN - 1] = 0
                    }
                    //右上角[posM-1, posN]
                    if (inRange(posM - 1, 0, m - 2) && inRange(posN, 0, n - 2)) {
                        rectangleMarks[posM - 1][posN] = 0
                    }
                    //左下角[posM, posN1-1]
                    if (inRange(posM, 0, m - 2) && inRange(posN - 1, 0, n - 2)) {
                        rectangleMarks[posM][posN - 1] = 0
                    }
                    //右下角[posM, posN]
                    if (inRange(posM, 0, m - 2) && inRange(posN, 0, n - 2)) {
                        rectangleMarks[posM][posN] = 0
                    }
                }
            }
        }
        rectangleMarks
    }

    /**
     * 根据标记矩阵求最大矩形，返回[最小行标 最大行标 最小列标 最大列标 最大面积]
     */
    static def maximalRectangle(int[][] matrix) {
        def m = matrix.length, n = matrix[0].length
        int[][] left = new int[m][n]
        (0..<m).each { i ->
            (0..<n).each { j ->
                if (matrix[i][j] == 1) {
                    left[i][j] = (!j ? 0 : left[i][j - 1]) + 1
                }
            }
        }

        def minC = -1, maxC = -1, minR = -1, maxR = -1, ret = 0
        (0..<n).each { j ->  // 对于每一列，使用基于柱状图的方法
            int[] up = new int[m]
            int[] down = new int[m]

            Deque<Integer> stack = new LinkedList<>()   //链栈
            for (int i = 0; i < m; i++) {
                while (!stack.isEmpty() && left[stack.peek()][j] >= left[i][j]) {
                    stack.pop()
                }
                up[i] = stack.isEmpty() ? -1 : stack.peek()
                stack.push(i)
            }
            stack.clear()
            for (int i = m - 1; i >= 0; i--) {
                while (!stack.isEmpty() && left[stack.peek()][j] >= left[i][j]) {
                    stack.pop()
                }
                down[i] = stack.isEmpty() ? m : stack.peek()
                stack.push(i)
            }
            for (int i = 0; i < m; i++) {
                int height = down[i] - up[i] - 1
                int area = height * left[i][j]
                if (area > ret) {
                    ret = area
                    minC = up[i] + 1
                    maxC = down[i] - 1
                    minR = j - left[i][j] + 1
                    maxR = j
                }
            }
        }
        return [minC, maxC, minR, maxR, ret]
    }

    private static def ceil(d, scale) {
        def n = (1 / scale) as int
        (BigDecimal) Math.ceil(d * n) / n
    }

    private static def floor(d, scale) {
        def n = (1 / scale) as int
        (BigDecimal) Math.floor(d * n) / n
    }

    private static def dist(Point p1, Point p2) {
        distHarversine(p1.lat, p1.lon, p2.lat, p2.lon)
    }

    /**
     * 这里使用Harversine公式来计算两个经纬度之间的距离（单位：米）.
     *
     * @param lat1 纬度1
     * @param lon1 经度1
     * @param lat2 纬度2
     * @param lon2 经度2
     * @return 距离 单位：米
     */
    private static double distHarversine(double lat1, double lon1, double lat2, double lon2) {
        double hsinX = Math.sin(Math.toRadians(lon1 - lon2) * 0.5)
        double hsinY = Math.sin(Math.toRadians(lat1 - lat2) * 0.5)
        double h = hsinY * hsinY + (Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) * hsinX * hsinX)
        return 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)) * 6371001
    }

}