I found a screen to world coordinates function, but it's not entirely working out.

Does anything look suspicious with this code?

vertex_3 Get_Screen_Ray(vertex_3 eye, vertex_3 look_at, vertex_3 up, const float fov_degrees, const int x, const int y, const int screen_width, const int screen_height)
{
	vertex_3 E(eye.x, eye.y, eye.z);
	vertex_3 T(look_at.x, look_at.y, look_at.z);
	vertex_3 w(up.x, up.y, up.z);
	w.normalize();

	vertex_3 t = T - E;
	vertex_3 b = w.cross(t);

	vertex_3 t_n = t;
	t_n.normalize();

	vertex_3 b_n = b;
	b_n.normalize();

	vertex_3 v_n = t_n.cross(b_n);

	static const float pi = 4.0f * atanf(1.0f);

	float theta = fov_degrees * pi / 180.0f;

	float aspect = static_cast<float>(screen_width) / static_cast<float>(screen_height);
	float g_y = -tanf(theta / 2.0f);
	float g_x = g_y * aspect;

	vertex_3 q_x = b_n * (2.0f * g_x) / (static_cast<float>(screen_width) - 1.0f);
	vertex_3 q_y = v_n * (2.0f * g_y) / (static_cast<float>(screen_height) - 1.0f);

	vertex_3 p_1m = t_n - b_n * g_x - v_n * g_y;
	vertex_3 p_ij = p_1m + q_x * static_cast<float>(x) + q_y * static_cast<float>(y);

	return p_ij;
}

The whole code is at: https://github.com/sjhalayka/opengl4_stlview

I wonder why your snippet uses an up vector. There should be no need for that just to construct a world space ray form camera and mouse coords.

Maybe this helps:

	void WorldToScreen (float screen[3], float world[3])
	{
		float s[4];
		s[0] = ( world[0] * MVP[0][0] ) + ( world[1] * MVP[1][0] ) + ( world[2] * MVP[2][0]) + MVP[3][0];
		s[1] = ( world[0] * MVP[0][1] ) + ( world[1] * MVP[1][1] ) + ( world[2] * MVP[2][1]) + MVP[3][1];
		s[2] = ( world[0] * MVP[0][2] ) + ( world[1] * MVP[1][2] ) + ( world[2] * MVP[2][2]) + MVP[3][2];
		s[3] = ( world[0] * MVP[0][3] ) + ( world[1] * MVP[1][3] ) + ( world[2] * MVP[2][3]) + MVP[3][3];	

		screen[0] = s[0] / s[3] * viewportWidth/2 + viewportWidth/2;
		screen[1] = s[1] / s[3] * viewportHeight/2 + viewportHeight/2;
		screen[2] = s[2] / s[3];
	}
	
	void ScreenToWorld (float world[3], float screen[3])
	{
		float s[3];
		s[0] = 2 * screen[0] / viewportWidth - 1;
		s[1] = 2 * screen[1] / viewportHeight - 1;
		s[2] = screen[2];

		float p[4];
		p[0] = ( s[0] * MVPinv[0][0] ) + ( s[1] * MVPinv[1][0] ) + ( s[2] * MVPinv[2][0]) + MVPinv[3][0];
		p[1] = ( s[0] * MVPinv[0][1] ) + ( s[1] * MVPinv[1][1] ) + ( s[2] * MVPinv[2][1]) + MVPinv[3][1];
		p[2] = ( s[0] * MVPinv[0][2] ) + ( s[1] * MVPinv[1][2] ) + ( s[2] * MVPinv[2][2]) + MVPinv[3][2];
		p[3] = ( s[0] * MVPinv[0][3] ) + ( s[1] * MVPinv[1][3] ) + ( s[2] * MVPinv[2][3]) + MVPinv[3][3];	

		float w = 1 / p[3];
		world[0] = p[0] * w;
		world[1] = p[1] * w;
		world[2] = p[2] * w;
	}

	float WorldToScreenScale (float worldDistanceFromCameraPlane)
	{
		return worldDistanceFromCameraPlane * -2 / (projection[1][1] * viewportHeight);
	};	

I think i still use OpenGL conventions, so should work.

Given member variables are MVP matrix and its inverse, plus viewport dimensions. It takes fov from there. I think it should work also for ortho projections but not fully tested.

I use WorldToScreenScale so i can draw e.g. a sphere at given location if i want to have it a constant size in pixels. This way my transform gizmo has always the same size, no matter how distant the selected object is.

My screen coordinate also has z. If z is negative, point is behind the camera.

Thanks JoeJ, yet again.

I found another solution, which works like a charm too (uses GLM for matrix inverse):

vec3 screen_coords_to_world_coords(const int x, const int y, const int screen_width, const int screen_height)
{
	const float half_screen_width = screen_width / 2.0f;
	const float half_screen_height = screen_height / 2.0f;

	mat4 m_proj;

	for (glm::mat<4, 4, float, glm::packed_highp>::length_type i = 0; i < 4; i++)
		for (glm::mat<4, 4, float, glm::packed_highp>::length_type j = 0; j < 4; j++)
			m_proj[i][j] = main_camera.projection_mat[i * 4 + j];

	mat4 m_view;

	for (glm::mat<4, 4, float, glm::packed_highp>::length_type i = 0; i < 4; i++)
		for (glm::mat<4, 4, float, glm::packed_highp>::length_type j = 0; j < 4; j++)
			m_view[i][j] = main_camera.view_mat[i * 4 + j];

	mat4 inv_mat = inverse(m_proj * m_view);

	vec4 n((x - half_screen_width) / half_screen_width, -1 * (y - half_screen_height) / half_screen_height, -1, 1.0);
	vec4 f((x - half_screen_width) / half_screen_width, -1 * (y - half_screen_height) / half_screen_height, 1, 1.0);

	vec4 near_result = inv_mat * n;
	vec4 far_result = inv_mat * f;

	near_result /= near_result.w;
	far_result /= far_result.w;
	
	vec3 dir = vec3(far_result - near_result);

	return normalize(dir);
}