Figured it out (this is the author of the question). It was just a typo, I didn't calculate the plane's normal in the right direction... It's not like I said it to be (v0 X v1) it was accidentally (v1 X v0). You can see it in the code. So, if the plane is directed the other way, you can't project on it.

Because it's a silly mistake like that. For all you wanting to know the math behind this matrix projection on a plane, I'll try to explain it:

Before I begin, I assume a little rough knowledge of linear algebra. Not much is needed.

Let's say we have a plane with normal N, with a point Q on it. For every point P, (P-Q).N=0 (that's a dot product) if (and only if) P is on the plane, since the vectors (P-Q) and N are perpendicular.

Now let's assume we also have a (fixed spot light) point S. We would like to project a point P on our plane from that spot light. This means we would like to find the point R, that is somewhere on the line defined by the points P and S, and that is also on the plane. In other words, find a scalar t, such that S+t(P-S)=R , such that R is on the plane. (P-S) is the direction vector from the spot light through the point P. We "walk" on that vector a certain amount t, starting from the point S until we land on the plane.

From 2 paragraphs ago, we learned a nice trick to know if a point is on the plane or not. So if we apply that on R we get that R is on the plane if (and only if):

N.(R-Q)=0 N.R-N.Q=0 N.R=N.Q N.(S+t(P-S))=N.Q N.S+tN.(P-S)=N.Q t=(N.Q-N.S)/(N.(P-S))

Now, if we put this back in the definition of R:

R=S+(N.Q-N.S)*(1/(N.(P-S))*(P-S)

Let's define N.(P-S) to be k

kR=(N.(P-S))*S+(N.Q-N.S)*P-(N.Q-N.S)*S kR=(N.P)*S+(N.Q-N.S)*P-(N.Q-N.S)*S-(N.S)*S kR=(N.P)*S+(N.Q-N.S)*P-(N.Q)*S

Let's remind ourselves what we know, and what we don't know, and what we want to know. We know N and Q and S. P is given to us, and we would like to find R. In other words, we would like to express R given P, and using N, Q and S. Let's continue breaking this down a bit further,

kR=(N_x*P_x+N_y*P_y+N_z*P_z)*S+(N.Q-N.S)*P-(N.Q)*S

R is point, so let's define each of its coordinates, with the coordinates of P (and also S because we also have him on the right side of the equation there).

kR_x=[N_x*S_x+(N.Q-N.S)]P_x+[N_y*S_x]P_y+[N_z*S_x]P_z-(N.Q)*S_x kR_y=[N_x*S_y]P_x+[N_y*S_y+(N.Q-N.S)]P_y+[N_z*S_y]P_z-(N.Q)*S_y kR_z=[N_x*S_z]P_x+[N_y*S_z]P_y+[N_z*S_z+(N.Q-N.S))]P_z-(N.Q)*S_z

It may seem like we didn't get anything, since we got that k on our left side, and we still need to divide by it (and that k is defined by P none the less!). Not to worry, because OpenGL uses four element vectors and not three. That last fourth element is used for translation matrices and perspective depth interpolation. For our needs right now, all we should know is that openGL divides the coordinates of each vertex by it's fourth element. This means is that that dreaded k is our fourth element. we get:

R_w=k R_w=N.(P-S) R_w=N.P-N.S R_w=[N_x]P_x+[N_y]P_y+[N_z]P_z-N.S

Alright, we defined our R through P, using N, S and Q. Let's put that in a matrix M. We want:

M*P=R

So,

M= N_x*S_x + (N.Q-N.S), N_y*S_x, N_z*S_x, -(N.Q)*S_x N_x*S_y, N_y*S_y + (N.Q-N.S), N_z*S_y, -(N.Q)*S_y N_x*S_z, N_y*S_z, N_z*S_z + (N.Q-N.S), -(N.Q)*S_z N_x, N_y, N_z, -(N.S)

While staring at this, remember that since P is a point it's fourth element is 1. (!=0, to be precise but we can assume it's a device normalized vertex)

About the plane equation. The plane equation is a vector, which has its first three elements as the normal. And its fourth element is its signed distance from the origin. Another way to calculate the distance of the plane from the origin is:

Given a point Q on the plane with normal N, the plane's distance from the origin is |N.Q|

Pretty simple, right? This is correct since :

N.Q=|N|*|Q|*cos(N,Q)

and |N|=1 , giving us:

N.Q=|Q|*cos(N,Q)=|Q|*d/|Q|=d

where d is the distance to the plane from the origin. Or, it's the size of the vector N; and N's magnitude is the size of the distance to the plane from the origin. You can see that by drawing the plane, choosing some point Q on the plane, drawing the normal N reaching out to the plane from the origin, and looking at the triangle made from the lines made by the two vectors and the plane.

In the above matrix replace -N.Q with d (the last element in the plane equation) and you're done. (d = -N.Q). That matrix given point P, will project it onto the plane defined by N and Q, from the spot light defined by S.

Hope that teaches you something new. If I made a mistake, comment and I'll fix it.