**

1.演算法原理

**

《數字正射影像鑲嵌中色彩一致性處理的若干問題研究》

**

2.程式碼實現

**

#include <opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>  
#include<opencv2/imgproc/imgproc.hpp>  
#include <iostream>  

using namespace cv;
using namespace std;

Mat getFFTresultImg(Mat& completeI,CvSize srcSize)
{
    Mat planes[2
 
];
    split(completeI, planes);//把變換後的結果分割到各個陣列的兩頁中，方便後續操作 
    Mat magI;
    magnitude(planes[0], planes[1], magI);//求傅立葉變換各頻率的幅值，幅值放在第一頁中。  


    //傅立葉變換的幅度值範圍大到不適合在螢幕上顯示。高值在螢幕上顯示為白點，  
    //而低值為黑點，高低值的變化無法有效分辨。為了在螢幕上凸顯出高低變化的連續性，我們可以用對數尺度來替換線性尺度:  
    magI += 1;
    log(magI, magI);//取對數  
    magI = magI(Rect(0
 
, 0,srcSize.width,srcSize.height));//前邊對原始影象進行了擴充套件，這裡把對原始影象傅立葉變換取出，剔除擴充套件部分。  


    //這一步的目的仍然是為了顯示。 現在我們有了重分佈後的幅度圖，  
    //但是幅度值仍然超過可顯示範圍[0,1] 。我們使用 normalize() 函式將幅度歸一化到可顯示範圍。  
    normalize(magI, magI, 0, 1, CV_MINMAX);//傅立葉影象進行歸一化。  


    //重新分配象限，使（0,0）移動到影象中心，  
    //在《數字影象處理》中，傅立葉變換之前要對源影象乘以（-1）^(x+y)進行中心化。  
 

    //這是是對傅立葉變換結果進行中心化  
    int cx = magI.cols / 2;
    int cy = magI.rows / 2;

    Mat tmp;
    Mat q0(magI, Rect(0, 0, cx, cy));
    Mat q1(magI, Rect(cx, 0, cx, cy));
    Mat q2(magI, Rect(0, cy, cx, cy));
    Mat q3(magI, Rect(cx, cy, cx, cy));


    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
    return magI;
}

Mat FFT(Mat& src_gray)
{
    //Mat src_gray;
    //cvtColor(src, src_gray, CV_RGB2GRAY);//灰度影象做傅立葉變換  

    int m = getOptimalDFTSize(src_gray.rows);//2,3,5的倍數有更高效率的傅立葉轉換  
    int n = getOptimalDFTSize(src_gray.cols);

    Mat dst;
    ///把灰度影象放在左上角，在右邊和下邊擴充套件影象，擴充套件部分填充為0；  
    // 0, m - src_gray.rows, 0, n - src_gray.cols 上邊填充0行，下面填充m - src_gray.rows行
    copyMakeBorder(src_gray, dst, 0, m - src_gray.rows, 0, n - src_gray.cols, BORDER_CONSTANT, Scalar::all(0));
    //cout << dst.size() << endl;

    //新建一個兩頁的array，其中第一頁用擴充套件後的影象初始化，第二頁初始化為0  
    Mat planes[] = { Mat_<float>(dst), Mat::zeros(dst.size(), CV_32F) };
    Mat  completeI;
    merge(planes, 2, completeI);//把兩頁合成一個2通道的mat  

    //對上邊合成的mat進行傅立葉變換，支援原地操作，傅立葉變換結果為複數。通道1存的是實部，通道2存的是虛部。  
    dft(completeI, completeI);


    return completeI.clone();
}
//計算高斯濾波係數矩陣
Mat clcGLPFMat( Mat& mat,int D0)
{
    int width = mat.rows;
    int height = mat.cols;
    int M = width;
    int N = height;
    Mat mat_GLPF(mat.size(), CV_32FC1);

    Mat U, V;
    U.create(M, N, CV_32FC1);
    V.create(M, N, CV_32FC1);

    for (int u = 0; u < M; ++u)
    {
        for (int v = 0; v < N; ++v)
        {
            float tm1, tm2;
            tm1 = (float)((u > cvRound(M / 2)) ? u - M : u);
            tm2 = (float)((v > cvRound(N / 2)) ? v - N : v);

            U.at<float>(u, v) = tm1;
            V.at<float>(u, v) = tm2;
        }
    }


    for (int u = 0; u < M; ++u)
    {
        for (int v = 0; v < N; ++v)
        {
            float t1, t2;
            t1 = U.at<float>(u, v);
            t2 = V.at<float>(u, v);
            float Elem_D = t1*t1 + t2*t2;
            mat_GLPF.at<float>(u, v) = (float)(exp(-(Elem_D) / (2 * D0*D0)) / 2 / 3.1415 / (2 * D0*D0));
        }
    }
    Mat_<float>::iterator begainIt = mat_GLPF.begin<float>();
    Mat_<float>::iterator endIt = mat_GLPF.end<float>();
    float sumValue = 0;
    for (; begainIt != endIt; begainIt++)
    {
        sumValue += *begainIt;
    }
    mat_GLPF = mat_GLPF / sumValue;
    return mat_GLPF.clone();
}

Mat mask(Mat& plane)
{
    Mat FFTresult = FFT(plane);//傅立葉變換包含實部和虛部，分別放在兩個planes裡
    Mat planes[2];
    split(FFTresult, planes);
    imshow("FFTresult",getFFTresultImg(FFTresult, FFTresult.size()));

    Mat GLPFMatIM = clcGLPFMat(planes[0], 10);//高斯濾波係數
    Mat GLPFMatRE = clcGLPFMat(planes[1],10);
    planes[0] = GLPFMatIM.mul(planes[0]);
    planes[1] = GLPFMatRE.mul(planes[1]);
    Mat GLPFresult;
    merge(planes, 2, GLPFresult);       //實部虛部分別高斯濾波，然後合成到濾波結果
    imshow("FFTresultAfterFlit", getFFTresultImg(GLPFresult, GLPFresult.size()));

    Mat maskResult;
    dft(GLPFresult, maskResult, DFT_INVERSE + DFT_SCALE);//濾波結果做傅立葉反變換

    split(maskResult, planes);//把反變換後的結果分割到兩頁中，方便後續操作 
    Mat mask;
    magnitude(planes[0], planes[1], mask);//求傅立葉變換各頻率的幅值  
    return mask.clone();
}

int main(){
    Mat src = imread("g:\\Opencv\\lb.jpg");
    if (src.empty())
    {
        return-1;
    }
    Mat BGR[3];
    split(src,BGR);
    BGR[0] = mask(BGR[0]);
    BGR[1] = mask(BGR[1]);
    BGR[2] = mask(BGR[2]);
    Mat maskBGR;
    merge(BGR, 3, maskBGR);
    maskBGR = maskBGR(Rect(0, 0, src.size().width, src.size().height));
    cout << maskBGR.size()<<maskBGR.type();
    cout << src.size()<<src.type();
    src.convertTo(src, maskBGR.type(), 1 / 255.0);
    imshow("maskBGR", maskBGR);
    maskBGR = src - maskBGR+20.0/255.0;
    namedWindow("InputImage");
    imshow("InputImage", src);
    imshow("RESULT", maskBGR);

    waitKey();
    return 0;
}

**

3.實驗結果

**

高斯濾波前後的頻譜圖

原始影象、還原影象g(x,y)、結果影象
結果影象中沒有加灰度偏移量offest，影象亮度偏低