Projection centrale sur linesegments

Question

J'ai le problème suivant. Donné sont deux ensembles de données dans l'espace 2D. Les deux ensembles de données sont des points de mesure 2D à partir d'un point central donné dans l'exemple par (0,0). J'ai besoin de projeter centralement les points du premier ensemble de données (x1, y1) sur des segments de ligne définis par le second ensemble de données en utilisant une interpolation linéaire entre les 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 ont déjà mis en oeuvre le calcul de l'intersection de la ligne classique pour toutes les lignes qui relient le set1 avec le point source et les segments de ligne de la set2. Malheureusement, la boucle for Nasted n'est pas très efficace.

Existe-t-il un moyen d'accélérer l'exécution de cet algorithme? Peut-être une implémentation vectorisée?

Des idées? Merci d'avance.

---

OK. Laissez-moi redéfinir le problème.

I ont le code ci-dessous:

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

Les projets centraux coder les points dans set1 sur la courbe qui est définie par les points en set2 par interpolation linéaire.

Answer 1

J'ai réussi à résoudre mon problème. Voici la solution si quelqu'un en aurait besoin à l'avenir. Il fonctionne beaucoup mieux:

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