How the 'perspective division' should be done?

x and y should be divided by w, not z. If you're implementing a depth buffer, z is also divided by w, and the result is your depth buffer value.

Just to extend what Hodgman said. After you apply your projection matrix. all components are divided by w which clamps your x and y between -1 to 1 and z from 0 to 1 in DirectX and -1 to 1 in GL.

Just wanted to add in some links that may be useful, all tough not exactly relevant to SDL, I found this to be an easy to follow description of the perspective divide http://learnopengl.com/#!Getting-started/Coordinate-Systems , and if you wish to delve into the maths a bit more, the following is worth a look, https://unspecified.wordpress.com/2012/06/21/calculating-the-gluperspective-matrix-and-other-opengl-matrix-maths/ . As said above, it looks like your error is dividing by z, which is the wrong value and probably producing some bizare results based solely on where the items are placed on the scene.

Frist of all, thank you for your answers.

It still doesn't work when I divide by w, so I'm going to show some of the code I think is relevant:

//View matrix creation
void CreateViewMatrix(Mat4 m, Vec3 cameraTarget, Vec3 cameraPos, Vec3 up)
{
	Vec3 w;   
	Sub3(cameraTarget, cameraPos, w);
	Normalize3(w);

	Vec3 u;   
	Cross3Prod(up, w, u);
	Normalize3(u);

	Vec3 v;  
	Cross3Prod(w, u, v);

	m[0][0] = u[0];		m[0][1] = v[0];		m[0][2] = w[0];		m[0][3] = 0;
	m[1][0] = u[1];		m[1][1] = v[1];		m[1][2] = w[1];		m[1][3] = 0;
	m[2][0] = u[2];		m[2][1] = v[2];		m[2][2] = w[2];		m[2][3] = 0;

	m[3][0] = -Dot3Prod(cameraPos, u);
	m[3][1] = -Dot3Prod(cameraPos, v);
	m[3][2] = -Dot3Prod(cameraPos, w);
	m[3][3] = 1;
}

//Projection matrix creation
void CreateProjectionMatrix(Mat4 m, float nearPlane, float farPlane, float fieldOfView, float aspectRatio)
{
	m[0][0] = 1/(aspectRatio * tanf(fieldOfView / 2));
	m[0][1] = m[0][2] = m[0][3] = 0;
	m[1][1] = 1/tanf(fieldOfView / 2);
	m[1][0] = m[1][2] = m[1][3] = 0;
	m[2][2] = farPlane / (farPlane - nearPlane);
	m[2][3] = 1;
	m[2][0] = m[2][1] = 0;
	m[3][2] = -nearPlane * farPlane / (farPlane - nearPlane);
	m[3][0] = m[3][1] = m[3][3] = 0;

//Matrix multiplication. Works with Vec3 because the vertex data comes in that way, it works as it the fourth component was 1.
void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 r)
{
	r[0] = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0] + m[3][0];
	r[1] = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1] + m[3][1];
	r[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2] + m[3][2];
	r[3] = v[0] * m[0][3] + v[1] * m[1][3] + v[2] * m[2][3] + m[3][3];
}

//Convert vertex data in arr (every element is a Vec3) to pixel. The res[3] element is the w component of arr[j] * wvp
Vec4 res;
Vec3Mat4Product(arr[j], wvp, res);
float x = ((res[0] / res[3]) + SCREEN_WIDTH / 2) / SCREEN_WIDTH; 
float y = ((res[1] / res[3]) + SCREEN_HEIGHT / 2) / SCREEN_HEIGHT;
                
int px = floor(x * SCREEN_WIDTH);
int py = floor((1 - y) * SCREEN_HEIGHT);

putpixel(screenSurface, px, py, 0xFF0000);

The Vec3, Vec4, Mat4 are just alias for float[3], float[4] and float[4][4].

I don't show the vectors of the view matrix I used because I've tried with so many values that I think the problem isn't there.

The basic math functions like Normalize, DotProduct, etc were checked many times so I'm not showing their code.

Thanks again.

void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 r) change to void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 & r) or even better

void Vec3Mat4Product(Vec3 v, Mat4 m, float * r)

same for void CreateProjectionMatrix(Mat4 m, you need to pass reference or a pointer to the var because you wont write into the variable.

and

Vec4 res;

Vec3Mat4Product(arr[j], wvp, res);

//it looks like nonsense that below

float x = ((res[0] / res[3]) + SCREEN_WIDTH / 2) / SCREEN_WIDTH;
float y = ((res[1] / res[3]) + SCREEN_HEIGHT / 2) / SCREEN_HEIGHT;

int px = floor(x * SCREEN_WIDTH);
int py = floor((1 - y) * SCREEN_HEIGHT);

putpixel(screenSurface, px, py, 0xFF0000);

try that:

Vec4 res;

Vec3Mat4Product(arr[j], wvp, res);

res[0] = res[0] / res[3];

res[1] = res[1] / res[3];

res[2] = res[2] / res[3];

for (int i=0; i < 3; i++)

res = res * 0.5 + 0.5;

putpixel(screenSurface, int(SCREEN_WIDTH*res[0]+0.5), int(SCREEN_HEIGHT*res[1]+0.5), 0xFF0000);

despite wrong parameters in your functions i cant compare the code in them till feburary :P

void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 r) change to void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 & r) or even better
void Vec3Mat4Product(Vec3 v, Mat4 m, float * r)

same for void CreateProjectionMatrix(Mat4 m, you need to pass reference or a pointer to the var because you wont write into the variable.

The code I'm writing is in C, so I don't need to pass references or pointers to arrays because they already are pointers. I mean, Vec4 is a typedef of float[4] and when I pass res I'm passing &res[0] and the values of the elements change.

When I use this

Vec4 res;
Vec3Mat4Product(arr[j], wvp, res);

res[0] = res[0] / res[3];
res[1] = res[1] / res[3];
res[2] = res[2] / res[3];

for (int i=0; i < 3; i++)
res = res * 0.5 + 0.5;

putpixel(screenSurface, int(SCREEN_WIDTH*res[0]+0.5), int(SCREEN_HEIGHT*res[1]+0.5), 0xFF0000);

I get a segmentation fault because x or y coordinates of some pixels become negative.

The original code I used, the one you said looks like nonsense, it's from here: http://www.scratchapixel.com/lessons/3d-basic-rendering/computing-pixel-coordinates-of-3d-point/mathematics-computing-2d-coordinates-of-3d-points

ah yes exceuse me my bad:

for (int i=0; i < 3; i++) res = res * 2.0 + 1.0;

Try this:

http://www.glprogramming.com/red/appendixf.html

then

https://www.opengl.org/wiki/GluLookAt_code

make a mix of these two and then

MVP * vertex will work;

anyways

void Vec3Mat4Product(Vec3 v, Mat4 m, Vec4 r)
{
r[0] = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0] + m[3][0];
r[1] = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1] + m[3][1];
r[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2] + m[3][2];
r[3] = v[0] * m[0][3] + v[1] * m[1][3] + v[2] * m[2][3] + m[3][3];
}

that function looks like is all wrong you should write: r[0] = v[0] * m[0][0] + v[1] * m[0][1] + v[2] * m[0][2] + m[0][3]; etc

this looks like you want to do mvp * vertex where your matrix says it must be pvm * vertex;

so after fixing that to * 2.0 - 1.0;

try putting different order (Projection * View * World) * vertex; if that doesn;'t work then you im 99% sure you messed row major matrix with column major matrix order.

here is old code (dont mind comments there) of matrix class i have made i cant guarantee everything there is correct, but you want:

gllookat function

gluperspectivea function - you set fovY not horizontal fov keep that in mind so if you want hor fox you need to to something like lets say you want fov 90.0 then your fovY = 90.0 * (scrH / scrW);

multiplication of the (matrix class)

and t4dpoint operator*(const t3dpoint p) const in matrix class

this will give you ability to do v = ((WorldMatrix * ViewMatrix) * ProjectionMatrix) * vertex;

then you do v = v * 0.5 + 0.5;

and

putpixel(screenSurface, int(SCREEN_WIDTH*res[0]+0.5), int(SCREEN_HEIGHT*res[1]+0.5), 0xFF0000);

click on the spoiler to see code:

[spoiler]

//---------------------------------------------------------------------------

#ifndef DxcMatrixH
#define DxcMatrixH
#include "DxcMath.h"
#include "SysUtils.hpp"
#include "Dialogs.hpp"
#include "Math.hpp"
#include "Math.h"                                                                       
#include "gl/glew.h"
#include "gl/gl.h"
//#include "OPENGL_ATTACH.h"
#include "gl/glu.h"
//---------------------------------------------------------------------------

//Orders for opengl and math
//object transform: Scale * Rotate * Translate   opengl   T * R * S
//Shader: model * view * projection




  const int DimMat = 16;

						 

template <class T> class Matrix44 {
public:



//int numb[DimMat];
//int nula1[DimMat];
//int ipiv[DimMat];
//int indxr[DimMat];
//int indxc[DimMat];

T m[16];


~Matrix44()	{}

Matrix44()	{		LoadIdentity();		}

Matrix44(T a1, T a2, T a3, T a4,
					 T a5, T a6, T a7, T a8,
					 T a9, T a10, T a11, T a12,
					 T a13, T a14, T a15, T a16)
{


m[0] = a1;m[1] =a2; m[2] = a3; m[3] = a4;
 m[4] = a5;m[5] = a6; m[6] = a7; m[7] = a8;
 m[8] = a9;m[9] = a10; m[10] = a11; m[11] = a12;
 m[12] = a13;m[13] = a14; m[14] = a15; m[15] = a16;
}



t4dpoint<T> operator*(const t3dpoint<T> p) const
{
t4dpoint<T> vertexPos;
vertexPos.x = p.x;
vertexPos.y = p.y;
vertexPos.z = p.z;
vertexPos.w = 1.0;

t4dpoint<T> vertexClip;
t4dpoint<T> matrow;

matrow.x = m[0]; matrow.y = m[1]; matrow.z = m[2]; matrow.w = m[3];
vertexClip.x = dp4(matrow, vertexPos);

matrow.x = m[4]; matrow.y = m[5]; matrow.z = m[6]; matrow.w = m[7];
vertexClip.y = dp4(matrow, vertexPos);

matrow.x = m[8]; matrow.y = m[9]; matrow.z = m[10]; matrow.w = m[11];
vertexClip.z = dp4(matrow, vertexPos);

matrow.x = m[12]; matrow.y = m[13]; matrow.z = m[14]; matrow.w = m[15];
vertexClip.w = dp4(matrow, vertexPos);

vertexClip.x = vertexClip.x / vertexClip.w;
vertexClip.y = vertexClip.y / vertexClip.w;
vertexClip.z = vertexClip.z / vertexClip.w;


return vertexClip;
}


t4dpoint<T> operator^(const t3dpoint<T> p) const
{
t4dpoint<T> vertexPos;
vertexPos.x = p.x;
vertexPos.y = p.y;
vertexPos.z = p.z;
vertexPos.w = 1.0;

t4dpoint<T> vertexClip;
t4dpoint<T> matrow;

matrow.x = m[0]; matrow.y = m[1]; matrow.z = m[2]; matrow.w = m[3];
vertexClip.x = dp4(matrow, vertexPos);

matrow.x = m[4]; matrow.y = m[5]; matrow.z = m[6]; matrow.w = m[7];
vertexClip.y = dp4(matrow, vertexPos);

matrow.x = m[8]; matrow.y = m[9]; matrow.z = m[10]; matrow.w = m[11];
vertexClip.z = dp4(matrow, vertexPos);

matrow.x = m[12]; matrow.y = m[13]; matrow.z = m[14]; matrow.w = m[15];
vertexClip.w = dp4(matrow, vertexPos);

return vertexClip;
}
 
Matrix44<T> operator*(const T mat[16]) const
{
	return Matrix44<T>(

(m[0]*mat[0])+(m[4]*mat[1])+(m[8]*mat[2])+(m[12]*mat[3]),
(m[1]*mat[0])+(m[5]*mat[1])+(m[9]*mat[2])+(m[13]*mat[3]),
(m[2]*mat[0])+(m[6]*mat[1])+(m[10]*mat[2])+(m[14]*mat[3]),
(m[3]*mat[0])+(m[7]*mat[1])+(m[11]*mat[2])+(m[15]*mat[3]),

(m[0]*mat[4])+(m[4]*mat[5])+(m[8]*mat[6])+(m[12]*mat[7]),
(m[1]*mat[4])+(m[5]*mat[5])+(m[9]*mat[6])+(m[13]*mat[7]),
(m[2]*mat[4])+(m[6]*mat[5])+(m[10]*mat[6])+(m[14]*mat[7]),
(m[3]*mat[4])+(m[7]*mat[5])+(m[11]*mat[6])+(m[15]*mat[7]),

(m[0]*mat[8])+(m[4]*mat[9])+(m[8]*mat[10])+(m[12]*mat[11]),
(m[1]*mat[8])+(m[5]*mat[9])+(m[9]*mat[10])+(m[13]*mat[11]),
(m[2]*mat[8])+(m[6]*mat[9])+(m[10]*mat[10])+(m[14]*mat[11]),
(m[3]*mat[8])+(m[7]*mat[9])+(m[11]*mat[10])+(m[15]*mat[11]),

(m[0]*mat[12])+(m[4]*mat[13])+(m[8]*mat[14])+(m[12]*mat[15]),
(m[1]*mat[12])+(m[5]*mat[13])+(m[9]*mat[14])+(m[13]*mat[15]),
(m[2]*mat[12])+(m[6]*mat[13])+(m[10]*mat[14])+(m[14]*mat[15]),
(m[3]*mat[12])+(m[7]*mat[13])+(m[11]*mat[14])+(m[15]*mat[15]));
}


Matrix44<T> operator *(Matrix44<T> mat)
{
	return Matrix44<T>(

(m[0]*mat.m[0])+(m[4]*mat.m[1])+(m[8]*mat.m[2])+(m[12]*mat.m[3]),
(m[1]*mat.m[0])+(m[5]*mat.m[1])+(m[9]*mat.m[2])+(m[13]*mat.m[3]),
(m[2]*mat.m[0])+(m[6]*mat.m[1])+(m[10]*mat.m[2])+(m[14]*mat.m[3]),
(m[3]*mat.m[0])+(m[7]*mat.m[1])+(m[11]*mat.m[2])+(m[15]*mat.m[3]),

(m[0]*mat.m[4])+(m[4]*mat.m[5])+(m[8]*mat.m[6])+(m[12]*mat.m[7]),
(m[1]*mat.m[4])+(m[5]*mat.m[5])+(m[9]*mat.m[6])+(m[13]*mat.m[7]),
(m[2]*mat.m[4])+(m[6]*mat.m[5])+(m[10]*mat.m[6])+(m[14]*mat.m[7]),
(m[3]*mat.m[4])+(m[7]*mat.m[5])+(m[11]*mat.m[6])+(m[15]*mat.m[7]),

(m[0]*mat.m[8])+(m[4]*mat.m[9])+(m[8]*mat.m[10])+(m[12]*mat.m[11]),
(m[1]*mat.m[8])+(m[5]*mat.m[9])+(m[9]*mat.m[10])+(m[13]*mat.m[11]),
(m[2]*mat.m[8])+(m[6]*mat.m[9])+(m[10]*mat.m[10])+(m[14]*mat.m[11]),
(m[3]*mat.m[8])+(m[7]*mat.m[9])+(m[11]*mat.m[10])+(m[15]*mat.m[11]),

(m[0]*mat.m[12])+(m[4]*mat.m[13])+(m[8]*mat.m[14])+(m[12]*mat.m[15]),
(m[1]*mat.m[12])+(m[5]*mat.m[13])+(m[9]*mat.m[14])+(m[13]*mat.m[15]),
(m[2]*mat.m[12])+(m[6]*mat.m[13])+(m[10]*mat.m[14])+(m[14]*mat.m[15]),
(m[3]*mat.m[12])+(m[7]*mat.m[13])+(m[11]*mat.m[14])+(m[15]*mat.m[15]));
}
	

void operator =(Matrix44<T> mat)
{
	memcpy(m, mat.m, sizeof(T) * 16);
}

void operator =(Matrix44<T> * mat)
{
	memcpy(m, mat->m, sizeof(T) * 16);
}



void operator =(T kk[16])
{
	memcpy(m, kk, sizeof(T) * 16);
}

void operator = (t3dpoint<T> v)
{
	this->LoadIdentity();
	m[12] = v.x;
	m[13] = v.y;
	m[14] = v.z;
}

void RotateX(T angle)
{
	Matrix44<T> rot(1,          0,           0, 0,
		          0, cos(angle), -sin(angle), 0,
				  0, sin(angle),  cos(angle), 0,
				  0,          0,           0, 1);

	(*this) = (*this) * rot;
}

void RotateY(T angle)
{
	Matrix44<T> rot(cos(angle), 0, -sin(angle), 0,
		                 0, 1,           0, 0,
				sin(angle), 0,  cos(angle), 0,
				         0, 0,           0, 1);

	(*this) = (*this) * rot;
}

void RotateZ(T angle)
{
	Matrix44<T> rot(cos(angle), -sin(angle), 0, 0,
		          sin(angle),  cos(angle), 0, 0,
				           0,           0, 1, 0,
				           0,           0, 0, 1);

	(*this) = (*this) * rot;
}

void Translate(T x, T y, T z)
{
	Matrix44<T> tran( 1, 0, 0, x,
					  0, 1, 0, y,
					  0, 0, 1, z,
					  0, 0, 0, 1);
	(*this) = (*this) * tran;
}

void TranslateP(t3dpoint<T> p)
{
	Matrix44<T> tran( 1, 0, 0, p.x,
					  0, 1, 0, p.y,
					  0, 0, 1, p.z,
					  0, 0, 0, 1);
	(*this) = (*this) * tran;
}


void Scale(T x, T y, T z)      //OK
{
	Matrix44<T> sc( x, 0, 0, 0,
				   0, y, 0, 0,
				   0, 0, z, 0,
				   0, 0, 0, 1);

	(*this) = (*this) * sc;
}





void LoadIdentity()
{
m[0] = 1.0;m[1] = 0.0; m[2] = 0.0; m[3] = 0.0;
 m[4] = 0.0;m[5] = 1.0; m[6] = 0.0; m[7] = 0.0;
 m[8] = 0.0;m[9] = 0.0; m[10] = 1.0; m[11] = 0.0;
 m[12] = 0.0;m[13] = 0.0; m[14] = 0.0; m[15] = 1.0;
}

void LoadBiasIdentity()
{
	(*this) = Matrix44<T>( 	0.5, 	0.0, 	0.0, 	0.0,
							0.0, 	0.5, 	0.0, 	0.0,
							0.0, 	0.0, 	0.5, 	0.0,
							0.5, 	0.5, 	0.5, 	1.0);


                        Transpose();

}

void Return_GL_MATRIX(T * mat)
{
Matrix44<T> tmp;

tmp = (*this);

tmp.Transpose();

int i;

for (i=0; i<16; i++) mat[i] = tmp.m[i];


}




//I notice 2 things; OpenGL uses column major matrix notation, so the matrix elements looks like this;
//M[0] = s[0];
//M[4] = s[1];
//M[8] = s[2];

void LoadGLMatrix(T v[16])
{
//	[r1	u1	l1	px]
//	[r2	u2	l2	py]
//	[r3	u3	l3	pz]
//	[0	0	0	1]
	m[0] = v[ 0];         //  0  1  2  3
	m[1] = v[ 1];         //  4  5  6  7
	m[2] = v[ 2];         //  8  9  10 11
	m[3] = v[ 3];         //  12 13 14 15
	m[4] = v[ 4];
	m[5] = v[ 5];         //  0 4 8  12
	m[6] = v[ 6];         //  1 5 9  13
	m[7] = v[ 7];         //  2 6 10 14
	m[8] = v[ 8];         //  3 7 11 15
	m[9] = v[ 9];
	m[10] = v[10];
	m[11] = v[11];
	m[12] = v[12];
	m[13] = v[13];
	m[14] = v[14];
	m[15] = v[15];


	Transpose(); //for now...
}

t3dpoint<T> TransformVertex(t3dpoint<T> v)
{
	t3dpoint<T> result;

	result.x = m[0] * v.x + m[4] * v.y + m[8] * v.z + m[12];
	result.y = m[1] * v.x + m[5] * v.y + m[9] * v.z + m[13];
	result.z = m[2] * v.x + m[6] * v.y + m[10] * v.z + m[14];

	return result;
}

void Transpose()
{
T tmp[16];
memcpy(tmp, m, sizeof(T) * 16);
	m[0] = tmp[0];
	m[4] = tmp[1];
	m[8] = tmp[2];
	m[12] = tmp[3];
	m[1] = tmp[4];
	m[5] = tmp[5];
	m[9] = tmp[6];
	m[13] = tmp[7];
	m[2] = tmp[8];
	m[6] = tmp[9];
	m[10] = tmp[10];
	m[14] = tmp[11];
	m[3] = tmp[12];
	m[7] = tmp[13];
	m[11] = tmp[14];
	m[15] = tmp[15];
}




Matrix44<T> __fastcall Inverse()
{
 Matrix44<T> invOut;
  double inv[16], det;
    int i;

    inv[0] = m[5]  * m[10] * m[15] - 
             m[5]  * m[11] * m[14] - 
             m[9]  * m[6]  * m[15] + 
             m[9]  * m[7]  * m[14] +
             m[13] * m[6]  * m[11] - 
             m[13] * m[7]  * m[10];

    inv[4] = -m[4]  * m[10] * m[15] + 
              m[4]  * m[11] * m[14] + 
              m[8]  * m[6]  * m[15] - 
              m[8]  * m[7]  * m[14] - 
              m[12] * m[6]  * m[11] + 
              m[12] * m[7]  * m[10];

    inv[8] = m[4]  * m[9] * m[15] - 
             m[4]  * m[11] * m[13] - 
             m[8]  * m[5] * m[15] + 
             m[8]  * m[7] * m[13] + 
             m[12] * m[5] * m[11] - 
             m[12] * m[7] * m[9];

    inv[12] = -m[4]  * m[9] * m[14] + 
               m[4]  * m[10] * m[13] +
               m[8]  * m[5] * m[14] - 
               m[8]  * m[6] * m[13] - 
               m[12] * m[5] * m[10] + 
               m[12] * m[6] * m[9];

    inv[1] = -m[1]  * m[10] * m[15] + 
              m[1]  * m[11] * m[14] + 
              m[9]  * m[2] * m[15] - 
              m[9]  * m[3] * m[14] - 
              m[13] * m[2] * m[11] + 
              m[13] * m[3] * m[10];

    inv[5] = m[0]  * m[10] * m[15] - 
             m[0]  * m[11] * m[14] - 
             m[8]  * m[2] * m[15] + 
             m[8]  * m[3] * m[14] + 
             m[12] * m[2] * m[11] - 
             m[12] * m[3] * m[10];

    inv[9] = -m[0]  * m[9] * m[15] + 
              m[0]  * m[11] * m[13] + 
              m[8]  * m[1] * m[15] - 
              m[8]  * m[3] * m[13] - 
              m[12] * m[1] * m[11] + 
              m[12] * m[3] * m[9];

    inv[13] = m[0]  * m[9] * m[14] - 
              m[0]  * m[10] * m[13] - 
              m[8]  * m[1] * m[14] + 
              m[8]  * m[2] * m[13] + 
              m[12] * m[1] * m[10] - 
              m[12] * m[2] * m[9];

    inv[2] = m[1]  * m[6] * m[15] - 
             m[1]  * m[7] * m[14] - 
             m[5]  * m[2] * m[15] + 
             m[5]  * m[3] * m[14] + 
             m[13] * m[2] * m[7] - 
             m[13] * m[3] * m[6];

    inv[6] = -m[0]  * m[6] * m[15] + 
              m[0]  * m[7] * m[14] + 
              m[4]  * m[2] * m[15] - 
              m[4]  * m[3] * m[14] - 
              m[12] * m[2] * m[7] + 
              m[12] * m[3] * m[6];

    inv[10] = m[0]  * m[5] * m[15] - 
              m[0]  * m[7] * m[13] - 
              m[4]  * m[1] * m[15] + 
              m[4]  * m[3] * m[13] + 
              m[12] * m[1] * m[7] - 
              m[12] * m[3] * m[5];

    inv[14] = -m[0]  * m[5] * m[14] + 
               m[0]  * m[6] * m[13] + 
               m[4]  * m[1] * m[14] - 
               m[4]  * m[2] * m[13] - 
               m[12] * m[1] * m[6] + 
               m[12] * m[2] * m[5];

    inv[3] = -m[1] * m[6] * m[11] + 
              m[1] * m[7] * m[10] + 
              m[5] * m[2] * m[11] - 
              m[5] * m[3] * m[10] - 
              m[9] * m[2] * m[7] + 
              m[9] * m[3] * m[6];

    inv[7] = m[0] * m[6] * m[11] - 
             m[0] * m[7] * m[10] - 
             m[4] * m[2] * m[11] + 
             m[4] * m[3] * m[10] + 
             m[8] * m[2] * m[7] - 
             m[8] * m[3] * m[6];

    inv[11] = -m[0] * m[5] * m[11] + 
               m[0] * m[7] * m[9] + 
               m[4] * m[1] * m[11] - 
               m[4] * m[3] * m[9] - 
               m[8] * m[1] * m[7] + 
               m[8] * m[3] * m[5];

    inv[15] = m[0] * m[5] * m[10] - 
              m[0] * m[6] * m[9] - 
              m[4] * m[1] * m[10] + 
              m[4] * m[2] * m[9] + 
              m[8] * m[1] * m[6] - 
              m[8] * m[2] * m[5];

    det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];

    if (det == 0)
        return invOut;

    det = 1.0 / det;

    for (i = 0; i < 16; i++)
		invOut.m[i] = inv[i] * det;


		return invOut;

}





Matrix44<T> __fastcall InverseTranspose()
{
	Matrix44<T> result;

	T tmp[12];												//temporary pair storage
	T det;													//determinant
//
//	//calculate pairs for first 8 elements (cofactors)
	tmp[0] = m[10] * m[15];
	tmp[1] = m[11] * m[14];
	tmp[2] = m[9] * m[15];
	tmp[3] = m[11] * m[13];
	tmp[4] = m[9] * m[14];
	tmp[5] = m[10] * m[13];
	tmp[6] = m[8] * m[15];
	tmp[7] = m[11] * m[12];
	tmp[8] = m[8] * m[14];
	tmp[9] = m[10] * m[12];
	tmp[10] = m[8] * m[13];
	tmp[11] = m[9] * m[12];

	//calculate first 8 elements (cofactors)
	result.m[0] =		tmp[0]*m[5] + tmp[3]*m[6] + tmp[4]*m[7]
					-	tmp[1]*m[5] - tmp[2]*m[6] - tmp[5]*m[7];

	result.m[1]=		tmp[1]*m[4] + tmp[6]*m[6] + tmp[9]*m[7]
					-	tmp[0]*m[4] - tmp[7]*m[6] - tmp[8]*m[7];

	result.m[2]=		tmp[2]*m[4] + tmp[7]*m[5] + tmp[10]*m[7]
					-	tmp[3]*m[4] - tmp[6]*m[5] - tmp[11]*m[7];

	result.m[3]=		tmp[5]*m[4] + tmp[8]*m[5] + tmp[11]*m[6]
					-	tmp[4]*m[4] - tmp[9]*m[5] - tmp[10]*m[6];
//
	result.m[4]=		tmp[1]*m[1] + tmp[2]*m[2] + tmp[5]*m[3]
					-	tmp[0]*m[1] - tmp[3]*m[2] - tmp[4]*m[3];

	result.m[5]=		tmp[0]*m[0] + tmp[7]*m[2] + tmp[8]*m[3]
					-	tmp[1]*m[0] - tmp[6]*m[2] - tmp[9]*m[3];

	result.m[6]=		tmp[3]*m[0] + tmp[6]*m[1] + tmp[11]*m[3]
					-	tmp[2]*m[0] - tmp[7]*m[1] - tmp[10]*m[3];

	result.m[7]=		tmp[4]*m[0] + tmp[9]*m[1] + tmp[10]*m[2]
					-	tmp[5]*m[0] - tmp[8]*m[1] - tmp[11]*m[2];
//
	//calculate pairs for second 8 elements (cofactors)
	tmp[0] = m[2]*m[7];
	tmp[1] = m[3]*m[6];
	tmp[2] = m[1]*m[7];
	tmp[3] = m[3]*m[5];
	tmp[4] = m[1]*m[6];
	tmp[5] = m[2]*m[5];
	tmp[6] = m[0]*m[7];
	tmp[7] = m[3]*m[4];
	tmp[8] = m[0]*m[6];
	tmp[9] = m[2]*m[4];
	tmp[10] = m[0]*m[5];
	tmp[11] = m[1]*m[4];

//	//calculate second 8 elements (cofactors)
	result.m[8]=	tmp[0]*m[13] + tmp[3]*m[14] + tmp[4]*m[15]
					-	tmp[1]*m[13] - tmp[2]*m[14] - tmp[5]*m[15];

	result.m[9]=	tmp[1]*m[12] + tmp[6]*m[14] + tmp[9]*m[15]
					-	tmp[0]*m[12] - tmp[7]*m[14] - tmp[8]*m[15];

	result.m[10]=	tmp[2]*m[12] + tmp[7]*m[13] + tmp[10]*m[15]
					-	tmp[3]*m[12] - tmp[6]*m[13] - tmp[11]*m[15];

	result.m[11]=	tmp[5]*m[12] + tmp[8]*m[13] + tmp[11]*m[14]
					-	tmp[4]*m[12] - tmp[9]*m[13] - tmp[10]*m[14];

	result.m[12]=	tmp[2]*m[10] + tmp[5]*m[11] + tmp[1]*m[9]
					-	tmp[4]*m[11] - tmp[0]*m[9] - tmp[3]*m[10];

	result.m[13]=	tmp[8]*m[11] + tmp[0]*m[8] + tmp[7]*m[10]
					-	tmp[6]*m[10] - tmp[9]*m[11] - tmp[1]*m[8];

	result.m[14]=		tmp[6]*m[9] + tmp[11]*m[11] + tmp[3]*m[8]
					-	tmp[10]*m[11] - tmp[2]*m[8] - tmp[7]*m[9];

	result.m[15]=		tmp[10]*m[10] + tmp[4]*m[8] + tmp[9]*m[9]
					-	tmp[8]*m[9] - tmp[11]*m[10] - tmp[5]*m[8];

	// calculate determinant
	det	=	 m[0]*result.m[0]
			+m[1]*result.m[1]
			+m[2]*result.m[2]
			+m[3]*result.m[3];
//
	if(det==0.0)
	{
		result.LoadIdentity();
		return result;
	}



int i;
for (i=0; i < 16; i++)
result.m[i]= result.m[i] * ( 1.0f / det );

	return result;



}





};


































template <class T>
void gluPerspectiveA(Matrix44<double> & matrix, T fovy, T aspect, T zmin, T zmax)
{
 T xmin, xmax, ymin, ymax;
 ymax = zmin * tan(fovy * M_PI / 360.0);
 ymin = -ymax;
 xmin = ymin * aspect;
 xmax = ymax * aspect;
 //mglFrustum(m, xmin, xmax, ymin, ymax, zmin, zmax);
 //mglFrustum(Matrix44<T> & matrix, T l, T r, T b, T t, T n, T f)


matrix.m[0] = 	(2.0*zmin)/(xmax-xmin);
matrix.m[1] = 	0.0;
matrix.m[2] =  	(xmax + xmin) / (xmax - xmin);
matrix.m[3] =	0.0;

matrix.m[4] = 	0.0;
matrix.m[5] = 	(2.0*zmin) / (ymax - ymin);
matrix.m[6] =   (ymax + ymin) / (ymax - ymin);
matrix.m[7] =   0.0;

matrix.m[8] = 	0.0;
matrix.m[9] =  	0.0;
matrix.m[10] =	-(zmax + zmin) / (zmax-zmin);
matrix.m[11] =  (-2.0*zmax*zmin) / (zmax-zmin);

matrix.m[12] =  0.0;
matrix.m[13] =  0.0;
matrix.m[14] =  -1.0;
matrix.m[15] =  0.0;


}



template <class T>
mglFrustum(Matrix44<T> & matrix, T l, T r, T b, T t, T n, T f)
{

matrix.m[0] = 	(2.0*n)/(r-l);
matrix.m[1] = 	0.0;
matrix.m[2] =  	(r + l) / (r - l);
matrix.m[3] =	0.0;

matrix.m[4] = 	0.0;
matrix.m[5] = 	(2.0*n) / (t - b);
matrix.m[6] =   (t + b) / (t - b);
matrix.m[7] =   0.0;

matrix.m[8] = 	0.0;
matrix.m[9] =  	0.0;
matrix.m[10] =	-(f + n) / (f-n);
matrix.m[11] =  (-2.0*f*n) / (f-n);

matrix.m[12] =  0.0;
matrix.m[13] =  0.0;
matrix.m[14] =  -1.0;
matrix.m[15] =  0.0;

}
void gluOrthoA(Matrix44<double> & matrix, double Left, double Right, double Bottom, double Top, double znear, double zfar)
{

matrix.m[0] = 	2.0 / (Right - Left);
matrix.m[1] = 	0.0;
matrix.m[2] =  	0.0;
matrix.m[3] =	(Right + Left) / (Right - Left);

matrix.m[4] = 	0.0;
matrix.m[5] = 	2.0 / (Top - Bottom);
matrix.m[6] =   0.0;
matrix.m[7] =   (Top + Bottom) / (Top - Bottom);

matrix.m[8] = 	0.0;
matrix.m[9] =  	0.0;
matrix.m[10] =	-2.0 / (zfar - znear);
matrix.m[11] =  (zfar + znear) / (zfar - znear);

matrix.m[12] =  0.0;
matrix.m[13] =  0.0;
matrix.m[14] =  0.0;
matrix.m[15] =  1.0;
}

template <class T>
glLookAt(Matrix44<T> &matrix, t3dpoint<T> eyePosition3D,
				  t3dpoint<T> center3D, t3dpoint<T> upVector3D )
{
   t3dpoint<T>  forward, side, up;
   forward = Normalize( vectorAB(eyePosition3D, center3D) );
   side = Normalize( forward * upVector3D );
   up = side * forward;
  matrix.LoadIdentity();

	matrix.m[0] = side.x;
	matrix.m[1] = side.y;
	matrix.m[2] = side.z;

	matrix.m[4] = up.x;
	matrix.m[5] = up.y;
	matrix.m[6] = up.z;

	matrix.m[8] 	= -forward.x;
	matrix.m[9] 	= -forward.y;
	matrix.m[10] 	= -forward.z;



Matrix44<T> transgender;
transgender.Translate(-eyePosition3D.x, -eyePosition3D.y, -eyePosition3D.z);


matrix = transgender * matrix;
}

template <class T>
void SendMatrixToShader(T mat[16], int offset)
{
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset, float(mat[0]), float(mat[1]), float(mat[2]), float(mat[3]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset+1, float(mat[4]), float(mat[5]), float(mat[6]), float(mat[7]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset+2, float(mat[8]), float(mat[9]), float(mat[10]), float(mat[11]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset+3, float(mat[12]), float(mat[13]), float(mat[14]), float(mat[15]));
}



template <class T>
void SendParamToShader(t4dpoint<T> param, int offset)
{
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset, float(param.x), float(param.y), float(param.z), float(param.w));
}

template <class T>
void SendParamToShader(T a, T b, T c, T d, int offset)
{
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, offset, float(a), float(b), float(c), float(d));
}


template <class T>
void SendParamToFRAGShader(t4dpoint<T> param, int offset)
{
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, offset, float(param.x), float(param.y), float(param.z), float(param.w));
}

template <class T>
void SendParamToFRAGShader(T a, T b, T c, T d, int offset)
{
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, offset, float(a), float(b), float(c), float(d));
}


//Orders for opengl and math
//object transform: Scale * Rotate * Translate
//Shader: model * view * projection

//*****************************
//*****************************
//************NOTE*************
//*****************************
//*****************************

//this function must be replaced since opengl 4.5 doesn't use glLoadMatrix();

template <class T>
SetShaderMatrix(Matrix44<T> model, Matrix44<T> view, Matrix44<T> projection )
{

Matrix44<T> tmp;


tmp = (model * view) * projection;

glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 1, float(tmp.m[0]), float(tmp.m[1]), float(tmp.m[2]), float(tmp.m[3]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 2, float(tmp.m[4]), float(tmp.m[5]), float(tmp.m[6]), float(tmp.m[7]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 3, float(tmp.m[8]), float(tmp.m[9]), float(tmp.m[10]), float(tmp.m[11]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 4, float(tmp.m[12]), float(tmp.m[13]), float(tmp.m[14]), float(tmp.m[15]));

}


template <class T>
SendMatrixToShader(Matrix44<T> mat, int startoffset)
{
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, startoffset+0, float(mat.m[0]), float(mat.m[1]), float(mat.m[2]), float(mat.m[3]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, startoffset+1, float(mat.m[4]), float(mat.m[5]), float(mat.m[6]), float(mat.m[7]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, startoffset+2, float(mat.m[8]), float(mat.m[9]), float(mat.m[10]), float(mat.m[11]));
glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, startoffset+3, float(mat.m[12]), float(mat.m[13]), float(mat.m[14]), float(mat.m[15]));
}


template <class T>
SendMatrixToFRAGShader(Matrix44<T> mat, int startoffset)
{
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, startoffset+0, float(mat.m[0]), float(mat.m[1]), float(mat.m[2]), float(mat.m[3]));
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, startoffset+1, float(mat.m[4]), float(mat.m[5]), float(mat.m[6]), float(mat.m[7]));
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, startoffset+2, float(mat.m[8]), float(mat.m[9]), float(mat.m[10]), float(mat.m[11]));
glProgramLocalParameter4fARB(GL_FRAGMENT_PROGRAM_ARB, startoffset+3, float(mat.m[12]), float(mat.m[13]), float(mat.m[14]), float(mat.m[15]));
}


template <class T>
Matrix44<T> ReturnTranslateMatP(t3dpoint<T> p)
{
	Matrix44<T> tran( 1, 0, 0, p.x,
					  0, 1, 0, p.y,
					  0, 0, 1, p.z,
					  0, 0, 0, 1);
	return tran;
}


//
//template <class T>
//LoadMatrixGL( Matrix44<T> matrix )
//{
//double proj[16];
//double mod[16];
//glGetDoublev(GL_PROJECTION_MATRIX, proj);
//glGetDoublev(GL_MODELVIEW_MATRIX, mod);
//
//Matrix44<T> prj;
////prj = proj;
////prj.Transpose();
//Matrix44<T> model;
//
//
//Matrix44<T> world; //already identity
//Matrix44<T> tmp;
//
//
//tmp = world * matrix * prj;
//
//
//glEnable( GL_VERTEX_PROGRAM_ARB );
//glBindProgramARB( GL_VERTEX_PROGRAM_ARB, ms_vertexProgramID );
//
//glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 1, float(tmp.m[0]), float(tmp.m[1]), float(tmp.m[2]), float(tmp.m[3]));
//glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 2, float(tmp.m[4]), float(tmp.m[5]), float(tmp.m[6]), float(tmp.m[7]));
//glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 3, float(tmp.m[8]), float(tmp.m[9]), float(tmp.m[10]), float(tmp.m[11]));
//glProgramLocalParameter4fARB(GL_VERTEX_PROGRAM_ARB, 4, float(tmp.m[12]), float(tmp.m[13]), float(tmp.m[14]), float(tmp.m[15]));
//
//}



//SetMatrix( Matrix44<double> * world,  Matrix44<double> * view,  Matrix44<double> * proj)
//{
//		 Matrix44<double> tmp;
//		 tmp = (*world) * (*view)* (*proj);
//glLoadMatrixd(  tmp.m	);
//}





#endif



//void LoadBiasIdentity()
//{
//	(*this) = Matrix44<T>( 0.50, 0.0, 0.0, 0.0,
//						0.0, 0.50, 0.0, 0.0,
//						0.0, 0.0, 0.50, 0.0,
//						0.50, 0.50, 0.50, 1.0);
//}
//
//
//
//void Load1DBiasIdentity()
//{
//	(*this) = Matrix44<T>(	0.0, 0.0, 0.0, 0.0,
//						0.0, 0.0, 0.0, 0.0,
//						0.5, 0.0, 0.0, 0.0,
//						0.5, 0.0, 0.0, 1.0);
//}

[/spoiler]

The original code I used, the one you said looks like nonsense, it's from here: http://www.scratchapixel.com/lessons/3d-basic-rendering/computing-pixel-coordinates-of-3d-point/mathematics-computing-2d-coordinates-of-3d-points

Well, Wiredcat is right, particulary this:

Vec3Mat4Product(arr[j], wvp, res);
float x = ((res[0] / res[3]) + SCREEN_WIDTH / 2) / SCREEN_WIDTH;
float y = ((res[1] / res[3]) + SCREEN_HEIGHT / 2) / SCREEN_HEIGHT;

computes in x and y a value that is always nearly 0.5.

In case you want pixel coordinates (0,0=top,left) from screen space (0,0=middleofscreensquare) you would need to perform something like:

float x = (res[0] / res[3]); // screen space x
float y = (res[1] / res[3]); // screen space y

int xPix=floor((x*0.5+0.5) * SCREEN_WIDTH); // pixel coordinate/index x

int yPix=floor((-1.0*y*0.5+0.5)* SCREEN_HEIGHT);// pixel coordinate/index y