Central projection on linesegments

Question

I have the following problem. Given are two data sets in 2D space. The two data sets are 2D measurement points from a central point given in the example by (0,0). I need to central project the points of the first data set (x1, y1) on line segments defined by the second data set using linear interpolation between the points.

import numpy as np
import matplotlib.pyplot as plt

x1 = np.arange(-50, 50, 1)
y1 = 50+np.random.rand(100)*5

x2 = np.arange(-20, 30, 50.0/320)
y2 = 30+np.random.rand(320)*0.5

plt.plot(x1, y1, '*',   
         x2, y2, 'x',
         0.0, 0.0, 'o') 

plt.show()

I have already implemented the classic line intersection calculation for all lines connecting the set1 with the source point and the line segments of the set2. Unfortunately the nasted for loop is not very efficient.

Is there a way to get this algorithm run faster? Maybe a vectorized implementation?

Any ideas? Thank you in advance.

OK. Let me redefine the problem.

I have the following code:

import numpy as np
import matplotlib.pyplot as plt
import time

set1 = np.zeros((100, 2))
set2 = np.zeros((320, 2))
set3 = np.zeros((100, 2))

# p1
set1[:, 0] = np.arange(-50, 50, 1)
set1[:, 1] = 50+np.random.binomial(5, 0.4, size=100)

# p2 and p3
set2[:, 0] = np.arange(-20, 50, 70.0/320)
set2[:, 1] = 30+np.random.binomial(8, 0.25, size=320)

# p0
sp = np.array([0.0, 0.0])

tstamp = time.time()

# building line direction vectors 
s1 = set1                   # set 1 is already the direction vector as sp=[0,0]
s2 = set2[1:] - set2[0:-1]  # set 2 direction vector (p3-p2)

projected_profile = np.zeros((100, 2))

# project set1 on set2
for i in range(np.size(s1)/2):
    intersect_points = np.zeros((100, 2))
    ts = np.zeros(100)
    ind1 = 0
    for j in range(np.size(s2)/2):
        # calculate line intersection
        div = s1[i, 0] * s2[j, 1] - s2[j, 0] * s1[i, 1]
        s = (s1[i, 1] * set2[j, 0] - s1[i, 0] * set2[j, 1]) / div
        t = (s2[j, 1] * set2[j, 0] - s2[j, 0] * set2[j, 1]) / div

        # check wether we are still on the line segments
        if (s>=0 and s<=1 and t>=0 and t <=1):
            intersect_points[ind1, :] = t * s1[i]
            ts[ind1] = t
            ind1 += 1

    # take the intersection with maximal distance from source point (sp)
    if ts.sum()>0:
        projected_profile[i, :] = intersect_points[np.argmax(ts), :]

print time.time()-tstamp

plt.plot(set1[:, 0], set1[:, 1], '*',
         set2[:, 0], set2[:, 1], '-',
         projected_profile[:, 0], projected_profile[:, 1], 'x',
         sp[0], sp[1], 'o')

plt.show()

The code central projects the points in set1 on the curve which is defined by the points in set2 using linear interpolation.

Under line segments I mean that in my measurement the second data set is a 2D section of a surface (of course not in this random data set). If I would connect the points with lines I would get the discrete section of my surface. — bdvd, Mar 01 '16 at 19:21
So the convex hull of your second dataset is where you want to project your points? — Juan Leni, Mar 01 '16 at 21:08
Could you clarify how you define the "line segments" in set2? In the example there are only points. I am nor sure how you are using the data — Juan Leni, Mar 01 '16 at 21:34
Hope that the edited question describes more clear what I meant. Sorry for the incorrect phrasing. — bdvd, Mar 02 '16 at 18:34

score 0 · Accepted Answer · answered Mar 03 '16 at 16:04

I managed to solve my problem. Here is the solution if somebody would need it in the future. It runs much better:

import numpy as np
import matplotlib.pyplot as plt
import time

set1 = np.zeros((100, 2))
set2 = np.zeros((320, 2))
set3 = np.zeros((100, 2))

# p1
set1[:, 0] = np.arange(-50, 50, 1)
set1[:, 1] = 50+np.random.binomial(5, 0.4, size=100)

# p2 and p3
set2[:, 0] = np.arange(-20, 50, 70.0/320)
set2[:, 1] = 30+np.random.binomial(8, 0.25, size=320)

# p0
sp = np.array([0.0, 0.0])

tstamp = time.time()

# building line direction vectors 
s1 = set1                   # set 1 is already the direction vector as sp=[0,0]
s2 = set2[1:] - set2[0:-1]  # set 2 direction vector (p3-p2)

num_master = np.size(s1, axis=0)
num_measure = np.size(s2, axis=0)

# calculate intersection
div = np.transpose(np.repeat([s1[:, 0]], num_measure, axis=0)) * s2[:, 1] - \
      np.transpose(np.transpose(np.repeat([s2[:, 0]], num_master, axis=0)) * s1[:, 1])

s = np.transpose(np.repeat([s1[:, 1]], num_measure, axis=0)) * set2[:-1, 0] - \
    np.transpose(np.repeat([s1[:, 0]], num_measure, axis=0)) * set2[:-1, 1]
s = s / div

t = s2[:, 1] * set2[:-1, 0] - s2[:, 0] * set2[:-1, 1]
t = t / div

# get results by masking invalid results
mask = np.bitwise_and(np.bitwise_and(s>=0, s<=1),
                      np.bitwise_and(t>=0, t<=1))

# mask indices
ind1 = mask.sum(1)>0
t[np.invert(mask)] = 0
ind2 = np.argmax(t[ind1], axis=1)

# calculate result
projected_profile = s1[ind1] * np.transpose(np.repeat([t[ind1, ind2]], 2, axis=0))

print time.time()-tstamp

plt.plot(set1[:, 0], set1[:, 1], '*',
         set2[:, 0], set2[:, 1], '-',
         projected_profile[:, 0], projected_profile[:, 1], 'x',
         sp[0], sp[1], 'o')

plt.show()

Central projection on linesegments

1 Answers1