#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

// 计算平均值
double mean(vector<double> *v) {
    double sum = 0;
    for (int i = 0; i < (*v).size(); i++) {
        sum += (*v)[i];
    }
    return sum / (*v).size();
}

// 计算协方差
double covariance(vector<double> *x, vector<double> *y) {
    double x_mean = mean(x);
    double y_mean = mean(y);
    double sum = 0;
    for (int i = 0; i < (*x).size(); i++) {
        sum += ((*x)[i] - x_mean) * ((*y)[i] - y_mean);
    }
    return sum / (*x).size();
}

// 计算相关系数
double correlation(vector<double> *x, vector<double> *y) {
    double cov = covariance(x, y);
    double x_std = sqrt(covariance(x, x));
    double y_std = sqrt(covariance(y, y));
    double res = 0;
    if (x_std != 0 && y_std != 0)
        res = cov / (x_std * y_std);
    return res;
}

//二维向量转成一维向量
vector<double> vector2to1(vector<vector<double>> *a) {
    int n = (*a).size(); // a的行数
    int m = (*a)[0].size(); // a的列数
    vector<double> b;

    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            b.push_back((*a)[i][j]);
    return b;
}
vector<vector<double>> sub_x;
//截取二维向量指定的区域
vector<vector<double>>* vector_sub(vector<vector<double>> *a, int start_row, int start_col, int end_row, int end_col) {
    sub_x.clear();
    for (int i = start_row; i < end_row; i++) {
        vector<double> row;
        for (int j = start_col; j < end_col; j++) {
            row.push_back((*a)[i][j]);
        }
        sub_x.push_back(row);
    }
    return &sub_x;
}

// 计算两个矩阵的归一化互相关系数
vector<vector<double>> correlation2(vector<vector<double>> x, vector<vector<double>> y) {
    int N1 = x.size(); // x的行数
    int N2 = x[0].size(); // x的列数
    int M1 = y.size(); // y的行数
    int M2 = y[0].size(); // y的列数
    int N = N1 + M1 - 1; // z的行数
    int M = N2 + M2 - 1; // z的列数
    vector<vector<double>> z(N, vector<double>(M, 0)); // 初始化z

    //补0的行
    for (int i = 0; i < N1 - 1; i++) {
        y.insert(y.begin(), vector<double>(M2, 0));//在y开始插入1行10元素
        y.push_back(vector<double>(M2, 0));//在y尾部追加1行0元素
    }
    //补0的列
    for (int i = 0; i < N2 - 1; i++) {
        for (int i = 0; i < y.size(); i++) {
            y[i].insert(y[i].begin(), 0);//在头部补0
            y[i].push_back(0);//在尾部补0
        }
    }
    vector<double> x1 = vector2to1(&x);
    vector<double> y1;
    for (int n = 0; n < N; n++)
        for (int m = 0; m < M; m++) {
            y1 = vector2to1(vector_sub(&y, n, m, n + N1, m + N2));
            z[n][m] = correlation(&x1, &y1);
        }

    return z;
}

int main() {

    vector<vector<double>> img1;
    Mat image1 = imread("C:/Users/MingYi-LZQ/Desktop/3.tif", IMREAD_UNCHANGED);
    if (image1.empty()) {
        cout << "Could not open or find the image" << endl;
        return -1;
    }

    img1.resize(image1.rows, vector<double>(image1.cols));
    for (int i = 0; i < image1.rows; i++) {
        for (int j = 0; j < image1.cols; j++) {
            img1[i][j] = (double)image1.at<uchar>(i, j);
        }
    }

    vector<vector<double>> img2;
    Mat image2 = imread("C:/Users/MingYi-LZQ/Desktop/4.tif", IMREAD_UNCHANGED);
    if (image2.empty()) {
        cout << "Could not open or find the image" << endl;
        return -1;
    }

    img2.resize(image2.rows, vector<double>(image2.cols));
    for (int i = 0; i < image2.rows; i++) {
        for (int j = 0; j < image2.cols; j++) {
            img2[i][j] = (double)image2.at<uchar>(i, j);
        }
    }
    vector<vector<double>> z = correlation2(img1, img2);

    double max_corr = 0;
    int x_peak = 0;
    int y_peak = 0;
    vector<vector<double>> offset;

    for (int i = 0; i < z.size(); i++) {
        for (int j = 0; j < z[0].size(); j++) {
            if (z[i][j] > max_corr) {
                max_corr = z[i][j];//找到最大的相关系数
                x_peak = j + 1;//最大相关系数对应的x坐标
                y_peak = i + 1;//最大相关系数对应的y坐标
            }
        }
    }
    //计算出实际x y的偏移量
    vector<double> temp;
    temp.push_back(x_peak - img1.size());
    temp.push_back(y_peak - img1[0].size());
    offset.push_back(temp);
    return 0;
}