How to simulate Unity Pinhole Camera from its intrinsic parameters

We have a real camera for which its intrinsic parameters are estimated, based on the calibration toolbox of Bouget: http://www.vision.caltech.edu/bouguetj/calib_doc/
We assume that there is no a non-linear optical distortion and estimate parameters:
Focal length: The focal length in pixels is stored in the 2x1 vector fc.
Principal point: The principal point coordinates are stored in the 2x1 vector cc
as described in http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/parameters.html
According to the classic theory presented in the book of Zisserman “Multiple View Geometry in Computer Vision Second Edition” (chapter 6), the projection matrix in this case is given by:
P = [fc(1), 0, cc(1), 0;
0, fc(2), cc(2), 0;
0, 0, 1, 0]
This matrix is 3x4 and it acts on the homogenous coordinates (X,Y,Z,1)’ to get projective homogenous coordinates; (x,y, z)’ = P*(X,Y,Z,1)'. The corresponding image pixel coordinates are then calculated based on xp= x/z and yp=y/z.
Does anyone can help us to simulate exactly the same camera in the Unity. We would like to know how the classic matrix P above is related to Camera.projectionMatrix that is of dimensionality [4x4] and we would like to save the obtained rendered scene that has the same number of pixels as images from the real camera.
Or another way around, what is the mathematical expression that relates between (XYZ)-coordinates in the camera coordinate system and (xp,yp)- projective pixel coordinates and Camera.projectionMatrix of the Unity?

The same question was posted by us in http://answers.unity3d.com/questions/814701/how-to-simulate-unity-pinhole-camera-from-its-intr.html

We ourselves found a close answer (not exact) to our question second question, i.e.
“What is the mathematical expression that relates between (XYZ)-coordinates in the camera coordinate system and (xp,yp)- projective pixel coordinates and Camera.projectionMatrix of the Unity?”
In general transform consists of 3 steps:

1. Matrix multiplication using 4-dimensional homogenous coordinates:
ps = projMat*[XYZ 1]
2. Projective division: xs = ps(1)/ps(4); ys = ps(2)/ps(4)
3. Transform to pixel coordinates:
xp = W*(xs + 1)0.5;
yp = H
(1 - (ys + 1) * 0.5);

According to the same reference the procedure is the same all the time, but there can be slight differences. It could be nice to have exact description from Unity as are available from scratchpixel.

For example our test using the formulas give slightly different results in Unity and using the formulas above. Here is an image showing the difference when simulating cube of size 1 meter from z=3 meters. Red points are calculated corners using Unity Projective Matrix and mathematics above, white rectangular is a Unity image. So there is some discrepancy and it could be really nice where from the difference comes.

Hi, How did you do it?
Thanks