版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/qq_35759050/article/details/53819667
QT:5.5.1 编译器:mcvs2012 OpenCV:2.4.9
算法写在Qt上一个按钮上,相当于一个main函数,主要是变化的源图像与目标图像的深度,通道的匹配容易报错。
void MainWindow::on_pushButton_clicked()
{
double pi = 3.14156;
IplImage *img = cvLoadImage("D:/3.jpg"); //读取本地图片
IplImage *imgGray = cvCreateImage(cvGetSize(img),IPL_DEPTH_8U,1); //创建灰度图头
cvCvtColor(img,imgGray,CV_RGB2GRAY); //转为灰度图
IplImage *imgSobel = cvCreateImage(cvGetSize(img),IPL_DEPTH_8U,1); //图像在x,y方向上做两次微分
cvSobel(imgGray,imgSobel,2,2,5);
IplImage *imgSobel_64F = cvCreateImage(cvGetSize(img),IPL_DEPTH_64F,1); //矩阵深度由8U转为64F
cvConvertScale(imgSobel, imgSobel_64F, 1.0/255, 0);
IplImage *imgPhase = cvCreateImage(cvGetSize(img),IPL_DEPTH_64F,1); //做正方向离散余弦变化
cvDCT(imgSobel_64F,imgPhase,CV_DXT_FORWARD);
for(int i=0;i<imgPhase->height;i++) //遍历图像
for(int j=0;j<imgPhase->width;j++)
{
double k1 = 2*cos((i-1)*pi/(imgPhase->height));
double k2 = 2*cos((j-1)*pi/(imgPhase->width));
double k = k1 + k2 - 4;
CvScalar k0 = cvScalar(cvGet2D(imgPhase,i,j).val[0]/k);
cvSet2D(imgPhase,i,j,k0); //对不同元素做相应加减
}
IplImage *imgPhaseDst = cvCreateImage(cvGetSize(imgPhase),IPL_DEPTH_64F,1);
CvScalar tmp = cvScalar((cvGet2D(imgPhase,1,1).val[0]-cvGet2D(imgPhase,1,2).val[0]-cvGet2D(imgPhase,2,1).val[0])/2); //第一个元素特殊处理
cvSet2D(imgPhaseDst,1,1,tmp);
cvDCT(imgPhase,imgPhaseDst,CV_DXT_INVERSE); //反离散余弦变化
cvNamedWindow("Phase"); //显示解包裹相位图
cvShowImage("Phase",imgPhaseDst);
}
下面是MATLAB代码,写在子函数中:
%% ***********************横向剪切最小二乘解包裹算法函数*************
function Un_Pha = Unwrap_Alg(Phi)
%**************判断输入图片转化为单通道(灰度图)***********
X = size(Phi);
Phi = double(rgb2gray(Phi));
%********************对包裹相位求一阶偏微分**************
[m,n] = size(Phi);
Phidx = zeros(m,n);
phidy = zeros(m,n);
Phidx(1:m-1,:) = angle(exp(j*(Phi(2:m,:)-Phi(1:m-1,:)))); %在X方向微分
phidy(:,1:n-1) = angle(exp(j*(Phi(:,2:n)-Phi(:,1:n-1)))); %在Y方向微分
%********************对包裹相位求二阶偏微分**************
Rou3 = zeros(m,n);
Rou3dx = zeros(m,n);
Rou3dy = zeros(m,n);
Rou3dx(1:m-1,:) = angle(exp(j*(Phidx(2:m,:)-Phidx(1:m-1,:))));%在X方向二次微分
Rou3dy(:,1:n-1) = angle(exp(j*(phidy(:,2:n)-phidy(:,1:n-1))));%在Y方向二次微分
Rou3 = Rou3dx + Rou3dy; %二阶的二次微分值相加
%***********************DCT求解泊松方程********************
% 横向剪切最小二乘解包裹实际是求解离散泊松方程
tic %时间计时开始
PP3 = dct2(Rou3); %离散泊松变换
for ii=1:m
for jj=1:n
k1=2*cos((ii-1)*pi/(m));
k2=2*cos((jj-1)*pi/(n));
KK = k1+k2-4;
PH3(ii,jj) = PP3(ii,jj)/KK; %离散泊松变换后,乘以一个因子
end
end
PH3(1,1) = -(PH3(1,2) + PH3(2,1) - PP3(1,1))/2;
Un_Pha = idct2(PH3); %离散泊松反变换
toc %时间计时结束
Un_Pha = Un_Pha(1:m,1:n);
Mat类型的,图片的行数列数必须是偶数:
void MainWindow::on_pushButton_3_clicked()
{
double pi = 3.14156;
Mat img = imread("D:/1.jpg",1);
Mat imgGray;
cvtColor(img,imgGray,CV_BGR2GRAY);
Mat imgSobel;
Sobel(imgGray,imgSobel,CV_8U,2,2,3,1,0,BORDER_DEFAULT);
Mat imgPhase,imgPhaseTmp;
imgSobel.convertTo(imgPhaseTmp,CV_64FC1);
cv::dct(imgPhaseTmp,imgPhase,CV_DXT_FORWARD);
for(int i=1;i<(imgPhase.rows);i++)
for(int j=1;j<(imgPhase.cols);j++)
{
double k1 = 2*cos((i-1)*pi/(imgPhase.rows));
double k2 = 2*cos((j-1)*pi/(imgPhase.cols));
double k = k1 + k2 -4;
imgPhase.at<uchar>(i,j) /= k;
}
imgPhase.at<uchar>(1,1) = (imgPhase.at<uchar>(1,1)-imgPhase.at<uchar>(1,2)-imgPhase.at<uchar>(2,1))/2;
Mat imgPhaseDst;
cv::dct(imgPhase,imgPhaseDst,CV_DXT_INVERSE);
namedWindow("src",1);
imshow("src",imgPhaseDst);
}
以上用opencv的代码有一些考虑不周到的地方,改进一下,还是有些问题,矩阵存在负数问题,利用Mat_加以解决,DCT离散余弦变换的问题,当图像矩阵是double型,dct正变换,再变换,图形只剩下边界,原因不知,正在解决。
double pi = 3.141593;
Mat imgGray = imread("D:/Pic/4.jpg",0);
//imshow("gray",imgGray);
//Mat imgGray;
//cvtColor(img,imgGray,CV_BGR2GRAY);
Mat imgGray1,imgGray2,imgSobelX1,imgSobelY1,imgSobelX1Tmp,imgSobelY1Tmp,imgSobelX2,imgSobelY2,imgSobelXY;
imgGray1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgGray2.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelX1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelY1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelX1Tmp.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelY1Tmp.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelX2.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelY2.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelXY.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
double xSin,yCos;
for(int i=0;i<imgGray1.rows;i++)
for(int j=0;j<imgGray1.cols;j++)
{
imgGray1.at<uchar>(i,j) = imgGray.at<uchar>(i+1,j)-imgGray.at<uchar>(i,j);
imgGray2.at<uchar>(i,j) = imgGray.at<uchar>(i,j+1)-imgGray.at<uchar>(i,j);
}
//imshow("imgGray1",imgGray1);
for(int i=0;i<imgGray1.rows;i++)
for(int j=0;j<imgGray1.cols;j++)
{
xSin = sin(imgGray1.at<uchar>(i,j));
yCos = cos(imgGray1.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX1.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX1.at<uchar>(i,j) = 0 - pi/2;
}
else
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = pi/2;
else
imgSobelX1.at<uchar>(i,j) = 0;
}
//imshow("x",imgSobelX1);
for(int i=0;i<imgGray2.rows;i++)
for(int j=0;j<imgGray2.cols;j++)
{
xSin = sin(imgGray2.at<uchar>(i,j));
yCos = cos(imgGray2.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY1.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY1.at<uchar>(i,j) = 0-pi/2;
}
else
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = pi/2;
else
imgSobelY1.at<uchar>(i,j) = 0;
}
//imshow("y",imgSobelY1);
////////////////////////2次微分//////////////////////////////
for(int i=0;i<imgSobelX1Tmp.rows;i++)
for(int j=0;j<imgSobelX1Tmp.cols;j++)
{
imgSobelX1Tmp.at<uchar>(i,j) = imgSobelX1.at<uchar>(i+1,j)-imgSobelX1.at<uchar>(i,j);
imgSobelY1Tmp.at<uchar>(i,j) = imgSobelY1.at<uchar>(i,j+1)-imgSobelY1.at<uchar>(i,j);
}
for(int i=0;i<imgSobelX1Tmp.rows;i++)
for(int j=0;j<imgSobelX1Tmp.cols;j++)
{
xSin = sin(imgSobelX1Tmp.at<uchar>(i,j));
yCos = cos(imgSobelX1Tmp.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX2.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = -atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX2.at<uchar>(i,j) = -pi/2;
}
else
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = pi/2;
else
imgSobelX2.at<uchar>(i,j) = 0;
}
for(int i=0;i<imgSobelY1Tmp.rows;i++)
for(int j=0;j<imgSobelY1Tmp.cols;j++)
{
xSin = sin(imgSobelY1Tmp.at<uchar>(i,j));
yCos = cos(imgSobelY1Tmp.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY2.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY2.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY2.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY2.at<uchar>(i,j) = 0-pi/2;
}
else
if(xSin>0)
imgSobelY2.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = pi/2;
else
imgSobelY2.at<uchar>(i,j) = 0;
}
imgSobelXY = imgSobelX2 + imgSobelY2;
//imshow("xy",imgSobelXY);
Mat imgPhase,imgPhaseTmp,imgPhaseTmp1;
imgSobelXY.convertTo(imgPhaseTmp,CV_64FC1);
cv::dct(imgPhaseTmp,imgPhaseTmp1,CV_DXT_FORWARD);
imgPhase.create(imgPhaseTmp.rows,imgPhaseTmp.cols,imgPhaseTmp.type());
//imshow("imgPhase",imgPhaseTmp1);
for(int i=0;i<(imgPhaseTmp1.rows);i++)
for(int j=0;j<(imgPhaseTmp1.cols);j++)
{
double k1 = 2*cos(i*pi/(imgPhaseTmp1.rows));
double k2 = 2*cos(j*pi/(imgPhaseTmp1.cols));
double k = k1 + k2 - 4; //k为负数
//qDebug()<<k;
imgPhase.at<uchar>(i,j) = imgPhaseTmp1.at<uchar>(i,j) / k;
//imgPhase.at<uchar>(i,j) = k;
//qDebug()<<imgPhase.at<uchar>(i,j);
}
//imshow("imgPhase1",imgPhase);
Mat imgPhaseDst;
cv::dct(imgPhase,imgPhaseDst,CV_DXT_INVERSE); //imgPhaseDst为最后的结果图
以下是Mat_类型的,此类型矩阵可以存在负数。
double pi = 3.141593;
Mat img = imread("D:/Pic/4.jpg",1);
Mat imgGray;
cvtColor(img,imgGray,CV_BGR2GRAY);
//imwrite("D:/5.jpg",imgGray);
Mat imgGray1,imgGray2,imgSobelX1,imgSobelY1,imgSobelX1Tmp,imgSobelY1Tmp,imgSobelX2,imgSobelY2,imgSobelXY;
imgGray1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgGray2.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelX1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelY1.create(imgGray.rows-1,imgGray.cols-1,imgGray.type());
imgSobelX1Tmp.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelY1Tmp.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelX2.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelY2.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
imgSobelXY.create(imgGray.rows-2,imgGray.cols-2,imgGray.type());
double xSin,yCos;
//Sobel(imgGray,imgSobelX1,CV_8U,1,0,7,1,0,BORDER_DEFAULT);
for(int i=0;i<imgGray1.rows;i++)
imgGray1.row(i) = imgGray.row(i+1)-imgGray.row(i);
for(int i=0;i<imgGray1.rows;i++)
for(int j=0;j<imgGray1.cols;j++)
{
xSin = sin(imgGray1.at<uchar>(i,j));
yCos = cos(imgGray1.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX1.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX1.at<uchar>(i,j) = 0-pi/2;
}
else
if(yCos>0)
imgSobelX1.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelX1.at<uchar>(i,j) = pi/2;
else
imgSobelX1.at<uchar>(i,j) = 0;
}
//Sobel(imgGray,imgSobelY1,CV_8U,0,1,7,1,0,BORDER_DEFAULT);
for(int i=0;i<imgGray2.cols;i++)
imgGray2.col(i) = imgGray.col(i+1)-imgGray.col(i);
for(int i=0;i<imgGray2.rows;i++)
for(int j=0;j<imgGray2.cols;j++)
{
xSin = sin(imgGray2.at<uchar>(i,j));
yCos = cos(imgGray2.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY1.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY1.at<uchar>(i,j) = 0-pi/2;
}
else
if(yCos>0)
imgSobelY1.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelY1.at<uchar>(i,j) = pi/2;
else
imgSobelY1.at<uchar>(i,j) = 0;
}
////////////////////////2次微分//////////////////////////////
//Sobel(imgSobelX1,imgSobelX2,CV_8U,1,0,7,1,0,BORDER_DEFAULT);
for(int i=0;i<imgSobelX1Tmp.rows;i++)
imgSobelX1Tmp.row(i) = imgSobelX1.row(i+1)-imgSobelX1.row(i);
for(int i=0;i<imgSobelX1Tmp.rows;i++)
for(int j=0;j<imgSobelX1Tmp.cols;j++)
{
xSin = sin(imgSobelX1Tmp.at<uchar>(i,j));
yCos = cos(imgSobelX1Tmp.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX2.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = -atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX2.at<uchar>(i,j) = -pi/2;
}
else
if(yCos>0)
imgSobelX2.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelX2.at<uchar>(i,j) = pi/2;
else
imgSobelX2.at<uchar>(i,j) = 0;
}
//Sobel(imgSobelY1,imgSobelY2,CV_8U,0,1,7,1,0,BORDER_DEFAULT);
for(int i=0;i<imgSobelY1Tmp.cols;i++)
imgSobelY1Tmp.col(i) = imgSobelY1.col(i+1)-imgSobelY1.col(i);
for(int i=0;i<imgSobelY1Tmp.rows;i++)
for(int j=0;j<imgSobelY1Tmp.cols;j++)
{
xSin = sin(imgSobelY1Tmp.at<uchar>(i,j));
yCos = cos(imgSobelY1Tmp.at<uchar>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY2.at<uchar>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY2.at<uchar>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY2.at<uchar>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY2.at<uchar>(i,j) = 0-pi/2;
}
else
if(xSin>0)
imgSobelY2.at<uchar>(i,j) = 0;
else if(yCos<0)
imgSobelY2.at<uchar>(i,j) = pi/2;
else
imgSobelY2.at<uchar>(i,j) = 0;
}
imgSobelXY = imgSobelX2 + imgSobelY2;
//namedWindow("x",1);
//namedWindow("y",1);
namedWindow("xy",1);
//imshow("x",imgSobelX2);
//imshow("y",imgSobelY2);
imshow("xy",imgSobelXY);
//for(int i=200;i<220;i++)
//for(int j=200;j<220;j++)
//qDebug()<<imgSobelXY.at<uchar>(i,j);
Mat imgPhase,imgPhaseTmp;
imgSobelXY.convertTo(imgPhaseTmp,CV_64FC1);
cv::dct(imgPhaseTmp,imgPhase,CV_DXT_FORWARD);
//namedWindow("phasePre",1);
//imshow("phasePre",imgPhase);
/*Mat imgPhase1;
imgPhase1.create(imgPhase.size(),imgPhase.type());
for(int i=0;i<imgPhase.rows;i++)
{
const uchar* inData=imgPhase.ptr<uchar>(i);
uchar* outData=imgPhase1.ptr<uchar>(i);
for(int j=0;j<imgPhase.cols;j++)
{
double k1 = 2*cos((i-1)*pi/(imgPhase.rows));
double k2 = 2*cos((j-1)*pi/(imgPhase.cols));
double k = k1 + k2 - 4;
outData[j]=inData[j]/k;
}
}*/
for(int i=1;i<(imgPhase.rows);i++)
for(int j=1;j<(imgPhase.cols);j++)
{
double k1 = 2*cos((i-1)*pi/(imgPhase.rows));
double k2 = 2*cos((j-1)*pi/(imgPhase.cols));
double k = k1 + k2 - 4;
//qDebug()<<k;
//qDebug()<<k2;
imgPhase.at<uchar>(i,j) /= k;
//qDebug()<<imgPhase.at<uchar>(i,j);
}
//imgPhase.at<uchar>(1,1) = (imgPhase.at<uchar>(1,1)-imgPhase.at<uchar>(1,2)-imgPhase.at<uchar>(2,1))/2;
//namedWindow("phasePre1",1);
//imshow("phasePre1",imgPhase);
Mat imgPhaseDst;
cv::dct(imgPhase,imgPhaseDst,CV_DXT_INVERSE);
namedWindow("phase",1);
imshow("phase",imgPhaseDst);
----------------------------------------------------------分割线(2017.5.9)----------------------------------------------------------------
终于在opencv环境下实现了最小二乘去包裹算法,贴上MATLAB(第一幅图)与opencv(第二幅图)的效果图:
总结一下出现的问题及解决的方法:
imshow()函数,显示矩阵中的负数时,会将负数置0显示,也就是一个黑色像素;
Mat矩阵运算时,如果出现负数,也会置0;
opencv离散余弦正变换的第一个点显示INF(无穷大),导致反变换矩阵都变成INF,用相邻点的值代替INF,总的来说,MATLAB与opencv的DCT变换效果不太一样。
处理opencv中的负数矩阵,利用Mat_型:
img = imread("D:/Pic/4.jpg",0);
Mat_<int> img1 = img;
opencv中好多自带的函数的输入矩阵需要double型,可直接用Mat_<double>型的矩阵作为输入。
为了便于观察,用qwtplot3d图形库对图像进行三维显示:
代码如下(调试注释的代码没有去掉):
void MainWindow::on_pushButton_4_clicked()
{
double pi = 3.141593;
Mat img = imread("D:/Pic/4.jpg",0);
Mat_<int> imgGray = img;
Mat_<int> imgGray1,imgGray2,imgSobelX1,imgSobelY1,imgSobelX1Tmp,imgSobelY1Tmp,imgSobelX2,imgSobelY2,imgSobelXY;
imgGray1.create(imgGray.rows-1,imgGray.cols-1);
imgGray2.create(imgGray.rows-1,imgGray.cols-1);
imgSobelX1.create(imgGray.rows-1,imgGray.cols-1);
imgSobelY1.create(imgGray.rows-1,imgGray.cols-1);
imgSobelX1Tmp.create(imgGray.rows-2,imgGray.cols-2);
imgSobelY1Tmp.create(imgGray.rows-2,imgGray.cols-2);
imgSobelX2.create(imgGray.rows-2,imgGray.cols-2);
imgSobelY2.create(imgGray.rows-2,imgGray.cols-2);
imgSobelXY.create(imgGray.rows-2,imgGray.cols-2);
float xSin,yCos;
for(int i=0;i<imgGray1.rows;i++)
for(int j=0;j<imgGray1.cols;j++)
{
imgGray1.at<int>(i,j) = imgGray.at<int>(i+1,j)-imgGray.at<int>(i,j);
imgGray2.at<int>(i,j) = imgGray.at<int>(i,j+1)-imgGray.at<int>(i,j);
}
//imshow("imgGray1",imgGray1);
for(int i=0;i<imgGray1.rows;i++)
for(int j=0;j<imgGray1.cols;j++)
{
xSin = sin(imgGray1.at<int>(i,j));
yCos = cos(imgGray1.at<int>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX1.at<int>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<int>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX1.at<int>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX1.at<int>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX1.at<int>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX1.at<int>(i,j) = 0 - pi/2;
}
else
if(yCos>0)
imgSobelX1.at<int>(i,j) = 0;
else if(yCos<0)
imgSobelX1.at<int>(i,j) = pi/2;
else
imgSobelX1.at<int>(i,j) = 0;
}
//imshow("x",imgSobelX1);
for(int i=0;i<imgGray2.rows;i++)
for(int j=0;j<imgGray2.cols;j++)
{
xSin = sin(imgGray2.at<int>(i,j));
yCos = cos(imgGray2.at<int>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY1.at<int>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<int>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY1.at<int>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY1.at<int>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY1.at<int>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY1.at<int>(i,j) = 0-pi/2;
}
else
if(yCos>0)
imgSobelY1.at<int>(i,j) = 0;
else if(yCos<0)
imgSobelY1.at<int>(i,j) = pi/2;
else
imgSobelY1.at<int>(i,j) = 0;
}
//imshow("y",imgSobelY1);
////////////////////////2次微分//////////////////////////////
for(int i=0;i<imgSobelX1Tmp.rows;i++)
for(int j=0;j<imgSobelX1Tmp.cols;j++)
{
imgSobelX1Tmp.at<int>(i,j) = imgSobelX1.at<int>(i+1,j)-imgSobelX1.at<int>(i,j);
imgSobelY1Tmp.at<int>(i,j) = imgSobelY1.at<int>(i,j+1)-imgSobelY1.at<int>(i,j);
}
for(int i=0;i<imgSobelX1Tmp.rows;i++)
for(int j=0;j<imgSobelX1Tmp.cols;j++)
{
xSin = sin(imgSobelX1Tmp.at<int>(i,j));
yCos = cos(imgSobelX1Tmp.at<int>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelX2.at<int>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<int>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelX2.at<int>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelX2.at<int>(i,j) = -atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelX2.at<int>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelX2.at<int>(i,j) = -pi/2;
}
else
if(yCos>0)
imgSobelX2.at<int>(i,j) = 0;
else if(yCos<0)
imgSobelX2.at<int>(i,j) = pi/2;
else
imgSobelX2.at<int>(i,j) = 0;
}
for(int i=0;i<imgSobelY1Tmp.rows;i++)
for(int j=0;j<imgSobelY1Tmp.cols;j++)
{
xSin = sin(imgSobelY1Tmp.at<int>(i,j));
yCos = cos(imgSobelY1Tmp.at<int>(i,j));
if(xSin>0)
{
if(yCos>0)
imgSobelY2.at<int>(i,j) = atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<int>(i,j) = pi - atan(abs(xSin/yCos));
else
imgSobelY2.at<int>(i,j) = pi/2;
}
else if(xSin<0)
{
if(yCos>0)
imgSobelY2.at<int>(i,j) = 0-atan(abs(xSin/yCos));
else if(yCos<0)
imgSobelY2.at<int>(i,j) = atan(abs(xSin/yCos)) - pi;
else
imgSobelY2.at<int>(i,j) = 0-pi/2;
}
else
if(xSin>0)
imgSobelY2.at<int>(i,j) = 0;
else if(yCos<0)
imgSobelY2.at<int>(i,j) = pi/2;
else
imgSobelY2.at<int>(i,j) = 0;
}
imgSobelXY = imgSobelX2 + imgSobelY2;
//imshow("xy",imgSobelXY);
/*for(int i=0;i<imgSobelXY.rows;i++)
for(int j=0;j<imgSobelXY.cols;j++)
{
qDebug()<<imgSobelXY.at<double>(i,j);
}*/
Mat_<double> imgSobelXYD(imgSobelXY.rows,imgSobelXY.cols,CV_64FC1);
imgSobelXY.convertTo(imgSobelXYD,CV_64FC1);
Mat_<double> imgPhase,imgPhaseTmp;
//imgSobelXY.convertTo(imgPhaseTmp,CV_64FC1);
cv::dct(imgSobelXYD,imgPhaseTmp,CV_DXT_FORWARD);
imgPhase.create(imgPhaseTmp.rows,imgPhaseTmp.cols);
//imshow("imgPhase",imgPhaseTmp);
for(int i=0;i<(imgPhaseTmp.rows);i++)
for(int j=0;j<(imgPhaseTmp.cols);j++)
{
double k1 = 2*cos(i*pi/(imgPhaseTmp.rows));
double k2 = 2*cos(j*pi/(imgPhaseTmp.cols));
double k = k1 + k2 - 4; //k为负数
//qDebug()<<k;
imgPhase.at<double>(i,j) = imgPhaseTmp.at<double>(i,j) / k;
if(imgPhase.at<double>(i,j) > 100000.0 )
imgPhase.at<double>(i,j) = imgPhase.at<double>(i,j-1);
//imgPhase.at<uchar>(i,j) = k;
//qDebug()<<imgPhase.at<uchar>(i,j);
}
//imshow("imgPhase",imgPhase);
//imgPhase.at<double>(0,0)=imgPhase.at<double>(0,1);
//for(int i=0;i<imgPhase.rows;i++)
// for(int j=0;j<imgPhase.cols;j++)
// qDebug()<<imgPhase.at<double>(i,j);
Mat_<double> imgPhaseDst;
cv::dct(imgPhase,imgPhaseDst,CV_DXT_INVERSE);
//imshow("img",imgPhaseDst);
//for(int i=0;i<imgPhaseDst.rows;i++)
//for(int j=0;j<imgPhaseDst.cols;j++)
//qDebug()<<imgPhaseDst.at<double>(i,j);
//////////////////////////////3D//////////////////////////////////
//Mat imgx(imgPhaseDst.rows, imgPhaseDst.cols, CV_64FC1);
//imgPhaseDst.convertTo(imgx, CV_64FC1); //imgx为要显示的double矩阵
Mat_<double> imgx = imgPhaseDst;
Plot* plot = new Plot(imgx);
QGridLayout *grid = new QGridLayout(ui->frame);
grid->addWidget(plot,0,0);
plot->setTitle("3D");
plot->coordinates()->axes[X1].setLabelString("X(um)"); //只能写在这
plot->coordinates()->axes[Y1].setLabelString("Y(um)"); //设置坐标轴标签
plot->coordinates()->axes[Z1].setLabelString("Z(nm)");
plot->coordinates()->axes[Z2].setLabelString("Z(nm)");
plot->coordinates()->axes[Z3].setLabelString("Z(nm)");
plot->coordinates()->axes[Z4].setLabelString("Z(nm)");
double start,stop; //设置颜色条
plot->coordinates()->axes[Z1].limits(start,stop);
plot->legend()->setLimits(start,stop);
// plot->legend()->setAutoScale(true);
// plot->resize(600,600);
plot->show();
grid->~QGridLayout();
}