27

I am trying to order an array of 3D coordinates by their order along a path. A sample:

points = np.array([[ 0.81127451,  0.22794118,  0.52009804],
                   [ 0.62986425,  0.4546003 ,  0.12971342],
                   [ 0.50666667,  0.41137255,  0.65215686],
                   [ 0.79526144,  0.58186275,  0.04738562],
                   [ 0.55163399,  0.49803922,  0.24117647],
                   [ 0.47385621,  0.64084967,  0.10653595]])

The points are in random order, but there is always a single path through them. I am finding the path with an adapted travelling salesman problem (TSP) soution, using the LKH solver (Helsgaun 2009). It involves two modifications:

  • Add a point at or near the origin. This finds the best starting point in every instance I've tackled so far. This was my idea, I have no other basis for it.
  • Add a point at a distance of zero from every point. This makes the solver find a route to the other end of the path. This idea was from this SO question.

Note that the TSP does not involve position, only the distances between nodes. So the solver does 'know' (or care) that I'm working in 3D. I just make a distance matrix like so:

import numpy as np
from scipy.spatial.distance import pdist, squareform

# Add a point near the origin.
points = np.vstack([[[0.25, 0, 0.5]], points])
dists = squareform(pdist(points, 'euclidean'))

# Normalize to int16s because the solver likes it.
d = 32767 * dists / np.sqrt(3)
d = d.astype(np.int16)

# Add a point that is zero units from every other point.
row, col = d.shape
d = np.insert(d, row, 0, axis=0)
d = np.insert(d, col, 0, axis=1)

I pass this to my fork of pytsp, which passes it to the LKH solver. And everything is fine... except when the path crosses itself. TSP solutions cannot have closed loops, so I always get the open loop shown on the right here:

Travelling salesman paths

Note that this is an analogous 2D version of my situation. Note also that the points are imperfectly aligned, even along the 'straight' bits.

So my question is: how can I help the solver preserve the direction of the path whenever possible? I've got two ill-formed ideas, but so far been unable to implement anything:

  • Use another metric instead of L2. But I don't think this can work, because at a given junction, there's nothing inherently different about the 'wrong' point. Its wrongness depends on the previous point. And we don't know yet which is the previous point (that's what we're trying to figure out). So I think this is no good.
  • Evaluate the local colinearity of every set of three points (e.g. using the determinant of every triple). Modulate the local '3D slope' (not sure what I mean) by this colinearity coefficient. Give every point another dimension expressing this local alignment. Now the norm will reflect local alignment and (hopefully) roughly colinear things will join up.

I have put these files on Dropbox:

Thank you for reading; any ideas appreciated.

Reference

K. Helsgaun, General k-opt submoves for the Lin-Kernighan TSP heuristic. Mathematical Programming Computation, 2009, doi: 10.1007/s12532-009-0004-6.

Community
  • 1
  • 1
Matt Hall
  • 7,614
  • 1
  • 23
  • 36
  • 3
    If you have one particular problem that you are trying to solve, and you can afford to do some iterative customized "intervention", you could run your code once, identify which points are supposed to be neighbors in the path, manually override their distances in the `dists` array to make them appear closer, rerun the algorithm, identify any additional points that should be neighbors, etc. Of course, if you will be repeatedly solving different instances of the problem and you need a solver that "just works" with no additional intervention, this idea is no help. – Warren Weckesser Aug 29 '16 at 04:21
  • @WarrenWeckesser Thank you for the thought. I tried using a BallTree to get local density, which does highlight these problem areas. In theory I could try rewiring things in dense spots. But yes — I'm still hoping there's a trick I can use to sidestep that kind of post hoc fiddling. For now, it's my fall-back. Cheers! – Matt Hall Aug 29 '16 at 11:07
  • 3
    This doesn't seem like a travelling salesman problem. You have an implicit objective function which is something different from the length of the path. Perhaps you could start by making your objective function explicit. Something like total length + sum of all magnitudes of direction changes. A direction change requires 3 successive points to compute (which is why I said that it doesn't seem like a TSP). Once you get your objective function clarified, you could use a general-purpose heuristic such as an evolutionary algorithm. – John Coleman Aug 29 '16 at 12:25
  • @JohnColeman Thank you for the input. I agree it feels like it's a step too far from the TSP to continue down that road. I've never solved a problem in the way you're describing before; if you can point me at a simple example, I'd appreciate it. Meanwhile I'll start up the Google. Onwards + upwards! – Matt Hall Aug 29 '16 at 13:00
  • Your formula for choosing a point needs to incorporate a weighted direction as well as distance so straight gives a lower score than turned. Think of the matrix as being a matrix of functions. The level of computation complexity goes up rapidly by a power of the size of the matrix. A two formula linear problem with the intersection being optimized to the lowest value of the intersection of these two curves. As distance vs direction is weighted in some way in your formula, say by a constant, to keep it from turning, at some distance it will turn which I think you would want. Hope this helps. – Tim Cederquist Aug 30 '16 at 21:13
  • No time for a full answer, but there seems to be rich literature on the topic by searching for 'curve reconstruction from point cloud'. I even found one that [is based on TSP](http://epubs.siam.org/doi/pdf/10.1137/S0097539700366115). – ojdo Aug 31 '16 at 14:32
  • Thank you @TimCederquist and ojdo for the remarks. As you both suggest, it seems there are ways to penalize turns. I need to do some more research. – Matt Hall Sep 01 '16 at 15:34

2 Answers2

1

Judging by the documentation on pytsp, the distance matrix doesn't have to be symmetric. This means that you could modify the L2 norm to incorporate information on a preferred direction into that matrix. Say you have a preferred direction for some pairs of points (i,j), then for each of these point you could divide dists[i,j] by (1+a) and multiply dists[j,i] by (1+a) to make that direction more favourable. This means that if your algorithm is sure to find the global optimum, you can force it to satisfy your preferred direction by taking a is sufficiently large.

Also, I'm not sure it's impossible to have closed loops in a solution where the distance matrix is taken from 3D data. It seems to me that the 'no closed loops' is a result (of the triangle inequality) specific to 2D.

KeithWM
  • 1,295
  • 10
  • 19
  • Thanks for this answer. I think you are spot on with that observation about the triangle inequality in 3D — in fact, the path does not perfectly cross. As for the asymmetric distance matrix... I had not thought about that either. I need to try it to see if it will help. The tricky thing is that I don't know in advance in which direction I want to traverse a given portion of the graph. Cheers! – Matt Hall Sep 01 '16 at 15:32
0

Perhaps if you were to include a mechanism for directional momentum that preferred nodes that fall along the existing trajectory - that may make the cost function optimize for what you're looking to achieve. I think this could be well-phrased as a local search problem and solved with a metaheuristic such as simulated annealing.

The advantage of solving the problem in this way is that you can consider the state of the entire solution when calculating the distance between two nodes.

Ryan Codrai
  • 172
  • 8