2D curve identification and smoothing

Question

My problem is as follows. I have 2D points representing longitude and latitude. I do not have any other information (no time-stamps). The number of points is limited, about some hundreds. The curve is not of the form Latitude = f(Longitude).

The task to achieve is to obtain a smooth version of the curve, i.e. a parametrization from t=0 to t=1 following the curve. As shown on the following example (desired output in red), I prefer to limit curvature of the predicted curve.

Are there general methods to solve this kind of problem? I already used some R packages and methods but with mixed results.

Here the example of points.

set.seed(1222)
N = 250
r = 7*runif(N)
x = r*cos(r) + rnorm(N, 0, 0.01)
y = x*r*sin(r) + rnorm(N, 0, 1)
M = data.frame(cbind(x,y))
M = M[-which(M[,1] > 4 & M[,1] < 5),]
plot(M)

For now, I tried different things:

1/ Doing a simple non-parametric smoothing do not give desired outputs.

# Take the mean over closest points
distance_M = as.matrix(dist(M))
M_new = M
eps = 5
for(i in 1:nrow(M)) {
  M_new[i,] = apply(M[which(distance_M[i,] < eps),], 2, mean)
}
plot(M)
lines(M_new, type = "p", col = "red")

2/ Manifold identification to transform points to a time-series. This works a little better, although it is very sensitive to the chosen parameter. I tested kPCA, LLE and isomap. Only isomap gives correct outputs.

## Order according to the 1st component of linear PCA
# (will not work for this example, even with a nonlinear kernel)
library(kernlab)
library(colorRamps)
kpc <- kpca(~.,data=data.frame(M), kernel="vanilladot", kpar = list(), features = 2)
out_kpc = predict(kpc,data.frame(M))[,1]
plot(M, col = blue2red(300)[cut(out_kpc,300,labels=FALSE)])

## Isomap: 
# * Global manifold learning
# * Need compute of the distance, so can be very long with > 2000 points
library(vegan)
library(colorRamps)
out_isomap = isomap(dist(M), k=5, ndim = 1)
plot(M, col=blue2red(300)[cut(out_isomap$points,300,labels=FALSE)])
# k < 5 : Data are fragmented
# 5 <= k <= 11 : OK
# k > 11 : Looks like PCA

## LLE:
# * Local manifold learning (may not be suitable for some applications)
# * Quicker than isomap
library(lle)
library(colorRamps)
out_lle = lle(M, m = 1, k=9)
plot(M, col=blue2red(300)[cut(out_lle$Y,300,labels=FALSE)])

Here is the result with isomap. From this result, I can try to smooth the obtained time series ('the obtained time series' = the curve shown in gray)...

3/ Other ideas...

• I found an article about "Geodesic Regression on Riemannian Manifolds", but I cannot find a related R package for it.

• For linear identification, orthonormal linear regression can do this task, because x and y are taken in account symmetrically. For non linear regression, I tried non linear PCA (in kernlab) and "Orthogonal Nonlinear Least-Squares Regression" (onls package) without results.

• Once the correct time series is identified, I have seen we may use "Kalman-Filter" to smooth the curve. I have shown many R examples with one output (one dimensional time series forecasting), but I cannot find a simple code example for this task (two dimensional smoothing).

Answer 1

Thank you for the comments which were helpful.

I don't have a general answer, and it could be interesting to know how to manage this problem in general (for example, how to detect roads from satellite mapping?).

For my data, I rotate the figure to obtain a relation of the form Latitude = f(Longitude). This is not very general but worked for my task. One way to automate this rotation is to do some PCA (see How to perform dimensionality reduction with PCA in R for example), however it will not work for the example plotted here.

After rotation, I can apply a regression (I lose the symmetry between Longitude and Latitude though). Finally, I get back to the original data.

Here is the code and the result:

rotate = function(x, theta) {
  exp_theta = cos(theta)+ 1i * sin(theta)
  out = ((x[1]+1i*x[2])*exp_theta)[[1]]
  return(c(Re(out), Im(out)))
}

theta = -0.15
k_gam = 15

rotate_M = t(apply(M, 1, rotate, theta = theta))
colnames(rotate_M) = c("x", "y")

## Smooth this rotated line with gam
library(mgcv)
smooth_line = gam(y ~ s(x, k = k_gam), data = data.frame(rotate_M))

## Predict for all the range
range_predict = range(rotate_M[,1])
to_predict = seq(from = range_predict[1], to = range_predict[2], length.out = 1000)
fitted_line = cbind(to_predict, predict(smooth_line, newdata = data.frame(x = to_predict)))

## Get back to the original space
out = t(apply(fitted_line, 1, rotate, theta = -theta))

## Output result
plot(M)
lines(out, col = "red")