학교 프로젝트로 시각장애인분들을 위한 횡단보도 보행 보조 어플리케이션을 제작하고 있다.

플랫폼은 안드로이드로 제작하고 있고, 지도 API는 네이버 지도를 이용한다.
횡단보도 데이터는 서울특별시 교통안전시설물관리시스템 T-GIS에서 네트워크 요청을 통해 가져올 수 있도록 했다.

TGIS에서 받아온 횡단보도 데이터를 초록색 폴리곤으로 나타내면 이렇게 보인다.

횡단보도 안내를 위해서 방향에 대한 정보도 필요한데, 제공되는 횡단보도 데이터에는 방향정보는 없다.
그래서 횡단보도의 방향을 정해야하는데, 위에 보시다시피 횡단보도 다각형이 변들도 많고 들쭉날쭉해서 복잡하다.

횡단보도가 사각형으로 만들어주면 이런식으로 방향을 추정해볼 수 있을 것 같다.
그래서 횡단보도를 사각형으로 만들기 위해, 도형을 단순화시킬 수 있는 알고리즘을 찾았다.

왼쪽부터 기본, Convex hull, Ramer-Douglas-Peucker, OpenCV minAreaRect, minimumOrientedBoundingBox
minimumOrientedBoundingBox가 가장 정확하다.

0. 네이버 지도에 폴리곤 그리기

val list = listOf(
            LatLng(37.56409886267795, 126.97816863257054),
            LatLng(37.56409886267800, 126.97816863257054),
            LatLng(37.56409886267795, 126.97816863257154)
        )
val polygon = PolygonOverlay().apply {
        coords = list; color = Color.WHITE
        outlineWidth = 5; outlineColor = Color.GREEN
    }
polygon.map = naverMap

폴리곤 오버레이를 생성하고, 옵션들을 원하는 대로 설정하고 map 속성에 네이버 지도 객체를 지정해주면 된다.
폴리곤 오버레이를 생성할 때, 위경도 좌표 리스트를 인자로 주면 위치를 설정할 수 있다.
단, 리스트에 좌표가 3개 미만이라면 오류가 발생하므로 주의해야 한다.

1. Convex Hull 알고리즘

Convex Hull은 어떤 다각형의 내부점들을 모두 포함하는 볼록다각형으로 만든다.

아래 같은 횡단보도라면


이렇게 될 것이다.

Convex Hull 의 순서

1. 좌표점들 중 가장 왼쪽 아래의 점을 찾는다.
2. 다른 점들의 리스트를 찾은 점 기준 시계 반대 방향으로 정렬한다. 두 점이 같은 각도인 경우, 가까운 것을 기준으로 정렬한다.
3. 기준점을 스택에 넣는다. 그 다음 정렬된 점들을 스택에 하나씩 넣으면서 Convex Hull인 점들을 찾는다.
   3-1. Convex Hull은 볼록다각형이니까, 점들의 각이 좌회전만 해야한다. 우회전하면 오목이니까.
   3-2. 그래서 정렬된 점들 중 넣을 차례인 점이 스택에 있는 점과 좌회전을 이루면 스택에 push한다. 우회전이면 스택 맨 위의 점을 pop하여 이전 점을 제거한다. 제거 후 남은 점들과 좌회전을 이루거나 스택에 점이 하나만 남으면 스택에 push해주고, 제거 후에도 우회전이면 pop하는 과정을 반복한다.
4. 완성

CCW 알고리즘

위에서 말한 시계 방향, 시계 반대 방향, 우회전, 좌회전은 CCW (Counter Clock Wise) 알고리즘을 이용해서 판별할 수 있다.
CCW 알고리즘은 점 3 개에 대하여 신발끈 공식을 이용했을 때, 양수이면 시계반대방향(우회전), 음수면 시계방향(좌회전), 0이면 직선 상에 있는 것이다.

구현

네이버 지도의 PolygonOverlay에 확장함수로 만들어서 사용할 수 있게했다.
typealias로 Pair<BigDecimal,BigDecimal>을 Point로 써주었다.
PolygonOverlay의 좌표들은 double 형태인데, 그대로 계산하면 double의 정확도 때문에 오차가 발생한다.
그러므로 BigDecimal을 이용해서 계산해준다.

typealias Point = Pair<BigDecimal, BigDecimal>

fun PolygonOverlay.getConvexHull(): PolygonOverlay {

    val arr = this.coords.map { Point(it.latitude.toBigDecimal(), it.longitude.toBigDecimal()) }
        .toMutableList()//좌표 리스트
    //1. 가장 왼쪽 아래 좌표 찾기
    var bottomPoint = Point(90.0.toBigDecimal(), 180.0.toBigDecimal())
    for (point in arr) {
        if (point.second < bottomPoint.second) {
            bottomPoint = point
        } else if (point.second == bottomPoint.second) {
            if (point.first < bottomPoint.first) {
                bottomPoint = point
            }
        }
    }
    //2. ccw 정렬
    arr.sortWith(kotlin.Comparator { o1, o2 ->
        val result = ccw(bottomPoint, o1, o2)
        if (result > 0) return@Comparator -1
        else if (result < 0) return@Comparator 1
        else {
            if (dist(bottomPoint, o1) > dist(bottomPoint, o2)) {
                return@Comparator 1
            }
        }
        return@Comparator -1
    })
	//3. convex hull 찾기
    val stack = Stack<Point>()
    stack.add(bottomPoint)
    for (i in 1 until arr.size) {
        while (stack.size > 1 && ccw(stack[stack.size - 2], stack[stack.size - 1], arr[i]) <= 0) {
            stack.pop()
        }
        stack.add(arr[i])
    }
    return PolygonOverlay(stack.toList().map { LatLng(it.first.toDouble(), it.second.toDouble()) })
}

fun ccw(p1: Point, p2: Point, p3: Point): Int {
    /*
    x1 x2 x3 x1
    y1 y2 y3 y1
     */
    val result =
        (p1.first * p2.second + p2.first * p3.second + p3.first * p1.second) - (p2.first * p1.second + p3.first * p2.second + p1.first * p3.second)
    if (result > 0.toBigDecimal()) return 1
    if (result < 0.toBigDecimal()) return -1
    return 0
}

fun dist(p1: Point, p2: Point): BigDecimal {
    return (p2.first - p1.first) * (p2.first - p1.first) + (p2.second - p1.second) * (p2.second - p1.second)
}

결과

왼쪽이 원본, 오른쪽이 convex hull

오목했던 다각형들이 대부분 사다리꼴 형태가 되었다.
여기까지만해도 각 변에 중점을 이용해 방향을 파악할 수 있을 것 같다.

그런데 가끔 육각형인 애들도 있다.
사각형으로 만들려면 어떻게 할까 하다가 Ramer-Douglas-Peucker 알고리즘을 찾았다.

2. Ramer-Douglas-Peucker 알고리즘

Ramer-Douglas-Peucker 알고리즘은 다각형을 단순하게 변환할 수 있는 알고리즘이다.

Ramer-Douglas-Peucker 알고리즘 순서

0. 한계치를 설정한다.
1. 먼저 좌표점들 중 서로 사이가 가장 먼 두 정점을 찾는다.
2. 1에서 찾은 가장 먼 두 좌표점이 이루는 선분(그림의 빨간선)과 나머지 정점들의 거리를 계산하여, 가장 먼 정점을 찾는다.
3. 2에서 찾은 가장 먼 정점의 거리(그림의 노란선)가 0번의 한계치 이상이라면, 해당 정점을 기준으로 선분을 둘로 나누고 또 다시 선분에서 가장 먼 정점을 찾는 작업을 재귀적으로 수행한다.
4. 한계치보다 작으면, 선분의 양끝을 꼭짓점에 추가하고 종료한다.

Rotating Calipers 알고리즘

위의 1번 과정에서 좌표점들 중 서로 사이가 가장 먼 두 정점은 어떻게 찾을까.
두 점을 잡고 일일이 계산하면 n^2 시간이 필요하지만, Rotating Calipers라는 알고리즘을 사용하면 n시간 안에 가능하다.
자세한 내용 (외부 링크)

코드

PolygonOverlay의 확장함수로 getSimplified라고 만들었다.
인자로 ep 한계치를 설정할 수 있다.

코드들은 구글링을 통해 구하고, BigDecimal 사용을 위해 적당히 수정해주었다.
또, BigDecimal의 루트 계산을 위해 Big-Math 라이브러리를 이용했다.

코드 출처

Rotating Caplipers 코드
Ramer-Douglas-Peucker 코드

//다각형 단순화
fun PolygonOverlay.getSimplified(ep: Double) = PolygonOverlay(RamerDouglasPeucker.douglasPeucker(
    this.rotatingCalipers().map { arrayOf(it.first, it.second) }, ep.toBigDecimal()
).map { LatLng(it[0].toDouble(), it[1].toDouble()) })

//rotating calipers를 이용해서 가장 먼 두점을 찾고, 정점 리스트 앞뒤에 배치하여 반환
fun PolygonOverlay.rotatingCalipers() : List<Point> {
    val convexHull = this.getConvexHull().coords.map { Point(it.latitude.toBigDecimal(),it.longitude.toBigDecimal()) }.toMutableList()
    var max_dist:BigDecimal = BigDecimal.ZERO
    var aIndex = 0
    var bIndex = 0
    var j = 1
    for (i in convexHull.indices) {
        val i_next: Int = (i + 1) % convexHull.size
        while (true) {
            val j_next: Int = (j + 1) % convexHull.size
            val bx = convexHull[i_next].first - convexHull[i].first // 왼쪽 벡터
            val by = convexHull[i_next].second - convexHull[i].second
            val cx = convexHull[j_next].first - convexHull[j].first // 오른쪽 벡터
            val cy = convexHull[j_next].second - convexHull[j].second
            val ccw = ccw(Point(BigDecimal.ZERO, BigDecimal.ZERO), Point(bx, by), Point(cx, cy))
            j = if (ccw > 0) {  // 반시계 방향이면 오른쪽에 있는 점을 다음으로
                j_next
            } else {    // 시계 방향이면 왼쪽에 있는 점을 다음으로
                break
            }
        }

        // 최대 거리 구하기
        if (dist(convexHull[i], convexHull[j]) > max_dist) {
            max_dist = dist(convexHull[i], convexHull[j])
            aIndex = i
            bIndex = j
        }
    }
    val temp = this.coords.map { Point(it.latitude.toBigDecimal(),it.longitude.toBigDecimal()) }.toMutableList()
    temp.remove(convexHull[aIndex])
    temp.remove(convexHull[bIndex])
    return (listOf(convexHull[aIndex]) + temp + listOf(convexHull[bIndex]))
}

object RamerDouglasPeucker {
    private fun sqr(x: BigDecimal): BigDecimal {
        return x * x
    }

    private fun distanceBetweenPoints(vx: BigDecimal, vy: BigDecimal, wx: BigDecimal, wy: BigDecimal): BigDecimal {
        return sqr(vx - wx) + sqr(vy - wy)
    }

    private fun distanceToSegmentSquared(
        px: BigDecimal,
        py: BigDecimal,
        vx: BigDecimal,
        vy: BigDecimal,
        wx: BigDecimal,
        wy: BigDecimal
    ): BigDecimal {
        val l2 = distanceBetweenPoints(vx, vy, wx, wy)
        if (l2 == BigDecimal.ZERO) return distanceBetweenPoints(px, py, vx, vy)
        val t = ((px - vx) * (wx - vx) + (py - vy) * (wy - vy)) / l2
        if (t < BigDecimal.ZERO) return distanceBetweenPoints(px, py, vx, vy)
        return if (t > BigDecimal.ONE) distanceBetweenPoints(px, py, wx, wy) else distanceBetweenPoints(
            px, py,
            vx + t * (wx - vx),
            vy + t * (wy - vy)
        )
    }

    private fun perpendicularDistance(
        px: BigDecimal,
        py: BigDecimal,
        vx: BigDecimal,
        vy: BigDecimal,
        wx: BigDecimal,
        wy: BigDecimal
    ): BigDecimal {
        return DefaultBigDecimalMath.sqrt(distanceToSegmentSquared(px, py, vx, vy, wx, wy))
    }

    private fun douglasPeucker(
        list: List<Array<BigDecimal>>,
        s: Int,
        e: Int,
        epsilon: BigDecimal,
        resultList: MutableList<Array<BigDecimal>>
    ) {
        // Find the point with the maximum distance
        var dmax = BigDecimal.ZERO
        var index = 0
        val end = e - 1
        for (i in s + 1 until end) {
            // Point
            val px = list[i][0]
            val py = list[i][1]
            // Start
            val vx = list[s][0]
            val vy = list[s][1]
            // End
            val wx = list[end][0]
            val wy = list[end][1]
            val d = perpendicularDistance(px, py, vx, vy, wx, wy)
            if (d > dmax) {
                index = i
                dmax = d
            }
        }
        // If max distance is greater than epsilon, recursively simplify
        if (dmax > epsilon) {
            // Recursive call
            douglasPeucker(list, s, index, epsilon, resultList)
            douglasPeucker(list, index, e, epsilon, resultList)
        } else {
            if (end - s > 0) {
                resultList.add(list[s])
                resultList.add(list[end])
            } else {
                resultList.add(list[s])
            }
        }
    }

    fun douglasPeucker(list: List<Array<BigDecimal>>, epsilon: BigDecimal): List<Array<BigDecimal>> {
        val resultList: MutableList<Array<BigDecimal>> = ArrayList()
        douglasPeucker(list, 0, list.size, epsilon, resultList)
        return resultList
    }
}

결과

무엇인가 잘못됐다. 어디가 잘못된걸까..
Ramer-Douglas-Peucker의 문제는 사각형을 보장하지도 않는다는 점과 설정해주는 한계치에 따라 원치 않는 결과를 낸다는 것이다. 삼각형이 나오거나, 그냥 직선하나로 끝나버리기도 한다.

3. OpenCV minAreaRect

Ramer-Douglas-Peucker 사용에 실패하고, 다른 방법을 찾던 중 OpenCV 라는 것을 찾았다.
컴퓨터비젼 라이브러리인데, 딱 내가 원하던 기능이 있었다.
바로 사용해봤다.!

Android OpenCV 설치 방법

Android OpenCV 설치하기

1. sdk 설치

Android Studio에서 쉬프트 두번 눌러서 sdk manager 검색, sdk tools 탭을 연다.
NDK와 CMake의 설치를 확인한다.

2. OpenCV sdk 다운로드

OpenCV Android SDK 4.5.5 링크
링크에서 파일을 다운받고 압축을 푼다.

3. 모듈 임포트

Import Module 메뉴를 열고

압축 푼 경로를 찾아준다. 그 다음 완료
그리고 Project Structure 메뉴로 간다.

저기 +버튼 누르고 opencv 모듈을 dependency로 추가해준다.

4. 모듈 로드하기

override fun onCreate() {
      super.onCreate()
      val isIntialized = OpenCVLoader.initDebug()
  }

어플리케이션이나 액티비티에서 모듈을 로드해준다.

코드

fun PolygonOverlay.minAreaRect(): PolygonOverlay {
    val points = this.coords.map { org.opencv.core.Point(it.latitude, it.longitude) }
    val mat = MatOfPoint2f()
    mat.fromList(points)
    val rectangle = Imgproc.minAreaRect(mat)
    val rectanglePoints = Array<org.opencv.core.Point?>(4) { null }
    rectangle.points(rectanglePoints)
    return PolygonOverlay(rectanglePoints.toList().map { LatLng(it!!.x, it.y) })
}

결과

왼쪽이 원본(초록색), 오른쪽이 처리된 데이터(빨간색)이다.
거의 맞는거 같다.
그런데 약간씩 어긋나 있다.
역시 double의 정확도..

4. minimumOrientedBoundingBox()

이걸 정확하게 찾고 싶다.

구글링을 해보니, rotateCalipers처럼 도형을 돌려서 찾을 수 있다고 한다.

순서

1. Convex Hull 을 찾는다. (위 그림의 빨간색)
2. Convex Hull의 각 정점의 각도를 찾는다.
3. Convex Hull을 회전시킨다.
4. 회전시킨 Convex Hull을 감싸는 사각형의 넓이를 구한다. (노란색 박스)
5. 2-4를 모든 점에 대해 반복하면서, 4의 사각형이 가장 좁은 정점을 기억한다.
6. 5에서 가장 좁은 사각형을 갖게 하는 정점을 찾았다면, 해당 정점에서 회전시킨 4번 사각형을 구하고 다시 원래 각도만큼 반대로 회전시키면 완성

코드

역시 확장함수로 만들었고, BigDecimal을 사용하였다.
Big-Math 라이브러리 필요

fun PolygonOverlay.orientedMinimumBoundingBox(): PolygonOverlay {
    val convexHull = this.getConvexHull().coords.map {
        Point(
            it.latitude.toBigDecimal(),
            it.longitude.toBigDecimal()
        )
    }.toTypedArray()
    val alignPointIndex = computeAlignPointIndex(convexHull)
    val r = computeAlignedPolygon(convexHull,alignPointIndex)
    val pob = computePointOfBound(r)
    val angle = computeEdgeAngleRad(convexHull,alignPointIndex)
    val omb = rotatePolygon(pob,angle,convexHull[alignPointIndex])
    return PolygonOverlay(omb.map { LatLng(it.first.toDouble(),it.second.toDouble()) })
}

fun computeAlignPointIndex(pg: Array<Point>): Int {
    var minArea = BigDecimal.TEN
    var minAreaIndex = -1
    for (i in pg.indices) {
        val r = computeAlignedPolygon(pg, i)
        val pob = computePointOfBound(r)
        val area = computeAreaOfPointOfBound(pob)
        if (area < minArea) {
            minArea = area
            minAreaIndex = i
        }
    }
    return minAreaIndex
}

fun computeAreaOfPointOfBound(pg:Array<Point>): BigDecimal {
    return (pg[3].first - pg[0].first) * (pg[3].second - pg[0].second)
}

fun computePointOfBound(pg: Array<Point>): Array<Point> {
    var minX = (90.0).toBigDecimal()
    var minY = (180.0).toBigDecimal()
    var maxX = (-90.0).toBigDecimal()
    var maxY = (-180.0).toBigDecimal()
    for (p in pg) {
        val x = p.first
        val y = p.second
        minX = if (minX > x) x else minX
        minY = if (minY > y) y else minY
        maxX = if (maxX < x) x else maxX
        maxY = if (maxY < y) y else maxY
    }
    return arrayOf(
        Point(minX, minY), Point(minX, maxY),
        Point(maxX, minY), Point(maxX, maxY)
    )
}

fun computeAlignedPolygon(pg: Array<Point>, i: Int): Array<Point> {
    val angle = computeEdgeAngleRad(pg, i)
    return rotatePolygon(pg, -angle, i)
}

fun computeEdgeAngleRad(
    points: Array<Point>, index: Int
): BigDecimal {
    val i1 = (index + 1) % points.size
    val p0 = points[index]
    val p1 = points[i1]
    val dx = p1.first - p0.first
    val dy = p1.second - p0.second
    return DefaultBigDecimalMath.atan2(dy, dx)
}

fun rotatePolygon(
    pg: Array<Point>,
    rotAngle: BigDecimal,
    i: Int
): Array<Point> {
    var x: BigDecimal
    var y: BigDecimal
    return Array(pg.size) {
        x = pg[it].first - pg[i].first
        y = pg[it].second - pg[i].second
        val nx =
            pg[i].first + (x * DefaultBigDecimalMath.cos(rotAngle) - y * DefaultBigDecimalMath.sin(
                rotAngle
            ))
        val ny =
            pg[i].second + (x * DefaultBigDecimalMath.sin(rotAngle) + y * DefaultBigDecimalMath.cos(
                rotAngle
            ))
        Point(nx, ny)
    }
}
fun rotatePolygon(
    pg: Array<Point>,
    rotAngle: BigDecimal,
    centroid: Point
): Array<Point> {
    var x: BigDecimal
    var y: BigDecimal
    return Array(pg.size) {
        x = pg[it].first - centroid.first
        y = pg[it].second - centroid.second
        val nx =
            centroid.first + (x * DefaultBigDecimalMath.cos(rotAngle) - y * DefaultBigDecimalMath.sin(
                rotAngle
            ))
        val ny =
            centroid.second + (x * DefaultBigDecimalMath.sin(rotAngle) + y * DefaultBigDecimalMath.cos(
                rotAngle
            ))
        Point(nx, ny)
    }
}


fun polygonCenterOfMass(pg: Array<Point>): Point {

    val N: Int = pg.size
    val polygon = Array<Point>(N) {
        Point(pg[it].first, pg[it].second)
    }
    var cx = BigDecimal.ZERO
    var cy = BigDecimal.ZERO
    val A = PolygonArea(polygon, N)
    var i: Int
    var j: Int
    var factor = BigDecimal.ZERO
    i = 0
    while (i < N) {
        j = (i + 1) % N
        factor = polygon[i].first * polygon[j].second - polygon[j].first * polygon[i].second
        cx += (polygon[i].first + polygon[j].first) * factor
        cy += (polygon[i].second + polygon[j].second) * factor
        i++
    }
    factor = 1.0.toBigDecimal() / (6.0.toBigDecimal() * A)
    cx *= factor
    cy *= factor
    return Point(
        cx.abs(),
        cy.abs()
    )
}


fun PolygonArea(polygon: Array<Point>, N: Int): BigDecimal {
    var i: Int
    var j: Int
    var area = BigDecimal.ZERO
    i = 0
    while (i < N) {
        j = (i + 1) % N
        area += polygon[i].first * polygon[j].second
        area -= polygon[i].second * polygon[j].first
        i++
    }
    area /= 2.0.toBigDecimal()
    return area.abs()
}

코드 출처 및 참고

https://stackoverflow.com/questions/23905127/oriented-minimum-bounding-box
http://www.java2s.com/example/java-utility-method/polygon-rotate/rotatepolygon-polygon-pg-double-rotangle-point-centroid-polygon-original-ee480.html

결과

가장 잘된다.
폴리곤은 이상하게 나오지만, 정점들이 사각형을 이루고 방향이 정확하다.

근데 좀 느리다. 좀 더 최적화할 수 있는지 생각해봐야한다.