0
102 Jan 11, 2010 at 11:56

Hi, consider the following problem: given a vector v=(x,y,z) located in a hyperplane. Given two reference vectors from that plane v1 and v2. Find a and b such that v = a*v1 + b*v2. How would you do that? a, b exist of course since any vector of a plane has a linear combination of 2 non-colinear vectors from the same plane. My maths are also not rusty and I generally know how to solve the system above, i.e.

x = a*x1 + b*x2
y = a*y1 + b*y2
z = a*z1 + b*z2

But what bothers me is this: to find a and b with the above system you need at most 2 of the 3 equations. But which two? Taking an equation out means you’re projecting the 3D vectors on a 2D plane, which might give a headache, if after the projection the vectors become colinear, which is very much possible. So before doing the above calculation one has to check on which plane the reference vectors are not colinear and then choose the right 2 equations (cooridnates). There’s just gotta be an easier and not that ugly way to do this I think, does anyone have an idea? Thank you!

#### 24 Replies

0
101 Jan 11, 2010 at 13:06

I’m not entirely sure and I have not tried it, but the pseudoinverse of the 3x2 matrix

[v1 v2]

should project a vector v into the coordinate system [v1,v2] ‘as good as possible’ and give you a and b.

For real-valued matrices the Moore-Penrose-Pseudoinverse is (M\^T * M)\^-1 * M\^T. (If the columns of M are linearly independent.)

0
102 Jan 11, 2010 at 13:41

Ah, interesting! Never knew about that property of the pseudo-inverse. I would definetely try it. By the way, before I sit down tonight and work out the maths, can you recommend me a good reference, describing these features of the pseudo-inverse?

edit: Are all types of pseudo-inverse matrices suitable or just the Moore-Penrose?

0
101 Jan 11, 2010 at 14:16

I just tried it in MATLAB with a few simple cases and it seems to work.

I don’t have an online reference ready, but I think about it that way:
The pseudoinverse is a way to ‘solve’ the overdetermined least squares problem

min norm(A*x - d)

x = [a;b]
A = [v1,v2]
d = v

One solution is then

x = pinv(A)*d

But again: Be careful, this is based on the fading memories of a math lecture I had years ago…

[EDIT: of course it is pinv(A)*d, not pinv(A)*b]

0
102 Jan 11, 2010 at 15:07

Ich danke dir herzlichst! :)

0
102 Jan 12, 2010 at 01:06

There’s a far more intuitive way to solve this.

Construct a vector v3, perpendicular to v1 and v2. This can be done using a cross product. Now you know that v has to be a linear combination of v1, v2, and v3, with the coefficient c for v3 in theory being 0:

x = a*x1 + b*x2 + c*x3
y = a*y1 + b*y2 + c*y3
z = a*z1 + b*z2 + c*z3

This is an ordinary system of three linear equations with three variables. Even if c isn’t exactly 0 due to numerical imprecision, you know a and b are the coefficients of the projection of v onto the plane.

0
101 Jan 12, 2010 at 10:32

You can also approach this problem geometrically. In order to get to v, you start at the origin, and then go in the direction of v1 for a certain amount (a), and from thereon in direction of v2 which points directly at v in another amount (:D. The point where you change directions is the intersection p of the line with direction v1 through the origin, and that of the line with direction v2 through v. p will then be a*v1, and (v-p) will be b*v2. Calculating a and b should be trivial.

0
102 Jan 12, 2010 at 11:11

@.oisyn

The point where you change directions is the intersection p of the line with direction v1 through the origin, and that of the line with direction v2 through v.

Problem is you can’t easily intersect two lines in 3D. You could however intersect a line and a plane, with the plane making use of v3, the cross product of v1 and v2. That makes things no easier than my solution though. The intersection point has to be computed using a linear system of equations, after which you still need to determine a and b.

0
102 Jan 12, 2010 at 11:21

@.oisyn:

I need to do this very quickly and many many BILLION times per second… The geometrical interpretation does not help in that account, but it was the way I came into these calculations in the first place.

@.nick:

The extended system seems to give a solution in this particular case, yes, but seems a bit computationally intensive, I don’t know, I have to benchmark it and I will. Also not sure if it’s ill-conditioned and/or stable. The way I was about to do it, before you posted, was by extracting the least squares approximation of the overdetermined system (as macnihist suggested), but not using the pseudoinverse (best solution, but computationally intensive). Instead I’m about to calculate it using QR decomposition. The decomposition itself is performed using Givens rotations (see the great book “Matrix Computations” by Gehe H. Golub et. al.) Let’s see if it works.

0
101 Jan 12, 2010 at 13:09

If you can cache the lengths of v1 and v2, or even better always have them normalised, cheap vector projection would do the trick. If I missed any mention of that property, well, don’t bother telling me :p

0
101 Jan 12, 2010 at 13:44

@Nick

Problem is you can’t easily intersect two lines in 3D. You could however intersect a line and a plane, with the plane making use of v3, the cross product of v1 and v2.

Yes, a line line intersection would obviously imply that one of the lines is converted into a plane.

That makes things no easier than my solution though.

I never said it was easier, it’s just another way of looking at the problem :D. It might be easier for the intuition though, as with most problems solved geometrically rather than algebraically.
@Mihail121

@.oisyn: I need to do this very quickly and many many BILLION times per second… The geometrical interpretation does not help in that account, but it was the way I came into these calculations in the first place.

Ah I see. How volatile are v, v1 and v2?

0
102 Jan 12, 2010 at 14:42

@coelurus

If you can cache the lengths of v1 and v2, or even better always have them normalised, cheap vector projection would do the trick. If I missed any mention of that property, well, don’t bother telling me :p

HANS, YOU’RE BACK!!!!!!!!!!!! I’m so glad to know you’re alive :) Yes, the reference vectors (base of the plane) are normalized.
@.oisyn

Ah I see. How volatile are v, v1 and v2?

I don’t know yet, but my (crude) guess is they are VERY volatile. We shall see after some experiments.

0
102 Jan 12, 2010 at 15:45

@Mihail121

The extended system seems to give a solution in this particular case, yes, but seems a bit computationally intensive…

It’s the least computationally intensive solution.

Also not sure if it’s ill-conditioned and/or stable.

It’s robust as long as v1 and v2 are not zero vectors and don’t coincide (i.e. they form a well-defined plane).

Yes, the reference vectors (base of the plane) are normalized.

In that case the determinant is 1, which makes computing the 3x3 matrix inverse cheaper. Note that in fact you only need 2 / 3 of the matrix inverse, since you’re not interested in the value of c.

0
101 Jan 12, 2010 at 15:58

Hi Mikhail, I have been lurking on these forums every now and then, never dared to post anything :p

To project vector A onto vector B, you do “A’ = (A.B) / (B.B) * B”. If B is normalised, you have just “A’ = A.B * B”.

0
102 Jan 12, 2010 at 18:44

@coelurus

To project vector A onto vector B, you do “A’ = (A.B) / (B.B) * B”. If B is normalised, you have just “A’ = A.B * B”.

That’s a perpendicular projection. I don’t think that’s of use here.

0
101 Jan 12, 2010 at 20:11

That was silly of me.
I have completely forgotten the instability behaviour with least squares fitting, would the normal solution really be a problem in this case?
QR decomposition definitely sounds overkill.

0
101 Jan 13, 2010 at 09:48

@Nick: Isn’t the determinant only 1 if v1 and v2 are also perpendicular?

Anyway:
Solving Nick’s system directly for a,b is probably the fastest (and best) solution, but after massaging the pseudoinverse stuff a bit further I think I have come up with a way that does not employ the ‘full construction’ of the pinv and offers a nice geometrical interpretation.
Of course this can also be seen as a solution to Nick’s linear system, depends on how you approach it.

The rows of the pseudoinverse are two normals onto which you project v:

a = dot(n1,v)
b = dot(n2,v)

you have three constraints on each normal:

dot(n1,v1) = 1
dot(n1,v2) = 0
dot(n1,cross(v1,v2)) = 0 (distances out of the plane should not count)

instead of solving for n1 and n2, in this special case you can construct them directly:

c1 = cross(v1,v2);
n1 = cross(v2,c1);
n1 *= 1/dot(v1,n1);

c2 = cross(v2,v1); // or -c1
n2 = cross(v1,c2);
n2 *= 1/dot(v2,n2);

Then you have, like written above:

a = dot(n1,v)
b = dot(n2,v)

I tried it in MATLAB with 100000 random vectors and it worked. (Let’s hope my program isn’t flawed…)

0
102 Jan 13, 2010 at 10:25

@coelurus

That was silly of me.
I have completely forgotten the instability behaviour with least squares fitting, would the normal solution really be a problem in this case?
QR decomposition definitely sounds overkill.

Yes, you were always silly :> Just don’t get too smart around that particle accelerator :D Ok, back to the topic, the normal equation you mean? I never really considered that as an option, not a single man seems to be fond of that, because of extreme instability. QR decomp. is not really an overkill if fast Givens or, hmm, “modified” Givens is used.

Of course I don’t have a problem to get the vectors normalized AND orthogonal and just transpose the coef. matrix to get the inverse, which will give me an answer instantly.

I was really hoping this thing has a small solution (1 mul, 2 adds or something), but that’s life and it’s complex.

@mac

Just tried it also, seems to work fine.

0
102 Jan 13, 2010 at 16:26

@macnihilist

Isn’t the determinant only 1 if v1 and v2 are also perpendicular?

Ah, right. The determinant of [v1 v2 (v1 x v2)] is (v1 x v2).(v1 x v2), which is the same as v3.v3. So it’s not 1 but still a very simple computation.

0
102 Jan 14, 2010 at 07:33

@Mihail121

I was really hoping this thing has a small solution (1 mul, 2 adds or something), but that’s life and it’s complex.

I’m curious why you have to do this in the first place. If v is known to be in the plane, it was constructed as a linear combination of v1 and v2 in the first place. So instead of working in world xyz coordinates, why not work in the plane’s coordinate system instead (ab coordinates)?

Can you reveal anything more about the actual application you need this for?

0
102 Jan 14, 2010 at 09:10

@Nick

I’m curious why you have to do this in the first place. If v is known to be in the plane, it was constructed as a linear combination of v1 and v2 in the first place. So instead of working in world xyz coordinates, why not work in the plane’s coordinate system instead (ab coordinates)? Can you reveal anything more about the actual application you need this for?

Sure I have nothing to hide, the project behind it is rather stupid and useless, but I got so involved in the maths it made me blind as usual :) A friend of mine got a quadcore so we wanted to brute force …ehm, test if they REALLY gave us a quadcore :) . We came up with a battery of sub-projects. One of them deals with breaking cyphers based on shift registers, another one deals with analysing data and yet another one we decided should be a software rasterizer, since I have some knowledge on the topic and need one for Silverlight anyway. I remember discussing an idea here with you (Nick et. al.), that I never found the time to test in practice, about how to crudely devide the rasterizer workload on multiple units. Can’t find the post now, but it went like this. You write a small rasterizer that fills only an ID-Z-buffer, meaning stores z-values as in a normal z-buffer, but also stores the ID of the polygon covered by that pixel. After you’ve pushed all geometry through it, you do this:

for each element Element of ID-Z-Buffer {
reverse project (Element.screen_x, Element.screen_y) using Element.z into
point.x,point.y,point.z
get triangle by Triangles[Element.ID]
v1 = triangle.reference_vector1
v2 = triangle.reference_vector2
v3 = triangle.reference_vector3
v = vector formed by point.x, point.y, point.z in the same plane as v1, v2
using the same reference point
// v1, v2, v3 precomputed, orthonormal
//
// now come the funny spot
//
calculate alpha, beta such that v=v1*alpha+v2*beta
//
// OK now comes the even funnier spot I'm not sure about at all :D
//
use alpha & beta with two precomputed vectors on the texture plane
covering the triangle plane to calculate (u1,v1)....(un,vn), the texture
coordinates of of all available textures for this point
use alpha & beta to calculate some other shading
use the available vectors for lighting & stuff
}


Now this is absurd, I know that, it’s perfectly clear to me that it’s a waste. It’s just meant to test the new hardware and we want to put it together fast. Why can we test the hardware with it? Easy, we devide the screen in four almost equal squares and use each core to rasterize a square. It should work in my eyes. The only thing I’m not sure about in the pseudocode above is using alpha and beta to extract (u,v), but it seems ok on paper. As I said, I’m so excited and blind I cannot think straight about other alternatives right now, sure they’ll be many.

0
167 Jan 14, 2010 at 18:39

Why not just go for full deferred shading? That seems to be what you’re heading toward, only in your system you have to do a bunch of extra per-pixel work to figure out triangle IDs, texture coordinates and crap…you could just rasterize out all the information for your shading equation and read it back in for the lighting pass. B)

Potentially even simpler, just write a normal forward shaded renderer with the screen split into four quadrants. This doesn’t load balance among the cores (unless you do something tricky like adjusting the size of the quadrants) but since this is just a benchmark anyway, you can choose the scene you render to be nicely balanced. :D

0
102 Jan 14, 2010 at 19:27

@Reedbeta

I want ALL of the information to be calculated in the very last loop, i.e. to be very easily deployable on a separate unit. I definetely do not want to calculate texture coordinates and other computations for pixels that are later going to be covered by other pixels. Besides, this mode of operation would allow for very simple scanline clipping than the Sutherland-Hodgeman & v-caches. The lighting is not my absolute priority anyhow.

0
167 Jan 14, 2010 at 19:46

What about storing UVs and material IDs at each pixel as you rasterize, instead of triangle IDs? You’d save yourself solving the plane equation later on, but would still be sampling the textures and doing all the lighting/shading work in the final pass.

Maybe I’ve misunderstood, but it seems like rasterizing the ID/Z-buffer is a preprocess before starting the actual benchmark? In which case, I’m not sure why you’re worried that much about computing pixels that will be later be covered by other pixels (in that phase)?

0
102 Jan 14, 2010 at 19:50

@Reedbeta

Maybe I’ve misunderstood, but it seems like rasterizing the ID/Z-buffer is a preprocess before starting the actual benchmark? In which case, I’m not sure why you’re worried that much about computing pixels that will be later be covered by other pixels (in that phase)?

Actually you’re right, but then again, not quite :) I am indeed benchmarking the last phase, but in that case it does not matter what I put in the buffers prior to it, does it? In fact, going with the ID/Z-buffer is more an act of laziness than a decision, supported by justificated claims and arguments. :)