Suppose you are given a set of points in R^3
that mostly live in a plane (e.g. X-Y plane + noise). We can use the SVD to find a basis for the points and project then into a 2D subspace:
function proj2D(points)
# https://math.stackexchange.com/a/99317
X = reduce(hcat, points)
μ = sum(X, dims=2) / size(X, 2)
Z = X .- μ
U = svd(Z).U
u = U[:, 1]
v = U[:, 2]
[[z ⋅ u, z ⋅ v] for z in eachcol(Z)]
end
The problem is that the basis (u, v)
is not necessarily oriented "upwards" with respect to the 3D space. For example, if the points lie in the X-Y plane, it may be the case that the basis (u,v)
produces a normal direction pointing downwards in the direction u × v = [0,0,-1]
.
Do you know how one could post-process the basis (u,v)
safely to make sure that it has a normal always pointing "upwards" according to the right-hand rule?
I guess the easiest thing to do would just be to do something like
u2 = copy(u)
v2 = copy(v]
for i in eachindex(u, v)
if (u[i] × v[i])[3] < 0
u2[i] = v[i]
v2[i] = u[i]
end
end
that is, just check if the z component is negative, and if so, switch their order
Can we generalize this idea to planes other than the X-Y plane?
For example, if the original points mostly live in Y-Z, then the right hand rule says that the "positive" basis is Y-Z, not Z-Y
yeah, let n
be the normal vector to that plane, then replace the if
statement with
if (u[i] × v[i]) ⋅ n < 0
and then swap the elements of the basis?
yeah
@Mason Protter the code u[i] and v[i] return scalars, can you please double check?
I don't have any runnable code here, so you're going to have to fill in the blanks. I though your u
and v
were a vector of vectors
I see now they're just one single vector each
so you don't need the loop at all then, just
if ((u × v) ⋅ n) < 0
u, v = v, u
end
to swap them
ok, and how do we get the right n
to use in the code?
given an arbitrary set of points, we don't know a priori if they live in X-Y, Y-Z or Z-X
oh you don't know the n
ahead of time? I might be mistaken but I'm not sure that makes sense then
Pretty sure you need the whole basis to determine if it's a left handed or right handed system
aren't left/right handed systems just conventions? iirc, to transform one to the other you just swap consistently swap x/z coordinates
whether a coordinate space is left/right handed is orthogonal to whether they point "upwards", as far as I understand it
i.e. if you have u
pointing in the y
direction and v
pointing in the x
direction, then u × v = [0, 0, -1]
. It's then arbitrary to say whether u
and v
are a left handed pair spanning an upwards oriented plane, or whether they're a right handed pair spanning a plane oriented downwards
right - whether it's left or right handed depends on your interpretation
In a N-dimensional space, you need all N
basis vectors to determine the handedness of the coordinate system
if all you want to do is make them point "upwards", the simplest way I know is to take the dotprodut between the vector in question and your up vector and flip the signs of your "up" coordinate if the dotproduct is negative, no?
if I understand your example correctly, the n
in @Mason Protters code is that desired "up" vector
The quantity of interest is going to be the pseudoscalar u₁ ∧ u₂ ∧ ... ∧ uₙ
which in the special case of 3D reduces to (u₁ × u₂) ⋅ u₃
yes
Yes, I am trying to basically fix this test here about orientation of polygonal chains that have 3D vertices: https://github.com/JuliaGeometry/Meshes.jl/blob/35de99a9b3b4887e9d70c67245647bb9f1ffb0f6/test/polytopes.jl#L145
What I am doing is projecting the 3D points on the 2D plane that best approximates them and checking if the orientation is clockwise or counter-clockwise.
The problem is that the projection of 3D points onto the plane sometimes has a "negative" basis and the orientation of the chain is reverted.
As you can see in the test, I just added a z coordinate to all points and the orientation changed.
that sounds like you're "viewing" the object from the front in one case, and the back in the other
exactly, and that is not my intent, I want to always view the object from the "same upwards side" according to the right-hand rule.
Perhaps it is safe to always use n = [1,1,1]
? It is pointing towards the positive orthant.
I will proceed and fix n=[0,0,1]
in this specific application.
what you deem "up" is entirely dependent on your conventions
popular choices are [0 0 1]
and [0 1 0]
thank you @Sukera and @Mason Protter for the help :smile:
Júlio Hoffimann said:
Perhaps it is safe to always use
n = [1,1,1]
? It is pointing towards the positive orthant.
what if u × v
is orthogonal to that?
good point, I will just stick to 0,0,1 which is the convention we need in our use case
(the orthogonality issue with your "up" can happen in every case)
though it sounds like those are the cases you don't need to worry about for your application, since that would mean the vector is already in your preferred plane
exactly
one thing I think we forgot to mention - if you have an object with some points, and some of those points lie on the "wrong" side of your plane, you don't actually know whether you're viewing the object from the front or the back without knowing which side of the object is the front. Also, if you only swap the coordinates for some of the points instead of all of the points, you'll end up "folding" the far side of your object "upwards"
Last updated: Nov 06 2024 at 04:40 UTC