What is an efficient Algorithm to draw a straight line from points?

Question

I have some cross section survey data, which is like this:

X Y Z pt_dt
2970969.635 359725.0088 83.4242 1-x
2970968.278 359722.2182 83.2591 1-x
2970941.771 359670.127  83.0655 1-x
2970961.369 359708.6424 83.4785 1-x
2970950.888 359688.048  85.2994 1-x
2970955.084 359696.2922 83.6479 1-x
2970953.47  359693.1165 84.9628 1-x
2970948.813 359683.977  83.1451 1-x
2970931.904 359650.7413 83.1042 1-x

When I use the X & Y to make points, I can see the points on the map like this: enter image description here

What I need to do now, is to get the points in linear order (i.e. from left to right, or right to left), so that I can create a 3D line, which can serve as the Cross sectional profile.

My current algorithm, is to calculate the distance between every pair of points(For one given cross section, which is indicated by pt_dt). The pair which has the maximum distance, are the outermost points. I then arrange the points based on distance from one of the outermost points. I have implemented this algorithm in python, and it works.

This algorithm is O(n2); My entire data consists of about 5000 such surveys, so any complexity that I can reduce will be beneficial.

So Is there a more efficient way, or algorithm to do this?

If your points are on an almost straight line, max/min X and Y should give you the outermost points without distance calculation of all pairs. — AndreJ, Aug 27 '14 at 06:04
@Andre A lexicographic sort would do the same but be more robust and reliable. Once the endpoints (X0,Y0) and (X1,Y1) are found, all points (X,Y) can then be sorted by their dot products with the direction vector; that is, by the values (X-X0)(X1-X0) + (Y-Y0)(Y1-Y0). — whuber, Aug 27 '14 at 14:40
@whuber: Can you explain what 'lexicographic sort' would mean in this context? I've only ever used this sorting for strings. — Devdatta Tengshe, Aug 28 '14 at 02:47
@Devdatta Sorry for being so terse: I was referring to a lexicographic sort on the x and y coordinates. (Such sorting is often a first step in many algorithms of computational geometry.) It is accomplished simply by sorting on one of the coordinates and then using the second coordinate to resolve any ties. This sort is built in to many programming environments. Finding the endpoints needs only O(n*log(n)) time and O(n) space (or O(n) time with a clever algorithm). It works well here because it automatically handles points that are lined up perfectly vertically or perfectly horizontally. — whuber, Aug 28 '14 at 14:16

score 4 · Accepted Answer · edited Apr 13 '17 at 12:33

Taking into account Andre's comment, and for the sake of learning I modified another user's answer to create a 3D Line from your sample data. I hope you or someone else might find it useful:

enter image description here

X,Y,Z,pt_dt
2970969.635,359725.0088,83.4242,1-x
2970968.278,359722.2182,83.2591,1-x
2970941.771,359670.127,83.0655,1-x
2970961.369,359708.6424,83.4785,1-x
2970950.888,359688.048,85.2994,1-x
2970955.084,359696.2922,83.6479,1-x
2970953.47,359693.1165,84.9628,1-x
2970948.813,359683.977,83.1451,1-x
2970931.904,359650.7413,83.1042,1-x

import csv
import numpy as np
try:
    from osgeo import ogr,osr
except ImportError:
    import ogr,osr

SHP_FILENAME = "output.shp"
CSV_FILENAME = "exp_data.csv"

r = csv.reader(open(CSV_FILENAME, 'r'), delimiter=',', quotechar=None)
header = dict(((str, i) for i, str in enumerate(r.next())))
#print header

# load data rows into memory
rows = [row for row in r]


# sort by x if larger than y. otherwise sort by y
np_rows = np.array(rows)
np_rows = np.delete(np_rows,header['pt_dt'],axis=1)
np_rows = np_rows.astype(np.float)

rows.sort(key=lambda x: x[0] if np_rows.ptp(axis=0)[0] > np_rows.ptp(axis=0)[1] else x[1])


# Create new shapefile
ogr.UseExceptions()
sr = osr.SpatialReference()
sr.ImportFromEPSG(900913)
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(SHP_FILENAME)

layer = ds.CreateLayer("data", sr, ogr.wkbLineString25D)

# Create a new line geometry
line = ogr.Geometry(type=ogr.wkbLineString25D)
line.SetCoordinateDimension(3)



# Add GPS points to line
lon_idx, lat_idx ,Z_idx = header['X'], header['Y'], header['Z']
for row in rows:
    line.AddPoint(float(row[lon_idx]),float(row[lat_idx]),float(row[Z_idx]))

# Add line as a new feature to the shapefile
feature = ogr.Feature(feature_def=layer.GetLayerDefn())
feature.SetGeometryDirectly(line)

layer.CreateFeature(feature)

# Cleanup
feature.Destroy()
ds.Destroy()

If the cross section line is North-South with some tolerance, your sorting will fail. You better look at maxX-minX and maxY-minY, and take the axis of the larger extent for sorting. — AndreJ, Aug 27 '14 at 09:28

What is an efficient Algorithm to draw a straight line from points?

1 Answers1