MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# Thread Subject: Point closest to a set of lines

 Subject: Point closest to a set of lines Date: 12 Aug, 2012 16:00:57 Message: 1 of 20 Hi, I have a set of lines (lets keep it simple to 2D). In an ideal case these should intersect at som point. Each line indicates the direction of arrival of a signal on some source in space. But the medium in which the wave travels may not be homogenous, so there may not be such a point. But there should be a point closest to all the lines, somewhere in the neighbourhood of the ideal point. Having the equations for the lines, how do I set up the (I suppose least square solution) in matlab? Thanks for listening, Kamran
 Subject: Point closest to a set of lines From: Nasser M. Abbasi Date: 12 Aug, 2012 16:12:33 Message: 2 of 20 On 8/12/2012 11:00 AM, kamran.iranpour@gmail.com wrote: > Hi, > Having the equations for the lines, how do I set up the (I suppose least square solution) in matlab? set up an Ax=b system and then write x=A\b --Nasser
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 12 Aug, 2012 16:25:21 Message: 3 of 20 On Aug 12, 6:12pm, "Nasser M. Abbasi" wrote: > On 8/12/2012 11:00 AM, kamran.iranp...@gmail.com wrote: > > > Hi, > > Having the equations for the lines, how do I set up the (I suppose least square solution) in matlab? > > set up an Ax=b system and then write > > x=A\b > > --Nasser Thanks Nasser, Is it this simple? Does "\" give the least square solution in Matlab. Well, thanks indeed. Kamran
 Subject: Point closest to a set of lines From: John D'Errico Date: 12 Aug, 2012 17:07:07 Message: 4 of 20 Kamran Iranpour wrote in message <08ae2046-ec05-4c20-9699-52b3eaf2473a@g1g2000vba.googlegroups.com>... > On Aug 12, 6:12pm, "Nasser M. Abbasi" wrote: > > On 8/12/2012 11:00 AM, kamran.iranp...@gmail.com wrote: > > > > > Hi, > > > Having the equations for the lines, how do I set up the (I suppose least square solution) in matlab? > > > > set up an Ax=b system and then write > > > > x=A\b > > > > --Nasser > > Thanks Nasser, > Is it this simple? Does "\" give the least square solution in Matlab. > Well, thanks indeed. > > Kamran Yes, A\b does give the least squares solution. John
 Subject: Point closest to a set of lines From: Bruno Luong Date: 12 Aug, 2012 19:17:07 Message: 5 of 20 As Nasser and John have hinted you, A \ b gives the least-square solution of the system Ax = b i.e., the solution that minimizes the least-square sum: sum (A*x-b).^2 . But that's only half pf the story. You need to form A and b careful so that this solution is the least-square in the geometrical sense. So be careful. If you can't figure it out, get back to us later. Bruno
 Subject: Point closest to a set of lines From: Bruno Luong Date: 13 Aug, 2012 07:41:07 Message: 6 of 20 Assuming you have n lines in R^m with a parametric form: Li = { P(:,i) + t * d(:,i) : t in R }, stored in two (m x n) arrays P and d. % For example 4 lines in R3: n = 4; d = randn(3,n); xyz = [1; 2; 3]; P = bsxfun(@plus,xyz,1e-2*randn(3,n)); % Add noise clear xyz % Here is the engine to find least-square intersection of those lines. [m n] = size(d); A = rand(n*(m-1),m); b = rand(n*(m-1),1); for k=1:n     idx = (k-1)*(m-1)+(1:m-1);     Q = null(d(:,k)')';     A(idx,:) = Q;     b(idx,1) = Q*P(:,k); end A \ b % Bruno
 Subject: Point closest to a set of lines From: Bruno Luong Date: 13 Aug, 2012 08:02:09 Message: 7 of 20 Sorry, typo correction of the code: > % Here is the engine to find least-square intersection of those lines. > [m n] = size(d); > A = zeros(n*(m-1),m); > b = zeros(n*(m-1),1); > for k=1:n > idx = (k-1)*(m-1)+(1:m-1); > Q = null(d(:,k)')'; > A(idx,:) = Q; > b(idx,1) = Q*P(:,k); > end > A \ b > > % Bruno
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 13 Aug, 2012 09:56:27 Message: 8 of 20 On Aug 13, 10:02am, "Bruno Luong" wrote: > Sorry, typo correction of the code: > > > > > > > > > % Here is the engine to find least-square intersection of those lines. > > [m n] = size(d); > > A = zeros(n*(m-1),m); > > b = zeros(n*(m-1),1); > > for k=1:n > > idx = (k-1)*(m-1)+(1:m-1); > > Q = null(d(:,k)')'; > > A(idx,:) = Q; > > b(idx,1) = Q*P(:,k); > > end > > A \ b > > > % Bruno Thanks indeed Bruno, my idea before seeking advice here was to set up the distance function from a point to a line, then differentiate the square distance with regards to the point coordinates and minimize the sum over "n" number of lines. Would you say it is the same thing ? Kamran
 Subject: Point closest to a set of lines From: Bruno Luong Date: 13 Aug, 2012 10:08:06 Message: 9 of 20 Kamran Iranpour wrote in message <2d0002ae-3bc1-4a7d-9115-2313d2470a10@w8g2000vbx.googlegroups.com>... > > Thanks indeed Bruno, > my idea before seeking advice here was to set up the distance function > from a point to a line, then differentiate the square distance with > regards to the point coordinates and minimize the sum over "n" number > of lines. Would you say it is the same thing ? > As long as it solves the same problem, then it's must be the same thing, no doubt. My code find point that minimizes the sum of the squared-distances to the lines. Bruno
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 13 Aug, 2012 13:01:56 Message: 10 of 20 On Aug 13, 12:08pm, "Bruno Luong" wrote: > Kamran Iranpour wrote in message <2d0002ae-3bc1-4a7d-9115-2313d2470...@w8g2000vbx.googlegroups.com>... > > > Thanks indeed Bruno, > > my idea before seeking advice here was to set up the distance function > > from a point to a line, then differentiate the square distance with > > regards to the point coordinates and minimize the sum over "n" number > > of lines. Would you say it is the same thing ? > > As long as it solves the same problem, then it's must be the same thing, no doubt. My code find point that minimizes the sum of the squared-distances to the lines. > > Bruno Bruno, Forgive me about the noise, but I don't quite understand the code. I tried, but the values I get does not seem correct. Could you tell me what xyz is and why is it equal to [1 2 3] in your example and why bxsfun and the plus ? I know I am acting stupid, but I am not well versed in matlab. I make three lines described in cartesian : Y1=-0.5.*X+3.5; Y2=-0.75.*X+2; Y3=0.833.*X+1.5; I set them up in vector form as in Li = { P(:,i) + t * d(:,i) : t in R }, and I expect the intersection to be at [-6.5 6 0]. But I land somewhere quite off the mark. Thanks
 Subject: Point closest to a set of lines From: Bruno Luong Date: 13 Aug, 2012 13:27:07 Message: 11 of 20 Kamran Iranpour wrote in message ... > and I expect the intersection to be at [-6.5 6 0]. Your expectation is off-mark. IMO, my post is clear enough and I don't have anything to add. You just need to follow it and eventually look at the matlab doc if the syntax escapes you. Bruno
 Subject: Point closest to a set of lines From: Torsten Date: 13 Aug, 2012 14:35:23 Message: 12 of 20 On 13 Aug., 15:01, Kamran Iranpour wrote: > On Aug 13, 12:08pm, "Bruno Luong" > wrote: > > > Kamran Iranpour wrote in message <2d0002ae-3bc1-4a7d-9115-2313d2470...@w8g2000vbx.googlegroups.com>... > > > > Thanks indeed Bruno, > > > my idea before seeking advice here was to set up the distance function > > > from a point to a line, then differentiate the square distance with > > > regards to the point coordinates and minimize the sum over "n" number > > > of lines. Would you say it is the same thing ? > > > As long as it solves the same problem, then it's must be the same thing, no doubt. My code find point that minimizes the sum of the squared-distances to the lines. > > > Bruno > > Bruno, > Forgive me about the noise, but I don't quite understand the code. I > tried, but the values I get does not seem correct. Could you tell me > what xyz is and why is it equal to [1 2 3] in your example and why > bxsfun and the plus ? I know I am acting stupid, but I am not well > versed in matlab. > I make three lines described in cartesian : > Y1=-0.5.*X+3.5; > Y2=-0.75.*X+2; > Y3=0.833.*X+1.5; > > I set them up in vector form as in Li = { P(:,i) + t * d(:,i) : t in > R }, > and I expect the intersection to be at [-6.5 6 0]. But I land > somewhere quite off the mark. > > Thanks I think you meant Y3=-0.833*x+1.5 ... Best wishes Torsten.
 Subject: Point closest to a set of lines From: Matt J Date: 13 Aug, 2012 14:50:05 Message: 13 of 20 Kamran Iranpour wrote in message ... > > I make three lines described in cartesian : > Y1=-0.5.*X+3.5; > Y2=-0.75.*X+2; > Y3=0.833.*X+1.5; =============== If you already have algebraic equations describing the lines, then you're already 90% finished. Just rearrange them in matrix form A=[1 1 1;.5, .75,.8333].'; %Note transpose b=[3.5,2,1.5].'; and solve x=A\b Bruno's code was mostly devoted to converting the lines from parametric form to algebraic form.
 Subject: Point closest to a set of lines From: Bruno Luong Date: 13 Aug, 2012 15:08:06 Message: 14 of 20 "Matt J" wrote in message ... > Kamran Iranpour wrote in message ... > Bruno's code was mostly devoted to converting the lines from parametric form to algebraic form. No, not only that, my code ensures the least-squares is on the *geometrical distances*, and this is not for fun. I though I emphasize enough in my previous posts, but obviously I failed. d = [1 1 1;      -0.5 -0.75 -0.833]; P = d; P(2,:) = P(2,:) + [3.5 2 1.5]; % Here is the engine to find least-square intersection of those lines. [m n] = size(d); A = rand(n*(m-1),m); b = rand(n*(m-1),1); for k=1:n     idx = (k-1)*(m-1)+(1:m-1);     Q = null(d(:,k)')';     A(idx,:) = Q;     b(idx,1) = Q*P(:,k); end xy = A \ b clf x = linspace(-10,10,2); P1=[-0.5 3.5]; P2=[-0.75 2]; P3=[-0.833 1.5]; plot(x,polyval(P1,x),'b-',...      x,polyval(P2,x),'g-.',...      x,polyval(P3,x),'k:',...      xy(1),xy(2),'or');  axis equal % Bruno
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 13 Aug, 2012 16:08:13 Message: 15 of 20 On Aug 13, 5:08pm, "Bruno Luong" wrote: > "Matt J" wrote in message ... > > Kamran Iranpour wrote in message ... > > Bruno's code was mostly devoted to converting the lines from parametric form to algebraic form. > > No, not only that, my code ensures the least-squares is on the *geometrical distances*, and this is not for fun. I though I emphasize enough in my previous posts, but obviously I failed. > > d = [1 1 1; > -0.5 -0.75 -0.833]; > P = d; > P(2,:) = P(2,:) + [3.5 2 1.5]; > > % Here is the engine to find least-square intersection of those lines. > [m n] = size(d); > A = rand(n*(m-1),m); > b = rand(n*(m-1),1); > for k=1:n > idx = (k-1)*(m-1)+(1:m-1); > Q = null(d(:,k)')'; > A(idx,:) = Q; > b(idx,1) = Q*P(:,k); > end > xy = A \ b > > clf > x = linspace(-10,10,2); > P1=[-0.5 3.5]; > P2=[-0.75 2]; > P3=[-0.833 1.5]; > plot(x,polyval(P1,x),'b-',... > x,polyval(P2,x),'g-.',... > x,polyval(P3,x),'k:',... > xy(1),xy(2),'or'); > axis equal > > % Bruno Thanks. I was doing some stupid "smart" coding that messed up the whole thing. You are absolutely correct. This answer is much closer to the correct answer when I test it on real data than using just A\b. But to tell the truth, I don't know the method you are using, I must go back to my text books. Thanks again. Kamran
 Subject: Point closest to a set of lines From: Matt J Date: 13 Aug, 2012 20:06:08 Message: 16 of 20 Kamran Iranpour wrote in message ... > > > Thanks. I was doing some stupid "smart" coding that messed up the > whole thing. You are absolutely correct. This answer is much closer to > the correct answer when I test it on real data than using just A\b. > But to tell the truth, I don't know the method you are using, I must > go back to my text books. ============== If I understand the approach, the orthogonal projection of a point x onto the i-th line can be written as an affine transform Pi*[x;y]+di where Pi is 2x2, di is 2x1, and i=1,2,3,... are indices of your different lines. You want to minimize the least squares cost function sum_i || Pi*[x;y]+di ||^2 which is equivalent to stacking all the Pi into a tall matrix A, all the di into a tall vector b, and finally doing xy=A\b
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 14 Aug, 2012 08:21:40 Message: 17 of 20 On Aug 13, 10:06pm, "Matt J " wrote: > Kamran Iranpour wrote in message ... > > > Thanks. I was doing some stupid "smart" coding that messed up the > > whole thing. You are absolutely correct. This answer is much closer to > > the correct answer when I test it on real data than using just A\b. > > But to tell the truth, I don't know the method you are using, I must > > go back to my text books. > > ============== > > If I understand the approach, the orthogonal projection of a point x onto the i-th line can be written as an affine transform > > Pi*[x;y]+di > > where Pi is 2x2, di is 2x1, and i=1,2,3,... are indices of your different lines. You want to minimize the least squares cost function > > sum_i || Pi*[x;y]+di ||^2 > > which is equivalent to stacking all the Pi into a tall matrix A, all the di into a tall vector b, and finally doing > > xy=A\b Thanks Matt for the explanation.
 Subject: Point closest to a set of lines From: Matt J Date: 14 Aug, 2012 10:21:07 Message: 18 of 20 Kamran Iranpour wrote in message <65268861-f66a-498f-ba35-0c5730f04208@e29g2000vbm.googlegroups.com>... > > Thanks Matt for the explanation. ============== Sure. I misstated a few things, though, so I've revised my remarks below > > If I understand the approach, the orthogonal projection of a point x onto the i-th line can be written as an affine transform > > > > Pi*[x;y]+di > > > > > where Pi is 2x2, di is 2x1, and i=1,2,3,... are indices of your different lines. ============= It's true that the orthogonal projection of xy to the line is an affine transformation, but what is actually more relevant here is that the perpendicular displacement from xy to the line is also an affine transformation: displacement  = x-(Pi*xy+di)  = (eye(2) -Pi)*xy - di  = Qi*xy - di Also, because this is 2D, we can obtain the (signed) geometric distance by taking the dot product of the displacement with the known unit vector ui perpendicular to the i-th line, which gives a simpler 1x2 affine transformation distance  = (ui*Qi)*xy- ui*di  =Ai*xy - bi The cost function to minimize is therefore sum_i || Ai*xy-bi ||^2 %sum of squared distances Stacking, we obtain A=[A1;A2;A3;...;An] b=[b1;b2;...;bn] xy=A\b
 Subject: Point closest to a set of lines From: Kamran Iranpour Date: 14 Aug, 2012 19:23:37 Message: 19 of 20 On Aug 14, 12:21pm, "Matt J " wrote: > Kamran Iranpour wrote in message <65268861-f66a-498f-ba35-0c5730f04...@e29g2000vbm.googlegroups.com>... > > > Thanks Matt for the explanation. > > ============== > > Sure. I misstated a few things, though, so I've revised my remarks below > > > > If I understand the approach, the orthogonal projection of a point x onto the i-th line can be written as an affine transform > > > > Pi*[x;y]+di > > > > where Pi is 2x2, di is 2x1, and i=1,2,3,... are indices of your different lines. > > ============= > > It's true that the orthogonal projection of xy to the line is an affine transformation, but what is actually more relevant here is that the perpendicular displacement from xy to the line is also an affine transformation: > > displacement > = x-(Pi*xy+di) > = (eye(2) -Pi)*xy - di > = Qi*xy - di > > Also, because this is 2D, we can obtain the (signed) geometric distance by taking the dot product of the displacement with the known unit vector ui perpendicular to the i-th line, which gives a simpler 1x2 affine transformation > > distance > = (ui*Qi)*xy- ui*di > =Ai*xy - bi > > The cost function to minimize is therefore > > sum_i || Ai*xy-bi ||^2 %sum of squared distances > > Stacking, we obtain > > A=[A1;A2;A3;...;An] > b=[b1;b2;...;bn] > > xy=A\b Thanks again Matt. If I understand correctly, the "A\b" on the set of line equations finds the best (in least squares sense) vertical distance (ie. only 'y' component is optimized) but not the shortest distance (orthogonal projection), right? Kamran
 Subject: Point closest to a set of lines From: Bruno Luong Date: 14 Aug, 2012 19:35:08 Message: 20 of 20 Kamran Iranpour wrote in message > Thanks again Matt. If I understand correctly, the "A\b" on the set of > line equations finds the best (in least squares sense) vertical > distance (ie. only 'y' component is optimized) Yes, that's what I wrote in my post #5 > but not the shortest > distance (orthogonal projection), right? The shortest distance is called http://en.wikipedia.org/wiki/Total_least_squares Bruno