#include <vector>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <string>

#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;


Scalar getMSSIM(const Mat& i1, const Mat& i2)
{
    const double C1 = 6.5025, C2 = 58.5225;
    int d = CV_32F;

    Mat I1, I2;
    i1.convertTo(I1, d);
    i2.convertTo(I2, d);

    Mat I2_2 = I2.mul(I2);
    Mat I1_2 = I1.mul(I1);
    Mat I1_I2 = I1.mul(I2);


    Mat mu1, mu2;
    GaussianBlur(I1, mu1, Size(11, 11), 1.5);
    GaussianBlur(I2, mu2, Size(11, 11), 1.5);

    Mat mu1_2 = mu1.mul(mu1);
    Mat mu2_2 = mu2.mul(mu2);
    Mat mu1_mu2 = mu1.mul(mu2);

    Mat sigma1_2, sigma2_2, sigma12;

    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= mu1_2;

    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= mu2_2;

    GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
    sigma12 -= mu1_mu2;

    Mat t1, t2, t3;

    t1 = 2 * mu1_mu2 + C1;
    t2 = 2 * sigma12 + C2;
    t3 = t1.mul(t2);

    t1 = mu1_2 + mu2_2 + C1;
    t2 = sigma1_2 + sigma2_2 + C2;
    t1 = t1.mul(t2);

    Mat ssim_map;
    divide(t3, t1, ssim_map);

    Scalar mssim = mean(ssim_map);
    return mssim;
}


double calcMatch(const Mat& lhs, const Mat& rhs) {
    int count = 0;
    for (int i = 0; i < lhs.rows; ++i) {
        for (int j = 0; j < lhs.cols; ++j) {
            if (lhs.at<uchar>(i, j) != rhs.at<uchar>(i, j)) {
                ++count;
            }
        }
    }
    return 1.0 - count * 1.0 / (lhs.size().area());
}

const size_t img_width = 2100;
const size_t img_height = 2970;


int main() {
    vector<string> files = { "P_20210309_064128_vHDR_Auto.jpg",
                             "P_20210410_170858_vHDR_On.jpg" ,
                             "P_20210410_170941_vHDR_On.jpg" ,
                             "P_20210410_171129_vHDR_On.jpg" ,
                             "P_20210410_171136_vHDR_On.jpg" ,
                             "P_20210410_171334_vHDR_On.jpg"
    };

    int k = files.size();

    //считывание изображений
    vector<Mat> src(k);
    for (int i = 0; i < k; ++i) {
        src[i] = imread(files[i]);
    }

    //подготовка изображений
    vector<Mat> prepared(k);
    for (int i = 0; i < k; ++i) {
        Mat gr, bl, di, bin;
        cvtColor(src[i], gr, COLOR_BGR2GRAY);
        GaussianBlur(gr, bl, { 3, 3 }, 0);
        dilate(bl, di, getStructuringElement(cv::MORPH_RECT, cv::Point(9, 9)));
        adaptiveThreshold(di, bin, 255, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 51, 5);
        medianBlur(bin, prepared[i], 3);
    }

    //выделение границ и линий
    for (int i = 0; i < k; ++i) {
        Mat edges;
        Canny(prepared[i], edges, 0, 300);

        vector<Vec4i> lines;
        HoughLinesP(edges, lines, 1, CV_PI / 180, 1, 1, 20);
        for (auto& l : lines) {
            line(edges, Point(l[0], l[1]), Point(l[2], l[3]), cv::Scalar(255, 0, 0), 3, 8);
        }
        edges.copyTo(prepared[i]);
    }

    //определение контуров документа
    vector<Mat> result(k);
    vector<vector<Point>> points(k);

    for (int i = 0; i < k; ++i) {
        medianBlur(prepared[i], prepared[i], 9);
        bitwise_not(prepared[i], prepared[i]);

        vector<vector<Point>> contours;
        vector<Vec4i> hierarchy;

        findContours(prepared[i], contours, hierarchy, RETR_LIST, CHAIN_APPROX_SIMPLE);
        sort(contours.begin(), contours.end(), [](vector<Point>& lhs, vector<Point>& rhs) {
            return fabs(contourArea(Mat(lhs))) < fabs(contourArea(Mat(rhs)));
        });
        vector<Point> mainContour = contours.back();
        double perimetr = arcLength(mainContour, true);

        Mat approx;
        approxPolyDP(mainContour, approx, 0.02 * perimetr, true);

        src[i].copyTo(result[i]);
        points[i] = approx;
        Point pt = points[i].back();
        for (int j = 0; j < points[i].size(); ++j) {
            line(result[i], pt, points[i][j], { 108, 73, 255 }, 8, LINE_8);
            pt = points[i][j];
        }
        imwrite("result" + to_string(i) + ".jpg", result[i]);
    }

    //совмещение документов
    vector<Point> frame = { {0, 0}, {0, img_height}, {img_width, img_height}, {img_width, 0} };
    vector<Mat> projections(k);

    for (int i = 0; i < k; ++i) {
        Mat perspective = findHomography(points[i], frame, RANSAC);
        warpPerspective(src[i], projections[i], perspective, Size{ img_width,img_height });
        imwrite("projection" + to_string(i) + ".jpg", projections[i]);
    }

    //бинаризация документов
    vector<Mat> binarization(k);
    for (int i = 0; i < k; ++i) {
        Mat temp;
        cvtColor(projections[i], temp, COLOR_BGR2GRAY);
        adaptiveThreshold(temp, binarization[i], 255, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 111, 15);
        imwrite("binarization" + to_string(i) + ".jpg", binarization[i]);
    }

    //оценка качества
    Mat etalon_input = imread("etalon.jpg"), gray, etalon;
    cvtColor(etalon_input, gray, COLOR_BGR2GRAY);
    adaptiveThreshold(gray, etalon, 255, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 111, 15);

    vector<Mat> difference(k);
    for (int i = 0; i < k; ++i) {
        bitwise_xor(etalon, binarization[i], difference[i]);
        imwrite("difference" + to_string(i) + ".jpg", difference[i]);

        cout << files[i] << ":" << endl;
        cout << "Match with etalon = " << calcMatch(etalon, binarization[i]) << endl;
        cout << "Structure similarity = " << getMSSIM(etalon, binarization[i])[0] << endl;
        cout << endl;
    }
    waitKey(0);
}

/*
*
### Задание
1. Реализовать выделение "границ и уголков" на фотографиях документа.
2. Реализовать детектирование наиболее правдоподобного четырехугольника изображения страницы документа на фото.
3. Используя эталонную геометрическую разметку из лабораторной 5 реализовать численную оценку качества выделения отдельных примитивов и финального результата для отдельных изображений и для набора.
4. Реализовать представление результатов выполнения лабораторной (иллюстрации, таблицы, графики и т.д.) и вставить в отчет.
*/