Центральная проекция на линиях

Question

У меня есть следующая проблема. Это два набора данных в 2D пространстве. Два набора данных представляют собой точки 2D измерения из центральной точки, приведенной в примере, по (0,0). Мне нужно центральное проектирование точек первого набора данных (x1, y1) на отрезках линии, определенных вторым набором данных, используя линейную интерполяцию между точками.

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

Я уже реализован расчет классической линии пересечения всех линий, соединяющих set1 с точкой источника и отрезками на set2. К сожалению, назойливый цикл не очень эффективен.

Есть ли способ получить этот алгоритм быстрее? Может быть, векторная реализация?

Любые идеи? Заранее спасибо.

---

ОК. Позвольте мне переопределить проблему.

У меня есть следующий код:

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

код центральные проекты точек в set1 на кривой, которая определяется точками в set2 с использованием линейной интерполяции.

Answer 1

Мне удалось решить мою проблему. Вот решение, если кому-то понадобится это в будущем. Он работает намного лучше:

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