5
$\begingroup$

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)

Desired final output 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)... enter image description here

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).

Improve this question
$\endgroup$
5
  • $\begingroup$ It seems that several unrelated methods are being tried arbitrarily $\endgroup$ Commented Apr 1, 2017 at 3:08
  • $\begingroup$ Yes it is the case, I tried to explore the problem and just posted what I tried, but I didn't find a correct path to the solution $\endgroup$ Commented Apr 1, 2017 at 3:15
  • $\begingroup$ This is an interesting problem with infinitely many solutions. There are some trivial symmetries (e.g. you can't distinguish forward vs. backward), but more troublesome ones too. For example, in your plot, say a trajectory hits the clump of points in the middle, then starts going around the loop. How many times does it go around before emerging again or stopping? Adjacent points might not be part of the same pass, but subsequent passes (or even doubling back). So, some strong constraints will be necessary to solve the problem. $\endgroup$ Commented Apr 1, 2017 at 6:32
  • $\begingroup$ Right now it sounds like you have a fuzzy idea of what a solution would look like. A good first step would be to figure out how to make that mathematically precise so it can be turned into an algorithm. Thinking about the physical process that generated the data could also help. For example, if you know some bounds on the speed/acceleration of the object (or other information about how it behaves) and the time interval between sensor readings, it may give you some useful constraints. $\endgroup$ Commented Apr 1, 2017 at 6:49
  • $\begingroup$ I think the best approach here is pre-smoothing before applying isomap. A few methods are possible. I would suggest fitting a Gaussian mixture model with say 20-50 components, then using the centers of the mixture components as input to the isomap. I have no idea how to restrict the curvature though. A very interesting problem. $\endgroup$ Commented Apr 1, 2017 at 7:11

1 Answer 1

Reset to default
2
$\begingroup$

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")

Fitted curve with the described method

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.