How to check the degree of parallelism
Mostra commenti meno recenti
How to check the degree parallelism of the result of intersecting planes in affine space (give the percentage of the parallelism) with tolerance to get an indicator of how much they are in the same orientation in the following case?
But from mathematical point of view not only commands.
The two objects that I want to check their parallelism are the intersecting object between plane 1 and plane 2 AND the intersected object between plane 1 and plane 3
the direction of two objects are expressed in terms of directions 1 and directions2
Note: the objects are intersecting in some points but I want an indicator that these two objects have almost the same directions!
P1=[396218721.065513 -205153846.153846 0 0 0 0
-205153846.153846 396218721.065513 -205153846.153846 0 0 0
0 -205153846.153846 396218721.065513 -205153846.153846 0 0
0 0 -205153846.153846 396218721.065513 -205153846.153846 0
0 0 0 -205153846.153846 396218721.065513 -205153846.153846
0 0 0 0 -205153846.153846 198109360.532757];
P2=[239737786.836378 -41030769.2307692 0 0 0 0
-41030769.2307692 239737786.836378 -205153846.153846 0 0 0
0 -205153846.153846 403860863.759455 -205153846.153846 0 0
0 0 -205153846.153846 403860863.759455 -205153846.153846 0
0 0 0 -205153846.153846 403860863.759455 -205153846.153846
0 0 0 0 -205153846.153846 201930431.879727];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq2=null([P2(2:6,:)-P2(1,:)])';
beq1=Aeq1*P1(1,:)';
beq2=Aeq2*P2(1,:)';
Aeq=[Aeq1;Aeq2];
beq=[beq1;beq2];
assert( rank([Aeq,beq])==rank(Aeq) , 'Hyperplanes do not intersect')
origin1 = pinv(Aeq)*beq
directions1 = null(Aeq)
P3=[376176938.886082 -184638461.538462 0 0 0 0
-184638461.538462 376176938.886082 -205153846.153846 0 0 0
0 -205153846.153846 396692323.501467 -205153846.153846 0 0
0 0 -205153846.153846 396692323.501467 -205153846.153846 0
0 0 0 -205153846.153846 396692323.501467 -205153846.153846
0 0 0 0 -205153846.153846 198346161.750733];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq3=null([P3(2:6,:)-P3(1,:)])';
beq1=Aeq1*P1(1,:)';
beq3=Aeq3*P3(1,:)';
Aeq22=[Aeq1;Aeq3];
beq22=[beq1;beq3];
assert( rank([Aeq22,beq22])==rank(Aeq) , 'Hyperplanes do not intersect')
origin2 = pinv(Aeq22)*beq22
directions2 = null(Aeq22)
4 Commenti
M
il 16 Lug 2023
M
il 16 Lug 2023
Walter Roberson
il 16 Lug 2023
@M sorry, my relevant classes were pretty much 40 years ago.
Bruno Luong
il 22 Lug 2023
This is buggy
assert( rank([Aeq22,beq22])==rank(Aeq) , 'Hyperplanes do not intersect')
The correct should be (Aeq replaced by Aeq22)
assert( rank([Aeq22,beq22])==rank(Aeq22) , 'Hyperplanes do not intersect')
Risposte (3)
The intersection of P1 and P2 is given by
x = origin1 + directions1*lambda
the intersection of P1 and P3 is given by
x = origin2 + directions2*mu
with arbitrary 4x1 vectors lambda and mu.
If the system
origin1 + directions1*lambda = origin2 + directions2*mu
has a solution for (lambda,mu), the two intersections cannot be parallel.
Let's try:
P1=[396218721.065513 -205153846.153846 0 0 0 0
-205153846.153846 396218721.065513 -205153846.153846 0 0 0
0 -205153846.153846 396218721.065513 -205153846.153846 0 0
0 0 -205153846.153846 396218721.065513 -205153846.153846 0
0 0 0 -205153846.153846 396218721.065513 -205153846.153846
0 0 0 0 -205153846.153846 198109360.532757];
P2=[239737786.836378 -41030769.2307692 0 0 0 0
-41030769.2307692 239737786.836378 -205153846.153846 0 0 0
0 -205153846.153846 403860863.759455 -205153846.153846 0 0
0 0 -205153846.153846 403860863.759455 -205153846.153846 0
0 0 0 -205153846.153846 403860863.759455 -205153846.153846
0 0 0 0 -205153846.153846 201930431.879727];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq2=null([P2(2:6,:)-P2(1,:)])';
beq1=Aeq1*P1(1,:)';
beq2=Aeq2*P2(1,:)';
Aeq=[Aeq1;Aeq2];
beq=[beq1;beq2];
assert( rank([Aeq,beq])==rank(Aeq) , 'Hyperplanes do not intersect')
origin1 = pinv(Aeq)*beq
directions1 = null(Aeq)
P3=[376176938.886082 -184638461.538462 0 0 0 0
-184638461.538462 376176938.886082 -205153846.153846 0 0 0
0 -205153846.153846 396692323.501467 -205153846.153846 0 0
0 0 -205153846.153846 396692323.501467 -205153846.153846 0
0 0 0 -205153846.153846 396692323.501467 -205153846.153846
0 0 0 0 -205153846.153846 198346161.750733];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq3=null([P3(2:6,:)-P3(1,:)])';
beq1=Aeq1*P1(1,:)';
beq3=Aeq3*P3(1,:)';
Aeq22=[Aeq1;Aeq3];
beq22=[beq1;beq3];
assert( rank([Aeq22,beq22])==rank(Aeq) , 'Hyperplanes do not intersect')
origin2 = pinv(Aeq22)*beq22
directions2 = null(Aeq22)
M = [-directions1, directions2];
rhs = origin1 - origin2;
lambdamu = pinv(M)*rhs;
lambda = lambdamu(1:4);
mu = lambdamu(5:8);
% Check that intersecting point lies in both affine spaces
intersecting_point_1 = origin1 + directions1*lambda
intersecting_point_2 = origin2 + directions2*mu
84 Commenti
I am sure that these intersecting object are parallel
This is contradictoy in itself.
Also, could you please elaborate
If the system
origin1 + directions1*lambda = origin2 + directions2*mu
has a solution for (lambda,mu), the two intersections cannot be parallel.
It means that if two affine subspaces have a common point, they cannot be parallel. And I don't know how you want to define "parallel" different from "non-intersecting" for arbitrary affine subspaces.
M
il 16 Lug 2023
Torsten
il 16 Lug 2023
I think this is not an easy problem. You can read about "Distances between affine subspaces" here:
M
il 16 Lug 2023
@Torsten to understand me more, the following are the lines that I got from 3d intersectiong plane,
Because these are lines we can know how much they are in the same orientation.
But how to express that the "directions1" and "directions2" have the same direction???
I think the distance is not a good solution of the direction!
Do we consider that "directions1" and "directions2" subspaces?
Yes. A1 = origin1 + span(directions1) and A2 = origin2 + span(directions2) are the two affine subspaces you have to compare.
Do you think that if the distance between the affine subspaces is small so they have the same orientation?
"Orientation" is very spongy. Linear (affine) subspaces don't have what is called "orientation" in mathematics. According to the presentation, at least for two vector spaces of the same dimension, their "distance" is computed according to slides 10-14. I did not read further on how this concept is generalized for affine spaces of the same dimension.
According to page 14, considering only the subspaces spanned by "directions1" and "directions2", we get
[U,S,V] = svd(directions1.'*directions2) %slice 14
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:4))) %s slice 12 and slice 14
thus a very small distance (near to being "parallel").
But you will have to check the fact that you must compare affine spaces, not the spanning vector spaces.
From my feeling it should not make a difference, but you should read this chaper in the presentation carefully.
The general solution for the intersection of the hyperplanes P1 and P2 is given by the solution of the linear system
Aeq*x = beq
with general solution
S = pinv(Aeq)*beq + null(Aeq)
(S = special solution of the inhomogeneous linear system + vector space of solutions of the homogeneous system)
The general solution for the intersection of the hyperplanes P1 and P3 is given by the solution of the linear system
Aeq22*x = beq22
with general solution
S' = pinv(Aeq22)*beq22 + null(Aeq22)
(S' = special solution of the inhomogeneous linear system + vector space of solutions of the homogeneous system)
I tried to calculate the distance between the affine spaces S= origin+directions but it always gives( 0+number i ) , the real part is always zero.
Example ?
Do you have another idea if thers is another solution other than calculating the distance? because in some cases it gives small distance but the two subspaces are not close to parallel!
Example ? And what makes you think the two subspaces are not close to parallel ?
Do you have another idea if thers is another solution other than calculating the distance?
The measure suggested in the link seems to do exactly what you want, doesn't it ?
Where did you read to add the origins to the directions ? This is certainly wrong.
Seems that for affine spaces, the concept from above has to be replaced by the "affine principle angles" from slide 28. Thus in the computation of the SVD, "directions" has to be replace by [directions,origin'/sqrt(1+norm(origin')^2);zeros(1,size(directions,2)),1/sqrt(1+norm(origin')^2)
But you have to find out how to compute "origin' " from "origin". My guess is that origin' is computed from [origin;1] which is then made orthonormal to A. But you should check this.
Torsten
il 20 Lug 2023
Look at my modified answer and the formula at the bottom of page 28 from
I think this is what the slides want to tell us about the distance of affine subspaces. I'm not certain, but if you want to test for parallelity and not for similarity of the affine subspaces, you won't have to consider the origin, have you ? So in my opinion, the first measure (only considering the vector spaces excluding the origin) was the correct one for your purpose.
P1=[396218721.065513 -205153846.153846 0 0 0 0
-205153846.153846 396218721.065513 -205153846.153846 0 0 0
0 -205153846.153846 396218721.065513 -205153846.153846 0 0
0 0 -205153846.153846 396218721.065513 -205153846.153846 0
0 0 0 -205153846.153846 396218721.065513 -205153846.153846
0 0 0 0 -205153846.153846 198109360.532757];
P2=[239737786.836378 -41030769.2307692 0 0 0 0
-41030769.2307692 239737786.836378 -205153846.153846 0 0 0
0 -205153846.153846 403860863.759455 -205153846.153846 0 0
0 0 -205153846.153846 403860863.759455 -205153846.153846 0
0 0 0 -205153846.153846 403860863.759455 -205153846.153846
0 0 0 0 -205153846.153846 201930431.879727];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq2=null([P2(2:6,:)-P2(1,:)])';
beq1=Aeq1*P1(1,:)';
beq2=Aeq2*P2(1,:)';
Aeq=[Aeq1;Aeq2];
beq=[beq1;beq2];
assert( rank([Aeq,beq])==rank(Aeq) , 'Hyperplanes do not intersect')
origin1 = pinv(Aeq)*beq
directions1 = null(Aeq)
P3=[376176938.886082 -184638461.538462 0 0 0 0
-184638461.538462 376176938.886082 -205153846.153846 0 0 0
0 -205153846.153846 396692323.501467 -205153846.153846 0 0
0 0 -205153846.153846 396692323.501467 -205153846.153846 0
0 0 0 -205153846.153846 396692323.501467 -205153846.153846
0 0 0 0 -205153846.153846 198346161.750733];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq3=null([P3(2:6,:)-P3(1,:)])';
beq1=Aeq1*P1(1,:)';
beq3=Aeq3*P3(1,:)';
Aeq22=[Aeq1;Aeq3];
beq22=[beq1;beq3];
assert( rank([Aeq22,beq22])==rank(Aeq) , 'Hyperplanes do not intersect')
origin2 = pinv(Aeq22)*beq22
directions2 = null(Aeq22)
directions1 = [directions1;zeros(1,size(directions1,2))];
origin1 = [origin1;1];
origin1_tilde = origin1;
for i = 1:size(directions1,2)
origin1_tilde = origin1_tilde - dot(directions1(:,i),origin1)*directions1(:,i);
end
origin1_tilde = origin1_tilde/norm(origin1_tilde);
directions1 = [directions1,origin1_tilde];
directions2 = [directions2;zeros(1,size(directions2,2))];
origin2 = [origin2;1];
origin2_tilde = origin2;
for i = 1:size(directions2,2)
origin2_tilde = origin2_tilde - dot(directions2(:,i),origin2)*directions2(:,i);
end
origin2_tilde = origin2_tilde/norm(origin2_tilde);
directions2 = [directions2,origin2_tilde];
[~,S,~] = svd (directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:size(S,1))))
M
il 21 Lug 2023
thanks but could you please clarify what did you mean in similarity of the affine subspaces?
Imagine 2 planes in 3d-space. As far as I understood what you asked is to decide whether the planes are parallel and the distance between them does not matter. What I mean by similar is that the planes are parallel and their distance is almost 0.
could you please clarify to me the following part of your code?Why did you compute it as this?
I computed the matrix [A,b0/sqrt(1+||b0||^2);0,1/sqrt(1+||b0||^2)] at the bottom of slide 28. All that is necessary is to apply Gram-Schmidt to make [origin1,1] orthonormal to directions1.' in order to get b0 from b.
All this is linear algebra first semester. If you have difficulties with these concepts, you should read an introductionary book on the subject first.
M
il 22 Lug 2023
could you please give me a reference that mentioned the concept of similarity? I didn't find it in the slides.
The subject of the slides (beginning at slide 25) is "distance of affine subspaces". In my opinion, this refers to what I called "similarity" (i.e. distance matters). That's why I wrote that it doesn't seem to be relevant for your case because you talked about a measure for "parallelity of affine subspaces". This concept is the subspace concept beginning with slide 3 (i.e. only taking the vectorspace into account).
Also could you please tell me where in the slides they added one in the end of the origin vector [origin1,1].
Slide 27. If the vectors in [directions1,origin1] are n-dimensional, they are embedded in (n+1)-dimensional space by [directions1,origin1;zeros(size(directions1,2),1),1].
To check whether the formulae make sense, I would take a plane in three-dimensional space written in two different ways
(x,y,z).' = origin1+directions1*(lambda,mu,zeta)
(x,y,z)' = origin2+directions2*(lambda,mu,zeta)
with directions1, directions2 being orthonormal and check whether their distance is 0 if you compute it following the code
directions1 = [directions1;zeros(1,size(directions1,2))];
origin1 = [origin1;1];
origin1_tilde = origin1;
for i = 1:size(directions1,2)
origin1_tilde = origin1_tilde - dot(directions1(:,i),origin1)*directions1(:,i);
end
origin1_tilde = origin1_tilde/norm(origin1_tilde);
directions1 = [directions1,origin1_tilde];
directions2 = [directions2;zeros(1,size(directions2,2))];
origin2 = [origin2;1];
origin2_tilde = origin2;
for i = 1:size(directions2,2)
origin2_tilde = origin2_tilde - dot(directions2(:,i),origin2)*directions2(:,i);
end
origin2_tilde = origin2_tilde/norm(origin2_tilde);
directions2 = [directions2,origin2_tilde];
[~,S,~] = svd (directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:size(S,1))))
M
il 22 Lug 2023
Could you please clarify how did you get this equation?
dimensional space by [directions1,origin1;zeros(size(directions1,2),1),1].
Slide 27. Instead of origin1, dimensions1 as affine subspace of IR^n, one works with the subspace in IR^(n+1) spanned by [dimension1;0] and [origin;0] + e_(n+1) where e_(n+1) is the (n+1)-th unit vector.
Example:
If origin1 = [5; 6] and directions1 = 1/sqrt(2) * [1; 1] spans an affine subspace of IR^2, one works with j(origin1+directions1) = span(1/sqrt(2) *[1; 1; 0], [[5 ;6 ;0]+[0 ;0 ;1]]) as a linear subspace of IR^3.
Torsten
il 22 Lug 2023
as I understand so far this is the fourth case in slide 27 "X c ∼= Gr(k + 1, n)"
If you mean slide 26, yes.
M
il 22 Lug 2023
M
il 22 Lug 2023
Torsten
il 22 Lug 2023
I looked at the result, and the result to compute the distance of two affine subspaces is the matrix product at the bottom of slide 28 and its singular value decomposition.
directions1 = [directions1;zeros(1,size(directions1,2))];
origin1 = [origin1;1];
origin1_tilde = origin1;
for i = 1:size(directions1,2)
origin1_tilde = origin1_tilde - dot(directions1(:,i),origin1)*directions1(:,i);
end
origin1_tilde = origin1_tilde/norm(origin1_tilde);
% The following "directions1" is the matrix
% [[A, b0/sqrt(1+norm(b0)^2)];[0,1/sqrt(1+norm(b0)^2)]
% at the bottom of slide 28
directions1 = [directions1,origin1_tilde];
directions2 = [directions2;zeros(1,size(directions2,2))];
origin2 = [origin2;1];
origin2_tilde = origin2;
for i = 1:size(directions2,2)
origin2_tilde = origin2_tilde - dot(directions2(:,i),origin2)*directions2(:,i);
end
origin2_tilde = origin2_tilde/norm(origin2_tilde);
% The following "directions2" is the matrix
% [[B, c0/sqrt(1+norm(c0)^2)];[0,1/sqrt(1+norm(c0)^2)]
% at the bottom of slide 28
directions2 = [directions2,origin2_tilde];
% The following svd computes the distance between the affine spaces,
% represented by "directions1" and "directions2"
[~,S,~] = svd (directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:size(S,1))))
After all the effort, maybe you can finally tell us what the underlying problem of your calculations is ?
Bruno Luong
il 23 Lug 2023
Modificato: Bruno Luong
il 23 Lug 2023
What bother me in the definition of distance from this pdf presentation paper, the distance is not invariant by translation. The world origin matter as showed in this simple example where I shift the global origin father and father, the distance get smaller instead of remains constant. May be that is desired property for the author of this paper, , some sort of distance observed by someone who seats at the world origin (of coordinates) but I'm not sure it is desired here by @M.
directions1 = [1;0];
origin1 = [0;0];
directions2 = [1;0];
origin2 = [0;1];
shifty = 2.^[0:20];
for i=1:length(shifty)
org = [0; shifty(i)];
d(i) = afinedist(directions1, origin1-org, directions2, origin2-org);
end
d
semilogx(shifty, d)
function distance = afinedist(directions1, origin1, directions2, origin2)
directions1 = [directions1;zeros(1,size(directions1,2))];
origin1 = [origin1;1];
origin1_tilde = origin1;
for i = 1:size(directions1,2)
origin1_tilde = origin1_tilde - dot(directions1(:,i),origin1)*directions1(:,i);
end
origin1_tilde = origin1_tilde/norm(origin1_tilde);
% The following "directions1" is the matrix
% [[A, b0/sqrt(1+norm(b0)^2)];[0,1/sqrt(1+norm(b0)^2)]
% at the bottom of slide 28
directions1 = [directions1,origin1_tilde];
directions2 = [directions2;zeros(1,size(directions2,2))];
origin2 = [origin2;1];
origin2_tilde = origin2;
for i = 1:size(directions2,2)
origin2_tilde = origin2_tilde - dot(directions2(:,i),origin2)*directions2(:,i);
end
origin2_tilde = origin2_tilde/norm(origin2_tilde);
% The following "directions2" is the matrix
% [[B, c0/sqrt(1+norm(c0)^2)];[0,1/sqrt(1+norm(c0)^2)]
% at the bottom of slide 28
directions2 = [directions2,origin2_tilde];
% The following svd computes the distance between the affine spaces,
% represented by "directions1" and "directions2"
[~,S,~] = svd (directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:min(size(S)))));
end
What bother me in the definition of distance from this pdf presentation paper, the distance is not invariant by translation.
Yes, that's counter-intuitive - seems to be some distance in the projective space. But as already mentionned, if @M wants to study the degree of parallelity of the affine spaces, I think considering only the underlying vector spaces that are added to the origins will suffice. For their distance, you gave "subspace" as equivalent MATLAB command to the svd solution in the presentation.
Maybe the method to determine the distance between two affine subspaces in IR^n discussed on slide 108 of
Bruno Luong
il 23 Lug 2023
Modificato: Bruno Luong
il 23 Lug 2023
This formula page 107-108 returns the distance projected on orthogonal to both afine space, but did not care how the afine space are oriented, which as I understand is the main preoccupation of @M.
The formula with SVD is fine, if one wants to remove the artefact of origin dependent, one can modify the formla of svd applied on the upper-left block and ignore the projectve part.
I believe this modification is different than MATLAB subspace. It seems MATLAB subspace considers only the largest singular value and not the whole spectrum. (But this is just an intuition, don't quote me on it)
I believe this modification is different than MATLAB subspace. It seems MATLAB subspace considers only the largest singular value and not the whole spectrum.
You mean
P1=[396218721.065513 -205153846.153846 0 0 0 0
-205153846.153846 396218721.065513 -205153846.153846 0 0 0
0 -205153846.153846 396218721.065513 -205153846.153846 0 0
0 0 -205153846.153846 396218721.065513 -205153846.153846 0
0 0 0 -205153846.153846 396218721.065513 -205153846.153846
0 0 0 0 -205153846.153846 198109360.532757];
P2=[239737786.836378 -41030769.2307692 0 0 0 0
-41030769.2307692 239737786.836378 -205153846.153846 0 0 0
0 -205153846.153846 403860863.759455 -205153846.153846 0 0
0 0 -205153846.153846 403860863.759455 -205153846.153846 0
0 0 0 -205153846.153846 403860863.759455 -205153846.153846
0 0 0 0 -205153846.153846 201930431.879727];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq2=null([P2(2:6,:)-P2(1,:)])';
beq1=Aeq1*P1(1,:)';
beq2=Aeq2*P2(1,:)';
Aeq=[Aeq1;Aeq2];
beq=[beq1;beq2];
assert( rank([Aeq,beq])==rank(Aeq) , 'Hyperplanes do not intersect');
origin1 = pinv(Aeq)*beq;
directions1 = null(Aeq);
P3=[376176938.886082 -184638461.538462 0 0 0 0
-184638461.538462 376176938.886082 -205153846.153846 0 0 0
0 -205153846.153846 396692323.501467 -205153846.153846 0 0
0 0 -205153846.153846 396692323.501467 -205153846.153846 0
0 0 0 -205153846.153846 396692323.501467 -205153846.153846
0 0 0 0 -205153846.153846 198346161.750733];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq3=null([P3(2:6,:)-P3(1,:)])';
beq1=Aeq1*P1(1,:)';
beq3=Aeq3*P3(1,:)';
Aeq22=[Aeq1;Aeq3];
beq22=[beq1;beq3];
assert( rank([Aeq22,beq22])==rank(Aeq) , 'Hyperplanes do not intersect');
origin2 = pinv(Aeq22)*beq22;
directions2 = null(Aeq22);
format long
[~,S,~] = svd(directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:size(S,1))))
subspace(directions1,directions2)
?
I want to measure the degree of parallelism to get how close the affine subspaces are (the directions+origins Or only directions).
You still didn't answer the main question:
Consider the two lines:
[x;y] = [0;0] + lambda*[0;1]
[x;y] = [5000;0] + lambda*[0;1] , lambda in IR
They are parallel, but not close. So what is your concern ? Should the calculation give a small or a high value for these two affine spaces ?
The calculation should give a small values these two affine spaces
But the Euclidean distance between them is 5000 ... This doesn't matter ? Only being parallel counts ?
but still it fails in some cases. Maybe I have to try the third method that you have suggested
@Bruno Luong is correct: this method is not suited in your case because it doesn't care about the similarity in orientation of the subspaces. E.g. you can get a distance of 0 although the suspaces are perpendicular.
Torsten
il 23 Lug 2023
I meant to say that the third method described on slide 108 of
is not suited for your purpose. I think you didn't test it yet, did you ?
I could imagine that one way to compare two affine subspaces could be to compare only directions as a first measure and the distance of the two points in the affine subspaces with minimum distance to the origin as a second measure. But I don't know how to combine these two measures because they happen on different scales (svd distance versus Euclidean distance).
Table 1 on page 19 of
gives a list of possible distances between affine subspaces. I cannot tell you in how far these metrics reflect both criteria ( orientation and distance ) and in how far they are useful for your purpose.
They are simply metrics on the set of affine subspaces of the same dimension with the three defining properties of a metric:
(1) d(x,y) >= 0 and d(x,y) = 0 <=> x = y
(2) d(x,y) = d(y,x)
(3) d(x,z) <= d(x,y) + d(y,z)
for x,y,z being affine subspaces of the same dimension.
Bruno Luong
il 24 Lug 2023
Modificato: Bruno Luong
il 24 Lug 2023
...
format long
[~,S,~] = svd(directions1.'*directions2);
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:size(S,1))))
subspace(directions1,directions2)
Yes, this makes clear the "parallel" between modify Grassmann's distance for only direction comparison (svd) and MATLAB subspace command. thanks.
Torsten
il 24 Lug 2023
When I tested the eqaution at the bottom of slide 28 it didnt give same result of directions1
I don't understand what you mean.
Also, the last column in directions1 is always [0 0 0 0 0 0 1] , How this can be an indication of the origin?
In case origin1_tilde(1:end-1) is the zero vector, your affine subspace contains the origin (i.e. origin1 is contained in span(directions1)). This means that the affine space equals the vector space.
Bruno Luong
il 31 Lug 2023
Modificato: Bruno Luong
il 31 Lug 2023
Before I move on for the record, I post here a function that compute Grassmann distance
function distance = GrassmannDist(A1, b1, A2, b2)
% Compute Grassmann distance of two afine spaces embeded in R^m of the same dimension of
% the form V = { A*x + b : x in R^n }
% where A is (m x n) matrix and b is (m x 1) vector.
% Usualy n <= dimesnion(V) = rank(A) < m.
pMat1 = GetpMat(A1, b1);
pMat2 = GetpMat(A2, b2);
s = svd(pMat1'*pMat2);
distance = norm(acos(min(s,1)));
end
function pMat = GetpMat(A, b)
A = orth(A);
b = b(:);
b0 = [b-A*(A'*b); 1];
b0 = b0/norm(b0);
A(end+1,end) = 0;
pMat = [A,b0];
end
PS: I hope @M could give a feedback or a thanks to Torsen for the lengthly thread, but no he just quits quitely
Torsten
il 31 Lug 2023
PS: I hope @M could give a feedback or a thanks to Torsen for the lengthly thread, but no he just quits quitely
He/she still has to digest https://arxiv.org/pdf/1807.10883.pdf
:-)
Torsten
il 13 Ago 2023
unfortunately I didn't find a robust solution till now, but your solution worked in some cases but in others it failed!
Well, that's life. But it's not my solution, but apparently the usual way on how to define the distance between affine subspaces of the same dimension.
M
il 13 Ago 2023
Torsten
il 13 Ago 2023
Yes it's helpful in finding the distances, but it dont give an indicator about the direction or orientation.
I think it's just the other way round as Bruno's example showed.
Also, maybe the cause of the failed cases came from that the origin in all my cases is zeros, so this makes the problem more difficult!
If the origin is zero in all your cases, why do you treat the concept of affine subspaces at all ? You have to compare two vector subspaces, and for vector subspaces the concept is pretty simple and has been elaborated at the beginning of our discussion.
Torsten
il 14 Ago 2023
Or you meant this
[U,S,V] = svd(directions1.'*directions2) %slice 14
distance = sqrt(sum(arrayfun(@(i)acos(S(i,i))^2,1:4))) %s slice 12 and slice 14
Of course.
And I vaguely remember that it was mentionned in the slices that affine and subspace approach both give the same distance if the spaces to compare contain the origin.
M
il 14 Ago 2023
Bruno Luong
il 14 Ago 2023
Modificato: Bruno Luong
il 14 Ago 2023
If the returned angle is 0 subspaces are exactly parallel, if the angle is small they are mostly parallel, if the angle is close to pi/2 (90 degree), subspace are mostly orthogonal.
It's not meaningful for you????
M
il 14 Ago 2023
Bruno Luong
il 14 Ago 2023
Modificato: Bruno Luong
il 14 Ago 2023
The "distance" in math have abstract meaning, it measures something that satisfied triangular inequality.
In the code
s = svd(pMat1'*pMat2);
distance = norm(acos(min(s,1)));
the Grassmann distance in this code returns result with angle unit, see the acos.
It has been showed manymany time by Torsen and me (see my answer) Grassmann distance and MATLAB subspace command give the SAME value
- when 2 affine spaces pass through the origin, or
- when we ignore the origin in the calculation of Grassmann formula.
This value can be interpreted as the angle between 2 subspaces so it is diretly provide the degree of parallelism between two subspaces, as you asks.
With 50+ posts I wonder why this is still not clear for you.
Bruno Luong
il 14 Ago 2023
Modificato: Bruno Luong
il 14 Ago 2023
You never explain why it doesn' work for you. You say both affine space pass through the origin then only the angle meaures how parallel they are, there is no other indicator.
If they are not goes through origin as you claim (sometime yes, somtime no, I can't keep up with your assumptions) then Grassmann distance combines the angles between 2 afine spaces AND their respective distances (true distance now) viewed from the origin.
If you prefer to separate in the affine S := { A*x+b } the orientation of span(A) and offset b, you can consider the projection of the origin to S, denoted by bp:
Q = orth(A)
bp = b-Q*(Q'*b)
Compare the orientations of two given afine spaces S1 a,d S2 using command
subspace(A1,A2) % angle in radian
and compare the offset by
norm(bp1-bp2) % distance as "viewed at the origin"
Here you have two indicators and not combined into one as Grassmann distance.
Now if that doesn't meet your need you have to explain why.
Torsten
il 14 Ago 2023
There are engineering proofs that it should be nearest object.
Then you must use the same definition of "nearest" in your computation as the engineers do.
Bruno Luong
il 14 Ago 2023
@M "In some cases it should return the smallest distance but it failed, in others it should return the nearest two or three object but it didn't."
Then simply use whatever the definition of distance and "nearest" used by this engineering proof.
Your origin question asks about "degree of parallelism" (which we already gives you the command subspace that answers exactly this question, no more no less).
I wonder now: what for YOU the relationship between parallelism and notions of distance/nearest?
And just for curiousity what this engineeing proof actually claims? This is something that you never bother to tell us (and you should) beside simply reportng "it doesn't work".
Funny, why you inquires us for a metric without the context and then wonder why it doesn't work? The engineening proof proves with the distance they are using, not with Grassmann distance, which is essentially an angle of the projective geometric consider in space of dimension 1-larger than the initial dimension.
M
il 15 Ago 2023
Walter Roberson
il 15 Ago 2023
These kinds of calculation seldom have unique solutions.
Bruno Luong
il 15 Ago 2023
Modificato: Bruno Luong
il 15 Ago 2023
The common study objects is "polytope". It can be embeded in afine space. There is nothing particular-interesting aspect about the embedding, simply to reduce the dimension.
Those objects are mathematically VERY different from the afine space.
This thread should be closed and we should move on, the orginal question of quantifying the parallelism of affine spaces has been answered very long ago.
The current discussion here sounds really like OT.
Bruno Luong
il 15 Ago 2023
Modificato: Bruno Luong
il 15 Ago 2023
It is not (interesting) to me. You still haven't give us a clue why it doesn't meet your need. All this blah blah about "distance" "nearest" "engineering proof" is just fuzzy verbal description like mud. There is no equation, not data, no code, no description about how you apply various solution proposed here, etc...
You might really want to measure a true distance between objects that reside on the affice space (my speculation), but never clearly state that to us.This is not a sincere discussion.
If you want to not disclose the details. Just do not ask question and hide weeks and work on your own. Do not expect any useful help here.
Torsten
il 15 Ago 2023
I mean the object or the space that is generated from the intersection(direction and origin) (we can't say plane right?)
The intersection of two vector spaces is a vector space; the intersection of two affine spaces is an affine space - so you have the same objects before and after your operations.
Bruno Luong
il 18 Ago 2023
Modificato: Bruno Luong
il 18 Ago 2023
The grassmann distance proposed in this answer and subspace command measure the angle between affine and subspace so they answer exactly the questin you asked.
The argument about distance and nearest is OT.
Torsten
il 18 Ago 2023
@M Imagine you ask someone for advice what you should eat and he advices you to take a cake. Then you tell the person that the cake doesn't meet your requirements, but you don't tell him that you suffer from diabetes. Do you think it makes sense continuing the discussion with this person ?
Bruno Luong
il 18 Ago 2023
Modificato: Bruno Luong
il 18 Ago 2023
This requirement is honestly is NOT in your original question which is for the record about mesuring degree of parallelism of two affine spaces (see your code).
Sahaj
il 16 Lug 2023
0 voti
Hi M.
You can use the following code to give the percentage of parallelism while including tolerance in determining if the lines can be considered parallel.
% Define the tolerance threshold
tolerance = 0.99; % Adjust this value according to your desired tolerance
% Determine the direction vectors
directions1 = null(Aeq);
directions2 = null(Aeq22);
% Compute the angle between the direction vectors
cosine_angle = dot(directions1(:, 1), directions2(:, 1)) / (norm(directions1(:, 1)) * norm(directions2(:, 1)));
angle_rad = acos(cosine_angle);
% Calculate the percentage of parallelism
angle_deg = rad2deg(angle_rad);
perpendicular_angle_deg = 90 - angle_deg;
percentage_parallelism = (1 - perpendicular_angle_deg / 90) * 100;
% Check if the lines are parallel based on the tolerance threshold
if cosine_angle >= tolerance
fprintf('Lines are parallel with a percentage of parallelism: %.2f%%\n', percentage_parallelism);
else
fprintf('Lines are not parallel.\n');
end
You can repeat the above steps to check the parallelism between plane1 and plane3.
1 Commento
Bruno Luong
il 20 Lug 2023
0 voti
1 Commento
Bruno Luong
il 22 Lug 2023
Modificato: Bruno Luong
il 22 Lug 2023
Tha angle between 2 afine spaces is 0.0303 rad or 1.7376 degrees
P1=[396218721.065513 -205153846.153846 0 0 0 0
-205153846.153846 396218721.065513 -205153846.153846 0 0 0
0 -205153846.153846 396218721.065513 -205153846.153846 0 0
0 0 -205153846.153846 396218721.065513 -205153846.153846 0
0 0 0 -205153846.153846 396218721.065513 -205153846.153846
0 0 0 0 -205153846.153846 198109360.532757];
P2=[239737786.836378 -41030769.2307692 0 0 0 0
-41030769.2307692 239737786.836378 -205153846.153846 0 0 0
0 -205153846.153846 403860863.759455 -205153846.153846 0 0
0 0 -205153846.153846 403860863.759455 -205153846.153846 0
0 0 0 -205153846.153846 403860863.759455 -205153846.153846
0 0 0 0 -205153846.153846 201930431.879727];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq2=null([P2(2:6,:)-P2(1,:)])';
beq1=Aeq1*P1(1,:)';
beq2=Aeq2*P2(1,:)';
Aeq=[Aeq1;Aeq2];
beq=[beq1;beq2];
assert( rank([Aeq,beq])==rank(Aeq) , 'Hyperplanes do not intersect')
origin1 = pinv(Aeq)*beq;
directions1 = null(Aeq)
P3=[376176938.886082 -184638461.538462 0 0 0 0
-184638461.538462 376176938.886082 -205153846.153846 0 0 0
0 -205153846.153846 396692323.501467 -205153846.153846 0 0
0 0 -205153846.153846 396692323.501467 -205153846.153846 0
0 0 0 -205153846.153846 396692323.501467 -205153846.153846
0 0 0 0 -205153846.153846 198346161.750733];
Aeq1=null([P1(2:6,:)-P1(1,:)])';
Aeq3=null([P3(2:6,:)-P3(1,:)])';
beq1=Aeq1*P1(1,:)';
beq3=Aeq3*P3(1,:)';
Aeq22=[Aeq1;Aeq3];
beq22=[beq1;beq3];
assert( rank([Aeq22,beq22])==rank(Aeq22) , 'Hyperplanes do not intersect')
origin2 = pinv(Aeq22)*beq22;
directions2 = null(Aeq22)
The angle between directions1 and directions2 (origins are irrelevant)
theta = subspace(directions1, directions2)
theta_degree = rad2deg(theta)
This is the same as the angle between the respective orthogonal subspaces
rad2deg(subspace(Aeq',Aeq22'))
Categorie
Scopri di più su Detection in Centro assistenza e File Exchange
il 16 Lug 2023
il 18 Ago 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
