基於OpenCV影象最小二乘相位解包裹演算法
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(); }