Projecting a point on a plane through a matrix

Question

I need to render some shadows in opengl, one way to do this is to render your object twice, a first time multiplying it by a special "shadow matrix" that flat your object on a plane generating the shadow and a second time for the object itself.

To create the shadow matrix, a 4x4 matrix in an 16-array by rows, I pass the plane, on which the shadow is going to lie, and the light position.

I created a little program in Java to test it, but its not working. I checked the passages tens of times, but I didnt find anything, so I guess it is a logic error

public class test {

    int x = 0;
    int y = 1;
    int z = 2;
    int w = 3;
    float floor[][] = {
        {-1.0f, -1.0f, 0.0f},
        {1.0f, -1.0f, 0.0f},
        {1.0f, 1.0f, 0.0f},
        {-1.0f, 1.0f, 0.0f}};
    private float shadow_floor[] = new float[16];
    float light_position[] = {0.0f, 0.0f, 10.0f, 0.0f};

    public test() {
        //Find floorplane based on thre known points
        float plane_floor[] = calculatePlane(floor[1], floor[2], floor[3]);

        //store shadowMatrix for floor
        shadow_floor = shadowMatrix(plane_floor, light_position);

        float[] point = new float[]{1.0f, 0.0f, 5.0f, 1.0f};

        float[] projectedPoint = pointFmatrixF(point, shadow_floor);

        System.out.println("plane: (" + plane_floor[x] + ", " + plane_floor[y] + ", "
                + plane_floor[z] + "," + plane_floor[w]);
        System.out.println("point: (" + point[x] + ", " + point[y] + ", " + point[z] + ", "
                + point[w] + ")");
        System.out.println("projectedPoint: (" + projectedPoint[x] + ", " + projectedPoint[y]
                + ", " + projectedPoint[z] + ", " + projectedPoint[w] + ")");
    }

    public static void main(String args[]) {
        test test = new test();
    }

    // make shadow matrix
    public float[] shadowMatrix(float plane[], float light_pos[]) {
        float shadow_mat[] = new float[16];
        float dot;


        dot = plane[x] * light_pos[x] + plane[y] * light_pos[y]
                + plane[z] * light_pos[z] + plane[w] * light_pos[w];

        shadow_mat[0] = dot - light_pos[x] * plane[x];
        shadow_mat[4] = -light_pos[x] * plane[y];
        shadow_mat[8] = -light_pos[x] * plane[z];
        shadow_mat[12] = -light_pos[x] * plane[3];

        shadow_mat[1] = -light_pos[y] * plane[x];
        shadow_mat[5] = dot - light_pos[y] * plane[y];
        shadow_mat[9] = -light_pos[y] * plane[z];
        shadow_mat[13] = -light_pos[y] * plane[w];

        shadow_mat[2] = -light_pos[z] * plane[x];
        shadow_mat[6] = -light_pos[z] * plane[y];
        shadow_mat[10] = dot - light_pos[z] * plane[z];
        shadow_mat[14] = -light_pos[z] * plane[w];

        shadow_mat[3] = -light_pos[w] * plane[x];
        shadow_mat[7] = -light_pos[w] * plane[y];
        shadow_mat[11] = -light_pos[w] * plane[z];
        shadow_mat[15] = dot - light_pos[w] * plane[w];

        return shadow_mat;
    }

    public float[] calculatePlane(float p1[], float p2[], float p3[]) {
        //Array for planlikningen
        float plane[] = new float[4];

        //Gitt to vektorer (tre punkter) i planet kan normalen regnes ut
        //Vi vil ha aboluttverdier
        plane[x] = (((p2[y] - p1[y]) * (p3[z] - p1[z])) - ((p2[z] - p1[z])
                * (p3[y] - p1[y])));
        plane[y] = (((p2[z] - p1[z]) * (p3[x] - p1[x])) - ((p2[x] - p1[x])
                * (p3[z] - p1[z])));
        plane[z] = (((p2[x] - p1[x]) * (p3[y] - p1[y])) - ((p2[y] - p1[y])
                * (p3[x] - p1[x])));
        plane[w] = -(plane[x] * p1[x] + plane[y] * p1[y] + plane[z] * p1[z]);


        return plane;
    }

    public float[] pointFmatrixF(float[] point, float[] matrix) {
        int x = 0;
        int y = 1;
        int z = 2;
        float[] transformedPoint = new float[4];

        transformedPoint[x] =
                matrix[0] * point[x]
                + matrix[4] * point[y]
                + matrix[8] * point[z]
                + matrix[12];
        transformedPoint[y] =
                matrix[1] * point[x]
                + matrix[5] * point[y]
                + matrix[9] * point[z]
                + matrix[13];
        transformedPoint[z] =
                matrix[2] * point[x]
                + matrix[6] * point[y]
                + matrix[10] * point[z]
                + matrix[14];
        transformedPoint[w] = 1;

        return transformedPoint;
    }
}

If my light is on (0, 0, 10) and an hypotetic point P1 (that could be a vertex of a triangle) is in (1, 0, 5) then the projection of P1 on the xy plane should be (2, 0, 0), but it is not..

This is the current output

plane: (0.0, -0.0, 4.0,-0.0
point: (1.0, 0.0, 5.0, 1.0)
projectedPoint: (40.0, 0.0, 0.0, 1.0)

Where am I wrong?

Answer 1

Much theory up front, the actual bug fixing is way down there. The theory covers the general approach of how you come up with such a projection matrix and servers as a general answer to this kind of question, while the bug fixing is pretty local and specific to your implementation.

Theory

The line

Suppose you have the (variable) coordinates of your input point $P$, and the (fixed) coordinates of your light source $L$. Then any point on the line between these two points can be expressed as a linear combination $\lambda P + \mu L$. And conversely any linear combination will represent a point on that line. So for the projected point $P'$ you get the condition $$P'=\lambda P+\mu L$$

The plane

Let $F\cdot Q$ represent the dot product between the coordinates of your plane $F$ and the coordinates of a point $Q$. A point $Q$ lies on that plane if and only if $F\cdot Q=0$. So for the projected point $P'$ you get the condition $$F\cdot P'=0$$

The intersection

Now you have to combine these two: find coefficients $\lambda$ and $\mu$ such that the second equation holds. There is a very elegant way to do this, known (to some) as Plücker's $\mu$.

\begin{align*} \lambda &= F\cdot L \\ \mu &= -F\cdot P \\ P' &= (F\cdot L)P - (F\cdot P)L \end{align*}

The point chosen in this fashion will automatically lie on the plane, therefore will be your desired projected point.

As a matrix

Now you simply have to split that definition into its individual components.

\begin{align*} P'_i &= \sum_{j=1}^4 F_jL_jP_i - \sum_{j=1}^4 F_jP_jL_i & \forall i&\in\{1,2,3,4\} \\ M_{ij} &= -F_jL_i + \begin{cases} F\cdot L &\text{if }i=j \\ 0 &\text{otherwise} \end{cases} & \forall i,j&\in\{1,2,3,4\} \\ M &= \begin{pmatrix} \lambda-F_1L_1 & -F_2L_1 & -F_3L_1 & -F_4L_1 \\ -F_1L_2 & \lambda-F_2L_2 & -F_3L_2 & -F_4L_2 \\ -F_1L_3 & -F_2L_3 & \lambda-F_3L_3 & -F_4L_3 \\ -F_1L_4 & -F_2L_4 & -F_3L_4 & \lambda-F_4L_4 \end{pmatrix} & \text{with }\lambda&=F\cdot L \end{align*}

Your code

So your matrix does look reasonable. Your dot corresponds to my $\lambda$, and all the rest looks sane as well. Your problem lies with the pointFmatrixF method. There you incorrectly assume that the last coordinate of the resulting point will be 1. In this case it will be 40, so after fixing that line you get projectedPoint: (40.0, 0.0, 0.0, 40.0) which is just another homogenous representantion of the point (1, 0, 0). Simply divide all other coordinates by the last coordinate.

Also notice that your pointFmatrixF does implicitely assume point[w] == 1. This is the case for your special input, but not true in general. So you probably should multiply by the fourth component as well.

Your question claims a light source at (0, 0, 10) which would be homogenized to (0, 0, 10, 1). Your code however has a light at (0, 0, 10, 0), which is a light infinitely far away in the z direction. Fixing that as well, you get projectedPoint: (40.0, 0.0, 0.0, 20.0) which is the (2, 0, 0) you expected.