Image center in the calibration matrix
<p>I have a <code>camera</code> and a <code>3D</code> object, I have calculated the <code>camera matrix</code>  as given in <a href="http://research.microsoft.com/en-us/um/people/zhang/Papers/Camera%20Calibration%20-%20book%20chapter.pdf" rel="nofollow noreferrer">this</a>. The image captured using the camera is having size <code>1600x1200.</code> But in the camera matrix I am not getting the value of image <code>centers</code> properly. Instead of <code>800, 600</code> I am getting some other values. What will be the possible reasons for this.</p>

<pre><code>  #include &lt;iostream&gt;
    #include "opencv2/imgproc/imgproc.hpp"
    #include "opencv2/highgui/highgui.hpp"
    #include &lt;math.h&gt; 
    using namespace std;
    using namespace cv;
    int main()
    {
      int numberofpoints=30;
      float p2d[30][2] =
  {{72, 169}, {72, 184}, {71, 264}, {96, 168}, {94, 261}, {94, 276}, {257, 133},
    {254, 322}, {254, 337}, {278, 132}, {278, 146}, {275, 321}, {439, 228},
    {439, 243}, {437, 328}, {435, 431}, {461, 226}, {459, 326}, {459, 342},
    {457, 427}, {457, 444}, {612, 196}, {609, 291}, {605, 392}, {605, 406},
    {634, 191}, {633, 206}, {630, 288}, {630, 303}, {627, 390}};



      float p3d[30][3] = 
{{0, 0, 98},  {0, 0, 94}, {0, 0, 75}, {4, 4, 98}, {4, 4, 75}, {4, 4, 71},
{0, 96, 75}, {0, 96, 25}, {0, 96, 21}, {4, 100, 75},
{4, 100, 71}, {4, 100, 25}, {96, 0, 98}, {96, 0, 94},
{96, 0, 75}, {96, 0, 50}, {100, 4, 98}, {100, 4, 75},
{100, 4, 71}, {100, 4, 50}, {100, 4, 46}, {96, 96, 75},
{96, 96, 50}, {96, 96, 25}, {96, 96, 21}, {100, 100, 75},
{100, 100, 71}, {100, 100, 50}, {100, 100, 46}, {100, 100, 25}};


      Mat Matrix2D= Mat(30, 2, CV_64FC1, &amp;p2d);
      Mat Matrix3D= Mat(30, 3, CV_64FC1, &amp;p3d);
      Mat X_3D=Matrix3D.col(0);
      Mat Y_3D=Matrix3D.col(1);
      Mat Z_3D=Matrix3D.col(2);
      Mat X_2D=Matrix2D.col(0);
      Mat Y_2D=Matrix2D.col(1);
      Mat z = Mat::zeros(numberofpoints,4, CV_64FC1);
      Mat o = Mat::ones(numberofpoints,1, CV_64FC1);
      Mat temp,temp1,C,D,AoddRows;
      hconcat( X_3D, Y_3D,temp); 
      hconcat(Z_3D ,o,temp1);
      hconcat(z, -X_2D.mul(X_3D),C);
      hconcat(-X_2D.mul(Y_3D),-X_2D.mul(Z_3D),D);
      hconcat(temp,temp1,AoddRows);
      hconcat(AoddRows,C,AoddRows);
      hconcat(AoddRows,D,AoddRows);
      hconcat(AoddRows,-X_2D,AoddRows);
      Mat A1,B1,C1,AevenRows;
      hconcat(z,Matrix3D,A1); 
      hconcat(o,-Y_2D.mul(X_3D) ,B1);
      hconcat(-Y_2D.mul(Y_3D),-Y_2D.mul(Z_3D),C1);
      hconcat(A1,B1,AevenRows);
      hconcat(AevenRows,C1,AevenRows);
      hconcat(AevenRows,-Y_2D,AevenRows);
      Mat A;                                      
      vconcat(AoddRows,AevenRows,A);              
      Mat  w, u, EigenVectors;
      SVD::compute(A, w, u, EigenVectors);
      Mat EigenVector_SmallestEigenvalue=EigenVectors.row(EigenVectors.rows-1);
      Mat P=Mat::zeros(3,4, CV_64FC1);
      int k=0;
      for(int i=0;i&lt;P.rows;i++)
      {
        for(int j=0;j&lt;P.cols;j++)
        {
          P.at&lt;double&gt;(i, j) = EigenVector_SmallestEigenvalue.at&lt;double&gt;(0,k);
          k++;
        }
      }
      double abs_lambda = sqrt(P.at&lt;double&gt;(2,0) * P.at&lt;double&gt;(2,0)  + P.at&lt;double&gt;(2,1) * P.at&lt;double&gt;(2,1) + P.at&lt;double&gt;(2,2) * P.at&lt;double&gt;(2,2));
      P = P / abs_lambda;
      Mat ProjectionMatrix = P;    
      Mat T= Mat::zeros(3,1, CV_64FC1);
      Mat R = Mat::zeros(3,3, CV_64FC1);
      Mat B = Mat::zeros(3,3, CV_64FC1);
      Mat b = Mat::zeros(3,1,CV_64FC1);
      P.row(0).colRange(0,3).copyTo(B.row(0));
      P.row(1).colRange(0,3).copyTo(B.row(1));
      P.row(2).colRange(0,3).copyTo(B.row(2));
      P.col(3).copyTo(b.col(0));
      Mat K=B*B.t();
      double Center_x = K.at&lt;double&gt;(0,2);
      double Center_y = K.at&lt;double&gt;(1,2);
      double beta = sqrt(K.at&lt;double&gt;(1,1)-Center_y*Center_y);
      double gemma = (K.at&lt;double&gt;(0,1)-Center_x*Center_y)/beta;
      double alpha = sqrt(K.at&lt;double&gt;(0,0)-Center_x*Center_x-gemma*gemma);

      Mat CameraMatrix=Mat::zeros(3,3, CV_64FC1);
      CameraMatrix.at&lt;double&gt;(0,0)=alpha;
      CameraMatrix.at&lt;double&gt;(0,1)=gemma;
      CameraMatrix.at&lt;double&gt;(0,2)=Center_x;
      CameraMatrix.at&lt;double&gt;(1,1)=beta;
      CameraMatrix.at&lt;double&gt;(1,2)=Center_y;
      CameraMatrix.at&lt;double&gt;(2,2)=1;
      R = CameraMatrix.inv()*B;
      T = CameraMatrix.inv()*b;
      Mat newMatrix2D=Mat::zeros(numberofpoints,2, CV_64FC1);
      for(int i=0;i&lt;numberofpoints;i++)
      {
        double num_u =  P.at&lt;double&gt;(0,0) * Matrix3D.at&lt;double&gt;(i,0) 
          + P.at&lt;double&gt;(0,1) * Matrix3D.at&lt;double&gt;(i,1) 
          + P.at&lt;double&gt;(0,2) * Matrix3D.at&lt;double&gt;(i,2) 
          + P.at&lt;double&gt;(0,3);

        double num_v =  P.at&lt;double&gt;(1,0) * Matrix3D.at&lt;double&gt;(i,0) 
          + P.at&lt;double&gt;(1,1) * Matrix3D.at&lt;double&gt;(i,1) 
          + P.at&lt;double&gt;(1,2) * Matrix3D.at&lt;double&gt;(i,2) 
          + P.at&lt;double&gt;(1,3);

        double den =    P.at&lt;double&gt;(2,0) * Matrix3D.at&lt;double&gt;(i,0) 
          + P.at&lt;double&gt;(2,1) * Matrix3D.at&lt;double&gt;(i,1) 
          + P.at&lt;double&gt;(2,2) * Matrix3D.at&lt;double&gt;(i,2) 
          + P.at&lt;double&gt;(2,3);

        newMatrix2D.at&lt;double&gt;(i,0)=num_u/den;
        newMatrix2D.at&lt;double&gt;(i,1)=num_v/den;
      } 

      Mat reprojMatrix=newMatrix2D;
      Mat errorDiff=reprojMatrix-Matrix2D;
      int size=errorDiff.rows;
      double sum=0;
      for(int i=0;i&lt;errorDiff.rows;i++)
      {
        sum=sum + sqrt(errorDiff.at&lt;double&gt;(i,0) * errorDiff.at&lt;double&gt;(i,0)
                       + errorDiff.at&lt;double&gt;(i,1) * errorDiff.at&lt;double&gt;(i,1));
      }
      sum=sum/size;
      cout&lt;&lt;"Average Error="&lt;&lt;sum&lt;&lt;endl;
      cout&lt;&lt;"Translation Matrix="&lt;&lt;T&lt;&lt;endl&lt;&lt;endl;
      cout&lt;&lt;"Rotational Matrix="&lt;&lt;R&lt;&lt;endl&lt;&lt;endl;
      cout&lt;&lt;"Camera Matrix="&lt;&lt;CameraMatrix&lt;&lt;endl;
      return 0;
    }
</code></pre>

<hr>

<p>SAMPLE 1:</p>

<pre><code>  const int numberofpoints = 30;
  float p2d[numberofpoints][2] =
    {{72, 169}, {72,  184}, {71, 264},  {96, 168},  {94, 261},  {94, 276}, {257, 133},
    {254, 322}, {254, 337}, {278, 132}, {278, 146}, {275, 321}, {439, 228},
    {439, 243}, {437, 328}, {435, 431}, {461, 226}, {459, 326}, {459, 342},
    {457, 427}, {457, 444}, {612, 196}, {609, 291}, {605, 392}, {605, 406},
    {634, 191}, {633, 206}, {630, 288}, {630, 303}, {627, 390}};

  float p3d[numberofpoints][3] =
    {{0, 0, 98},  {0, 0, 94},   {0, 0, 75},   {4, 4, 98}, {4, 4, 75}, {4, 4, 71},
    {0, 96, 75},  {0, 96, 25},  {0, 96, 21},  {4, 100, 75},
    {4, 100, 71}, {4, 100, 25}, {96, 0, 98},  {96, 0, 94},
    {96, 0, 75},  {96, 0, 50},  {100, 4, 98}, {100, 4, 75},
    {100, 4, 71}, {100, 4, 50}, {100, 4, 46}, {96, 96, 75},
    {96, 96, 50}, {96, 96, 25}, {96, 96, 21}, {100, 100, 75},
    {100, 100, 71}, {100, 100, 50}, {100, 100, 46}, {100, 100, 25}};
</code></pre>

<hr>

<p>SAMPLE 2:</p>

<pre><code> const int numberofpoints = 24;
  float p2d[24][2] =
  {{41, 291}, {41, 494},  {41, 509}, {64, 285}, {64, 303},
     {64, 491}, {195, 146},  {216, 144}, {216, 160}, {431, 337},
     {431, 441}, {431, 543}, {431, 558}, {452, 336}, {452, 349},
     {452, 438}, {452, 457}, {452, 539}, {568, 195}, {566, 291},
     {588, 190}, {587, 209}, {586, 289}, {586, 302}};
  float p3d[24][3] =
  {{0, 0, 75}, {0, 0, 25}, {0, 0, 21}, {4, 4, 75}, {4, 4, 71}, {4, 4, 25},
    {0, 96, 75 }, {4, 100, 75}, {4, 100, 71}, {96, 0, 75}, {96, 0, 50}, {96, 0, 25},
  {96, 0, 21}, {100, 4, 75}, {100, 4, 71}, {100, 4, 50}, {100, 4, 46}, {100, 4, 25}, {96,  96, 75 },  {96,  96, 50 }, {100, 100, 75}, {100, 100, 71}, {100, 100, 50}, {100, 100, 46}};
</code></pre>

<hr>

<p>SAMPLE 3:</p>

<pre><code> const int numberofpoints = 33;
  float p2d[numberofpoints][2] =
  {{45, 310},   {43, 412}, {41, 513}, {41, 530}, {68, 305},
    {68, 321}, {66, 405}, {66, 423}, {64, 509}, {201, 70}, {201, 84},
    {200, 155}, {198, 257}, {218, 259}, {218, 271}, {441, 364},
    {439, 464}, {437, 569}, {437, 586}, {462, 361}, {462, 376},
    {460, 462}, {460, 479}, {459, 565}, {583, 117}, {583, 131},
    {580, 215}, {576, 313}, {603, 113}, {600, 214}, {599, 229},
    {596, 311}, {596, 323}};

    float p3d[numberofpoints][3] =
   {{0, 0, 75},   {0, 0, 50}, {0, 0, 25}, {0, 0, 21}, {4, 4, 75}, {4, 4, 71},
    {4, 4, 50}, {4, 4, 46}, {4, 4, 25}, {0, 96, 98}, {0, 96, 94}, {0, 96, 75},
    {0, 96, 50}, {4, 100,75}, {4, 100,71}, {96, 0, 75}, {96, 0, 50}, {96, 0, 25},
    {96, 0, 21}, {100, 4,75}, {100, 4,71}, {100, 4,50}, {100, 4,46}, {100, 4,25},{96, 96,98}, {96, 96,94}, {96, 96,75}, {96, 96,50}, {100, 100, 98}, {100, 100, 75},{100, 100, 71}, {100, 100, 50}, {100, 100, 46}};
</code></pre>

<hr>