2016-03-01 43 views
0

我遇到以下問題。給定的是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() 

example profiles

我已經執行用於連接與源點和SET2的線段的SET1的所有行的經典線交點計算。不幸的是,nassed for循環效率不高。

有沒有辦法讓這個算法運行得更快?也許是矢量化的實現?

任何想法?先謝謝你。


好的。讓我重新定義這個問題。

我有以下代碼:這是由點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() 

enter image description here

代碼中央項目點集1。

+0

什麼是你描述的'線段'? – purpletentacle

+0

在線段下方,我的意思是在我的測量中,第二個數據集是表面的2D部分(當然不是在這個隨機數據集中)。如果我用線連接點,我會得到我的表面的離散部分。 – bdvd

+0

因此,第二個數據集的凸包是您想要投影點的位置? – purpletentacle

回答

0

我設法解決我的問題。如果將來有人需要這個解決方案,那麼這是一個解決方案。它運行得更好:

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