Мне нужно кодировать алгоритм сканирования Грэма, чтобы найти выпуклую оболочку заданного набора точек. Я много пробовал и не могу найти решение.
У меня есть следующий код для функции Грэма
РЕДАКТИРОВАТЬ: добавить комментарии, чтобы помочь понять поток кода
grahamScanConvexHull<-function(points){
# Points is a dataset with eeach row a point one column x and one column
# (We can assume that there are not three points collinear)
# Vertices is a simulated stack were the solution is stored
vertices <- data.frame(x=double(), y=double())
flag <- T # Flag used in case there two points with the least y
# Get the most left down_point and leave others in points
points <- quicksort_bottom_top(points, 1, nrow(points))
# if there are two points with the same lowest y we get the leftmost
if(points[1,]$y == points[2,]$y){
if(points[1,]$x < points[2,]$x){
point_0 <- points[1,]
points <- points[(2:nrow(points)),]
}else{
point_0 <- points[2,]
points <- points[c(1, 3:nrow(points)),]
}
}else{
point_0 <- points[1,]
points <- points[(2:nrow(points)),]
flag <- F
}
# Sort the points related to de most left point
points <- quicksort_angle(points, 1, nrow(points), point_0)
# This is used in case that there are two points with the lowest y be put first
# in the list of points (it will we the last in other case)
if(flag){
first_point <- points[nrow(points),]
points <- rbind(first_point, points[(1:(nrow(points)-1)),])
}
# Push the point with the lowest y (ledtmost) and the point with lest polar
#angle
vertices <- push(vertices, point_0)
vertices <- push(vertices, points[1,])
for(i in 1:nrow(points)){
# Get the top 2 points in the stack as vectors
next_to_top <- c(next_to_top(vertices)$x,next_to_top(vertices)$y)
top <- c(top(vertices)$x,top(vertices)$y)
# Get the oncheck point as a vector
point_i <- c(points[i,]$x,points[i,]$y)
while(nrow(vertices) >= 2 & (rightTurn(next_to_top, top, point_i) <= 0)){
vertices <- pop(vertices)
# Recalculate the top two points of the stack (as vector) after pop
next_to_top <- next_to_top(vertices)
next_to_top <- c(next_to_top$x, next_to_top$y)
top <- top(vertices)
top <- c(top$x, top$y)
}
# Push the onCheck point into the stack
vertices <- push(vertices, points[i,])
}
return(vertices)
}
Для управления стеком у меня есть код этой короткой вспомогательной функции
top <- function(stack){
return(stack[1,])
}
next_to_top <- function(stack){
return(stack[2,])
}
pop <- function(stack){
return(stack[-1,])
}
push <- function(stack, e){
if(e$x %in% stack$x & e$y %in% stack$y ){
if(which(e$x %in% stack$x) == which(e$y %in% stack$y)){
return(stack)
}else{
stack <- rbind(e, stack)
return(stack)
}
}else{
stack <- rbind(e, stack)
return(stack)
}
}
А для против часовой стрелки у меня есть следующая функция:
rightTurn <- function(A, B, C) {
val <- (B[1] - A[1]) * (C[2] - A[2]) -
(B[2] - A[2]) * (C[1] - A[1])
return(val)
}
Функции сортировки, которые я проверял много раз и, я вполне уверен, что они в порядке (также я использую их в течение длительного времени без каких-либо проблем). Я не знаю, что я делаю неправильно, поэтому любая помощь приветствуется.
РЕДАКТИРОВАТЬ:
Код quicksort_bottom_top ()
quicksort_bottom_top <- function(points, low, high) {
if (low < high) {
list_particion <- partition_bottom_top(points, low, high)
p <- list_particion[['p']]
points <- list_particion[['points']]
points <- quicksort_bottom_top(points, low, p - 1)
points <- quicksort_bottom_top(points, p + 1, high)
}
return(points)
}
partition_bottom_top <- function(points, low, high) {
pivot <- points[high, ]
i <- low
for (j in low:high) {
if (points[j, ]$y < pivot$y) {
points <- swap(points, i, j)
i <- i + 1
}
}
points <- swap(points, i, high)
return(list('points' = points, 'p' = i))
}
Код quicksort_angle ()
quicksort_angle<-function(points, low, high, reference_point){
if (low < high) {
list_particion <- partition_angle(points, low, high, reference_point)
p <- list_particion[['p']]
points <- list_particion[['points']]
points <- quicksort_angle(points, low, p - 1, reference_point)
points <- quicksort_angle(points, p + 1, high, reference_point)
}
return(points)
}
partition_angle<-function(points, low, high, reference_point){
pivot <- points[high,]
i <- low
for(j in low:high){
if(less_angle(points[j,], pivot, reference_point)){
points<-swap(points, i, j)
i<-i+1
}
}
points<-swap(points, i, high)
return(list('points'=points, 'p'=i))
}
less_angle<-function(check_point, pivot_point, reference_point){
# m_check = (check_point$y - reference_point$y) / (check_point$x - reference_point$x)
# m_pivot = (pivot_point$y - reference_point$y) / (pivot_point$x - reference_point$x)
m_check = (reference_point$y - check_point$y) / (reference_point$x - check_point$x)
m_pivot = (reference_point$y - pivot_point$y) / (reference_point$x - pivot_point$x)
if((m_pivot * m_check) < 0){
return(m_check > m_pivot)
}else{
return(m_check < m_pivot)
}
}
И пример использования (красный - это контрольная точка, а синий - отсортированные точки) 
И вспомогательная функция swap ()
swap <- function(points, i, j) {
temp <- points[i, ]
points[i, ] <- points[j, ]
points[j, ] <- temp
return(points)
}
РЕДАКТИРОВАТЬ: пример данных (dput(head(points, n = 10)))
structure(list(x = c(407.377968561836, 320.925913405139, 806.236479437212,
682.732644010102, 59.5188127646688, 195.601556062, 290.540224567521,
449.757491987199, 760.845693320502, 544.514419493731), y = c(255.027047196869,
63.905220010085, 922.746908532688, 956.717615367146, 244.059993790463,
498.689062974649, 986.184075321769, 605.362085530069, 254.609465040965,
711.789759260369)), row.names = c(NA, 10L), class = "data.frame")
РЕДАКТИРОВАТЬ: решения
Мое решение (vertices <- grahamScanConvexHull(points); dput(head(vertices, n = 10)))
structure(list(x = c(59.5188127646688, 320.925913405139), y = c(244.059993790463,
63.905220010085)), row.names = c(10L, 1L), class = "data.frame")
Функция библиотеки (vertices <- points[chull(points),];dput(head(vertices, n = 10)))
structure(list(x = c(760.845693320502, 320.925913405139, 59.5188127646688,
290.540224567521, 682.732644010102, 806.236479437212), y = c(254.609465040965,
63.905220010085, 244.059993790463, 986.184075321769, 956.717615367146,
922.746908532688)), row.names = c(9L, 2L, 5L, 7L, 4L, 3L), class = "data.frame")