Длина дуги кусочного сплайна с использованием R - PullRequest
Новая
       7

Длина дуги кусочного сплайна с использованием R

0 голосов
/ 26 апреля 2018

Я знаю, что есть много способов вычислить длину дуги кривой, но я ищу эффективный способ вычисления длины дуги кусочно-сплайнового через неравномерно распределенные точки.

Фактическая кривая I 'Я пытаюсь найти длину довольно сложной (контурная линия), поэтому вот быстрый пример с использованием круга, в котором фактическая длина дуги, как известно, равна 2*pi: 

# Generate "random" data
set.seed(50)
theta = seq(0, 2*pi, length.out = 50) + runif(50, -0.05, 0.05)
theta =  c(0, theta[theta >=0 & theta <= 2*pi], 2*pi)
data = data.frame(x = cos(theta), y = sin(theta))

# Bezier Curve fit
library("bezier")
bezierArcLength(data, t1=0, t2=1)$arc.length

# Calculate arc length using euclidean distance 
library("dplyr")
data$eucdist = sqrt((data$x - lag(data$x))^2 + (data$y - lag(data$y))^2)
print(paste("Euclidean distance:", sum(data$eucdist[-1])))
print(paste("Actual distance:", 2*pi))

# Output
Bezier distance: 5.864282
Euclidean distance: 6.2779
Actual distance: 6.2831

Ближайшая вещь, которую я нашелэто https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/arclength, но я должен был бы параметризовать свои данные, чтобы они были function(t) ...spline(data, t)..., чтобы использовать arclength.Я попробовал это, но подобранный сплайн проходил по середине круга, а не по окружности.

Еще одна альтернатива (безуспешно), которую я пробовал, - это подгонка кусочных сплайнов и определение длины каждого сплайна.

Любая помощь будет принята с благодарностью!

EDIT : Добавлен альтернативный метод с использованием пакета Bezier, но найденная длина дуги еще хуже, чем просто с помощью евклидова метода.

1 Ответ

0 голосов
/ 26 апреля 2018

Вместо ответов сообщества я собрал воедино решение, которое, кажется, работает для того, что я преследовал!Я оставлю свой код здесь на тот случай, если у кого-то возникнет тот же вопрос и он с этим столкнется. 

# Libraries
library("bezier")
library("pracma")
library("dplyr")
# Very slow for loops, sorry! Didn't write it as an apply function
output = data.frame()
for (i in 1:100) {
  # Generate "random" data
  # set.seed(50)
  theta = seq(0, 2*pi, length.out = 50) + runif(50, -0.1, 0.1)
  theta = sort(theta)
  theta =  c(0, theta[theta >=0 & theta <= 2*pi], 2*pi)
  data = data.frame(x = cos(theta), y = sin(theta))
  # Bezier Curve fit
  b = bezierArcLength(data, t1=0, t2=1)$arc.length
  # Pracma Piecewise cubic
  t = atan2(data$y, data$x)
  t = t + ifelse(t < 0, 2*pi, 0)
  csx <- cubicspline(t, data$x)
  csy <- cubicspline(t, data$y)
  dcsx = csx; dcsx$coefs = t(apply(csx$coefs, 1, polyder))
  dcsy = csy; dcsy$coefs = t(apply(csy$coefs, 1, polyder))
  ds <- function(t) sqrt(ppval(dcsx, t)^2 + ppval(dcsy, t)^2)
  s = integral(ds, t[1], t[length(t)])
  # Calculate arc length using euclidean distance 
  data$eucdist = sqrt((data$x - lag(data$x))^2 + (data$y - lag(data$y))^2)
  e = sum(data$eucdist[-1])
  # Use path distance as parametric variable
  data$d = c(0, cumsum(data$eucdist[-1]))
  csx <- cubicspline(data$d, data$x)
  csy <- cubicspline(data$d, data$y)
  dcsx = csx; dcsx$coefs = t(apply(csx$coefs, 1, polyder))
  dcsy = csy; dcsy$coefs = t(apply(csy$coefs, 1, polyder))
  ds <- function(t) sqrt(ppval(dcsx, t)^2 + ppval(dcsy, t)^2)
  d = integral(ds, data$d[1], data$d[nrow(data)])
  # Actual value
  a = 2*pi
  # Append to result
  output = rbind(
    output, 
    data.frame(bezier=b, cubic.spline=s, cubic.spline.error=(s-a)/a*100,
               euclidean.dist=e, euclidean.dist.error=(e-a)/a*100,
               dist.spline=d, dist.spline.error=(d-a)/a*100))
}
# Summary
apply(output, 2, mean)
# Summary output
          bezier         cubic.spline   cubic.spline.error       euclidean.dist euclidean.dist.error           dist.spline    dist.spline.error 
    5.857931e+00         6.283180e+00        -7.742975e-05         6.274913e+00        -1.316564e-01           6.283085683         -0.001585570

Я до сих пор не совсем понимаю, что делает bezierArcLength, но я очень доволен своим решением.используя cubicspline из пакета pracma, так как это намного точнее.

Другие решения по-прежнему приветствуются!

...