#include <opencv/cv.hpp>
#include <opencv2/core.hpp>
#include <opencv2/calib3d.hpp>
#include "transformation.hpp"
//https://www.cnblogs.com/shengguang/p/5932522.html
void HouseHolderQR(const cv::Mat &A, cv::Mat &Q, cv::Mat &R)
{
assert ( A.channels() == 1 );
assert ( A.rows >= A.cols );
auto sign = [](double value) { return value >= 0 ? 1: -1; };
const auto totalRows = A.rows;
const auto totalCols = A.cols;
R = A.clone();
Q = cv::Mat::eye ( totalRows, totalRows, A.type() );
for ( int col = 0; col < A.cols; ++ col )
{
cv::Mat matAROI = cv::Mat ( R, cv::Range ( col, totalRows ), cv::Range ( col, totalCols ) );
cv::Mat y = matAROI.col ( 0 );
auto yNorm = norm ( y );
cv::Mat e1 = cv::Mat::eye ( y.rows, 1, A.type() );
cv::Mat w = y + sign(y.at<double>(0,0)) * yNorm * e1;
cv::Mat v = w / norm( w );
cv::Mat vT; cv::transpose(v, vT );
cv::Mat I = cv::Mat::eye( matAROI.rows, matAROI.rows, A.type() );
cv::Mat I_2VVT = I - 2 * v * vT;
cv::Mat matH = cv::Mat::eye ( totalRows, totalRows, A.type() );
cv::Mat matHROI = cv::Mat(matH, cv::Range ( col, totalRows ), cv::Range ( col, totalRows ) );
I_2VVT.copyTo ( matHROI );
R = matH * R;
Q = Q * matH;
}
}
int main()
{
/*
Thera are two notations for Projection Matrix
1)P=K[R|t]
┌X┐
┌R,t┐ |Y|
X_Cam=R*X+t ==> Homogeneous Coordinate ==> X_Cam=|0,1| |Z|
└ ┘4x4└1┘4x1
┌X┐
┌ ┐ |Y|
image plane <--x=K[I|0]3x4 |R,t| |Z|
└0,1┘4x4 └1┘4x1
P=K[R|t]
2)P=KR[I|-X0]
┌X┐
┌R -Rc┐|Y|
X_Cam=R*[X_w-C] ==> Homogeneous Coordinate ==> | ||Z|
└0 1 ┘└1┘
┌X┐
┌R -Rc┐|Y|
image plane <--x=K[I|0]3x4 |0 1 ||Z|
└ ┘└1┘
P=K[R|-RC]= KR[I|-C]
(1) and (2) ==> t=-RC
*/
int numberOfPixelInHeight,numberOfPixelInWidth;
double heightOfSensor, widthOfSensor;
double focalLength=1.50;
double mx, my, U0, V0;
numberOfPixelInHeight=600;
numberOfPixelInWidth=800;
heightOfSensor=10;
widthOfSensor=10;
my=(numberOfPixelInHeight)/heightOfSensor ;
U0=(numberOfPixelInHeight)/2 ;
mx=(numberOfPixelInWidth)/widthOfSensor;
V0=(numberOfPixelInWidth)/2;
cv::Mat cameraMatrix = (cv::Mat_<double>(3,3) <<
focalLength*mx, 0, V0,
0,focalLength*my,U0,
0,0,1);
double tx,ty,tz,roll,pitch,yaw;
tx=1.0;
ty=2.1;
tz=-1.4;
cv::Mat translation = (cv::Mat_<double>(3,1) <<tx,ty,tz);
cv::Vec3d theta;
roll=+M_PI/4 ;
pitch=+M_PI/10;
yaw=-+M_PI/6;
theta[0]=roll;
theta[1]=pitch;
theta[2]=yaw;
cv::Mat rotation=eulerAnglesToRotationMatrix(theta);
//1)P=K[R|t]
cv::Mat R_T;
cv::hconcat(rotation, translation, R_T);
cv::Mat projectionMatrix=cameraMatrix*R_T;
cv::Mat calculatedCameraMatrix,calculatedRotation,calculatedTranslation;
cv::decomposeProjectionMatrix(projectionMatrix,calculatedCameraMatrix,calculatedRotation,calculatedTranslation);
std::cout<<"==============================Ground Truth==============================" <<std::endl;
std::cout<<"Rotation Matrix (Ground Truth)" <<std::endl;
std::cout<<rotation <<std::endl;
std::cout<<"Translation Matrix (Ground Truth)" <<std::endl;
std::cout<<translation <<std::endl;
std::cout<<"Camera Matrix (Ground Truth)" <<std::endl;
std::cout<<cameraMatrix <<std::endl;
std::cout<<"==============================Decomposing Using OpenCV==============================" <<std::endl;
std::cout<<"Computed Rotation Matrix (OpenCV)" <<std::endl;
std::cout<<calculatedRotation <<std::endl;
std::cout<<"Computed Translation Matrix (OpenCV)" <<std::endl;
//std::cout<<calculatedTranslation/calculatedTranslation.at<double>(3,0) <<std::endl;
cv::Mat tempT=(cv::Mat_<double>(3,1)<<
calculatedTranslation.at<double>(0,0)/calculatedTranslation.at<double>(3,0),
calculatedTranslation.at<double>(1,0)/calculatedTranslation.at<double>(3,0),
calculatedTranslation.at<double>(2,0)/calculatedTranslation.at<double>(3,0));
std::cout<<-rotation*tempT<<std::endl;
std::cout<<"Computed Camera Matrix (OpenCV)" <<std::endl;
std::cout<<calculatedCameraMatrix <<std::endl;
///////////////////////////////////////decomposing the projection matrix///////////////////////////////////////////
/*
P=KR[I|-X0]=[H_inf3x3|h3x1]
KR=H_inf3x3
1)X0
-KRX0=h3x1 => X0=-(KR)^-1*h3x1 ==>X0=-(H_inf3x3)^-1*h3x1
2)K,R
KR=H_inf3x3 =>(KR)^-1= H_inf3x3^-1 =>R^-1*K^-1=H_inf3x3^-1 | R^-1*K^-1=Q*R => R=Q^-1, K=R^-1
H_inf3x3^-1=Q*R |
*/
cv::Mat H_inf3x3=(cv::Mat_<double>(3,3)<< projectionMatrix.at<double>(0,0),projectionMatrix.at<double>(0,1),projectionMatrix.at<double>(0,2)
,projectionMatrix.at<double>(1,0),projectionMatrix.at<double>(1,1),projectionMatrix.at<double>(1,2)
,projectionMatrix.at<double>(2,0),projectionMatrix.at<double>(2,1),projectionMatrix.at<double>(2,2));
cv::Mat h3x1=(cv::Mat_<double>(3,1)<<projectionMatrix.at<double>(0,3),projectionMatrix.at<double>(1,3),projectionMatrix.at<double>(2,3) );
cv::Mat Q,R;
cv::Mat H_inf3x3_inv=H_inf3x3.inv();
//R=Q^-1, K=R^-1
HouseHolderQR(H_inf3x3_inv, Q, R);
cv::Mat K=R.inv();
std::cout<<"==============================Decomposing Using My Code==============================" <<std::endl;
//due to homogeneity we divide it by last element
std::cout<<"Estimated Camera Matrix\n"<<K/K.at<double>(2,2) <<std::endl;
cv::Mat rotationMatrix=Q.inv();
std::cout<<"Estimated Camera Rotation\n"<<rotationMatrix*-1 <<std::endl;
std::cout<<"Estimated Camera Translation" <<std::endl;
//t=-R*C, Q.inv()=R
std::cout<<-1*(-Q.inv()*(-H_inf3x3.inv()*h3x1 ))<<std::endl;