Stream: helpdesk (published)

Topic: SVD and orientation

Júlio Hoffimann (Jun 19 2023 at 19:11):

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?

Mason Protter (Jun 19 2023 at 19:20):

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

Mason Protter (Jun 19 2023 at 19:21):

that is, just check if the z component is negative, and if so, switch their order

Júlio Hoffimann (Jun 19 2023 at 19:21):

Can we generalize this idea to planes other than the X-Y plane?

Júlio Hoffimann (Jun 19 2023 at 19:22):

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

Mason Protter (Jun 19 2023 at 19:23):

yeah, let n be the normal vector to that plane, then replace the if statement with

if (u[i] × v[i]) ⋅ n < 0

Júlio Hoffimann (Jun 19 2023 at 19:24):

and then swap the elements of the basis?

Mason Protter (Jun 19 2023 at 19:24):

yeah

Júlio Hoffimann (Jun 19 2023 at 19:26):

@Mason Protter the code u[i] and v[i] return scalars, can you please double check?

Mason Protter (Jun 19 2023 at 19:27):

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

Mason Protter (Jun 19 2023 at 19:27):

I see now they're just one single vector each

Mason Protter (Jun 19 2023 at 19:28):

so you don't need the loop at all then, just

if ((u × v) ⋅ n) < 0
    u, v = v, u
end

to swap them

Júlio Hoffimann (Jun 19 2023 at 19:33):

ok, and how do we get the right n to use in the code?

Júlio Hoffimann (Jun 19 2023 at 19:33):

given an arbitrary set of points, we don't know a priori if they live in X-Y, Y-Z or Z-X

Mason Protter (Jun 19 2023 at 19:39):

oh you don't know the n ahead of time? I might be mistaken but I'm not sure that makes sense then

Mason Protter (Jun 19 2023 at 19:40):

Pretty sure you need the whole basis to determine if it's a left handed or right handed system

Sukera (Jun 19 2023 at 19:41):

aren't left/right handed systems just conventions? iirc, to transform one to the other you just swap consistently swap x/z coordinates

Sukera (Jun 19 2023 at 19:42):

whether a coordinate space is left/right handed is orthogonal to whether they point "upwards", as far as I understand it

Mason Protter (Jun 19 2023 at 19:42):

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

Sukera (Jun 19 2023 at 19:43):

right - whether it's left or right handed depends on your interpretation

Mason Protter (Jun 19 2023 at 19:43):

In a N-dimensional space, you need all N basis vectors to determine the handedness of the coordinate system

Sukera (Jun 19 2023 at 19:46):

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?

Sukera (Jun 19 2023 at 19:47):

if I understand your example correctly, the n in @Mason Protters code is that desired "up" vector

Mason Protter (Jun 19 2023 at 19:48):

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₃

Sukera (Jun 19 2023 at 19:48):

yes

Júlio Hoffimann (Jun 19 2023 at 19:48):

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

Júlio Hoffimann (Jun 19 2023 at 19:49):

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.

Júlio Hoffimann (Jun 19 2023 at 19:49):

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.

Júlio Hoffimann (Jun 19 2023 at 19:51):

As you can see in the test, I just added a z coordinate to all points and the orientation changed.

Sukera (Jun 19 2023 at 19:52):

that sounds like you're "viewing" the object from the front in one case, and the back in the other

Júlio Hoffimann (Jun 19 2023 at 19:55):

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.

Júlio Hoffimann (Jun 19 2023 at 20:01):

Perhaps it is safe to always use n = [1,1,1]? It is pointing towards the positive orthant.

Júlio Hoffimann (Jun 19 2023 at 20:21):

I will proceed and fix n=[0,0,1] in this specific application.

Sukera (Jun 19 2023 at 20:49):

what you deem "up" is entirely dependent on your conventions

Sukera (Jun 19 2023 at 20:49):

popular choices are [0 0 1] and [0 1 0]

Júlio Hoffimann (Jun 19 2023 at 20:52):

thank you @Sukera and @Mason Protter for the help :smile:

Mason Protter (Jun 19 2023 at 20:56):

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?

Júlio Hoffimann (Jun 19 2023 at 20:57):

good point, I will just stick to 0,0,1 which is the convention we need in our use case

Sukera (Jun 19 2023 at 20:58):

(the orthogonality issue with your "up" can happen in every case)

Sukera (Jun 19 2023 at 20:58):

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

Júlio Hoffimann (Jun 19 2023 at 21:00):

exactly

Sukera (Jun 20 2023 at 04:09):

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: Oct 02 2023 at 04:34 UTC