#include "opencv2/core.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/ml.hpp"

#include <algorithm>
#include <iostream>
#include <vector>

using namespace cv;
using namespace std;

const int SZ = 20;  // size of each digit is SZ x SZ
const int CLASS_N = 10;
const char* DIGITS_FN = "digits.png";

static void help()
{
    cout <<
    "\n"
    "SVM and KNearest digit recognition.\n"
    "\n"
    "Sample loads a dataset of handwritten digits from 'digits.png'.\n"
    "Then it trains a SVM and KNearest classifiers on it and evaluates\n"
    "their accuracy.\n"
    "\n"
    "Following preprocessing is applied to the dataset:\n"
    " - Moment-based image deskew (see deskew())\n"
    " - Digit images are split into 4 10x10 cells and 16-bin\n"
    "   histogram of oriented gradients is computed for each\n"
    "   cell\n"
    " - Transform histograms to space with Hellinger metric (see [1] (RootSIFT))\n"
    "\n"
    "\n"
    "[1] R. Arandjelovic, A. Zisserman\n"
    "    \"Three things everyone should know to improve object retrieval\"\n"
    "    http://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/arandjelovic12.pdf\n"
    "\n"
    "Usage:\n"
    "   ./digits\n" << endl;
}

static void split2d(const Mat& image, const Size cell_size, vector<Mat>& cells)
{
    int height = image.rows;
    int width = image.cols;

    int sx = cell_size.width;
    int sy = cell_size.height;

    cells.clear();

    for (int i = 0; i < height; i += sy)
    {
        for (int j = 0; j < width; j += sx)
        {
            cells.push_back(image(Rect(j, i, sx, sy)));
        }
    }
}

static void load_digits(const char* fn, vector<Mat>& digits, vector<int>& labels)
{
    digits.clear();
    labels.clear();

    String filename = samples::findFile(fn);

    cout << "Loading " << filename << " ..." << endl;

    Mat digits_img = imread(filename, IMREAD_GRAYSCALE);
    split2d(digits_img, Size(SZ, SZ), digits);

    for (int i = 0; i < CLASS_N; i++)
    {
        for (size_t j = 0; j < digits.size() / CLASS_N; j++)
        {
            labels.push_back(i);
        }
    }
}

static void deskew(const Mat& img, Mat& deskewed_img)
{
    Moments m = moments(img);

    if (abs(m.mu02) < 0.01)
    {
        deskewed_img = img.clone();
        return;
    }

    float skew = (float)(m.mu11 / m.mu02);
    float M_vals[2][3] = {{1, skew, -0.5f * SZ * skew}, {0, 1, 0}};
    Mat M(Size(3, 2), CV_32F);

    for (int i = 0; i < M.rows; i++)
    {
        for (int j = 0; j < M.cols; j++)
        {
            M.at<float>(i, j) = M_vals[i][j];
        }
    }

    warpAffine(img, deskewed_img, M, Size(SZ, SZ), WARP_INVERSE_MAP | INTER_LINEAR);
}

static void mosaic(const int width, const vector<Mat>& images, Mat& grid)
{
    int mat_width = SZ * width;
    int mat_height = SZ * (int)ceil((double)images.size() / width);

    if (!images.empty())
    {
        grid = Mat(Size(mat_width, mat_height), images[0].type());

        for (size_t i = 0; i < images.size(); i++)
        {
            Mat location_on_grid = grid(Rect(SZ * ((int)i % width), SZ * ((int)i / width), SZ, SZ));
            images[i].copyTo(location_on_grid);
        }
    }
}

static void evaluate_model(const vector<float>& predictions, const vector<Mat>& digits, const vector<int>& labels, Mat& mos)
{
    double err = 0;

    for (size_t i = 0; i < predictions.size(); i++)
    {
        if ((int)predictions[i] != labels[i])
        {
            err++;
        }
    }

    err /= predictions.size();

    cout << format("error: %.2f %%", err * 100) << endl;

    int confusion[10][10] = {};

    for (size_t i = 0; i < labels.size(); i++)
    {
        confusion[labels[i]][(int)predictions[i]]++;
    }

    cout << "confusion matrix:" << endl;
    for (int i = 0; i < 10; i++)
    {
        for (int j = 0; j < 10; j++)
        {
            cout << format("%2d ", confusion[i][j]);
        }
        cout << endl;
    }

    cout << endl;

    vector<Mat> vis;

    for (size_t i = 0; i < digits.size(); i++)
    {
        Mat img;
        cvtColor(digits[i], img, COLOR_GRAY2BGR);

        if ((int)predictions[i] != labels[i])
        {
            for (int j = 0; j < img.rows; j++)
            {
                for (int k = 0; k < img.cols; k++)
                {
                    img.at<Vec3b>(j, k)[0] = 0;
                    img.at<Vec3b>(j, k)[1] = 0;
                }
            }
        }

        vis.push_back(img);
    }

    mosaic(25, vis, mos);
}

static void bincount(const Mat& x, const Mat& weights, const int min_length, vector<double>& bins)
{
    double max_x_val = 0;
    minMaxLoc(x, NULL, &max_x_val);

    bins = vector<double>(max((int)max_x_val, min_length));

    for (int i = 0; i < x.rows; i++)
    {
        for (int j = 0; j < x.cols; j++)
        {
            bins[x.at<int>(i, j)] += weights.at<float>(i, j);
        }
    }
}

static void preprocess_hog(const vector<Mat>& digits, Mat& hog)
{
    int bin_n = 16;
    int half_cell = SZ / 2;
    double eps = 1e-7;

    hog = Mat(Size(4 * bin_n, (int)digits.size()), CV_32F);

    for (size_t img_index = 0; img_index < digits.size(); img_index++)
    {
        Mat gx;
        Sobel(digits[img_index], gx, CV_32F, 1, 0);

        Mat gy;
        Sobel(digits[img_index], gy, CV_32F, 0, 1);

        Mat mag;
        Mat ang;
        cartToPolar(gx, gy, mag, ang);

        Mat bin(ang.size(), CV_32S);

        for (int i = 0; i < ang.rows; i++)
        {
            for (int j = 0; j < ang.cols; j++)
            {
                bin.at<int>(i, j) = (int)(bin_n * ang.at<float>(i, j) / (2 * CV_PI));
            }
        }

        Mat bin_cells[] = {
            bin(Rect(0, 0, half_cell, half_cell)),
            bin(Rect(half_cell, 0, half_cell, half_cell)),
            bin(Rect(0, half_cell, half_cell, half_cell)),
            bin(Rect(half_cell, half_cell, half_cell, half_cell))
        };
        Mat mag_cells[] = {
            mag(Rect(0, 0, half_cell, half_cell)),
            mag(Rect(half_cell, 0, half_cell, half_cell)),
            mag(Rect(0, half_cell, half_cell, half_cell)),
            mag(Rect(half_cell, half_cell, half_cell, half_cell))
        };

        vector<double> hist;
        hist.reserve(4 * bin_n);

        for (int i = 0; i < 4; i++)
        {
            vector<double> partial_hist;
            bincount(bin_cells[i], mag_cells[i], bin_n, partial_hist);
            hist.insert(hist.end(), partial_hist.begin(), partial_hist.end());
        }

        // transform to Hellinger kernel
        double sum = 0;

        for (size_t i = 0; i < hist.size(); i++)
        {
            sum += hist[i];
        }

        for (size_t i = 0; i < hist.size(); i++)
        {
            hist[i] /= sum + eps;
            hist[i] = sqrt(hist[i]);
        }

        double hist_norm = norm(hist);

        for (size_t i = 0; i < hist.size(); i++)
        {
            hog.at<float>((int)img_index, (int)i) = (float)(hist[i] / (hist_norm + eps));
        }
    }
}

static void shuffle(vector<Mat>& digits, vector<int>& labels)
{
    vector<int> shuffled_indexes(digits.size());

    for (size_t i = 0; i < digits.size(); i++)
    {
        shuffled_indexes[i] = (int)i;
    }

    randShuffle(shuffled_indexes);

    vector<Mat> shuffled_digits(digits.size());
    vector<int> shuffled_labels(labels.size());

    for (size_t i = 0; i < shuffled_indexes.size(); i++)
    {
        shuffled_digits[shuffled_indexes[i]] = digits[i];
        shuffled_labels[shuffled_indexes[i]] = labels[i];
    }

    digits = shuffled_digits;
    labels = shuffled_labels;
}

int main()
{
    help();

    vector<Mat> digits;
    vector<int> labels;

    load_digits(DIGITS_FN, digits, labels);

    cout << "preprocessing..." << endl;

    // shuffle digits
    shuffle(digits, labels);

    vector<Mat> digits2;

    for (size_t i = 0; i < digits.size(); i++)
    {
        Mat deskewed_digit;
        deskew(digits[i], deskewed_digit);
        digits2.push_back(deskewed_digit);
    }

    Mat samples;

    preprocess_hog(digits2, samples);

    int train_n = (int)(0.9 * samples.rows);
    Mat test_set;

    vector<Mat> digits_test(digits2.begin() + train_n, digits2.end());
    mosaic(25, digits_test, test_set);
    imshow("test set", test_set);

    Mat samples_train = samples(Rect(0, 0, samples.cols, train_n));
    Mat samples_test = samples(Rect(0, train_n, samples.cols, samples.rows - train_n));
    vector<int> labels_train(labels.begin(), labels.begin() + train_n);
    vector<int> labels_test(labels.begin() + train_n, labels.end());

    Ptr<ml::KNearest> k_nearest;
    Ptr<ml::SVM> svm;
    vector<float> predictions;
    Mat vis;

    cout << "training KNearest..." << endl;
    k_nearest = ml::KNearest::create();
    k_nearest->train(samples_train, ml::ROW_SAMPLE, labels_train);

    // predict digits with KNearest
    k_nearest->findNearest(samples_test, 4, predictions);
    evaluate_model(predictions, digits_test, labels_test, vis);
    imshow("KNearest test", vis);
    k_nearest.release();

    cout << "training SVM..." << endl;
    svm = ml::SVM::create();
    svm->setGamma(5.383);
    svm->setC(2.67);
    svm->setKernel(ml::SVM::RBF);
    svm->setType(ml::SVM::C_SVC);
    svm->train(samples_train, ml::ROW_SAMPLE, labels_train);

    // predict digits with SVM
    svm->predict(samples_test, predictions);
    evaluate_model(predictions, digits_test, labels_test, vis);
    imshow("SVM test", vis);
    cout << "Saving SVM as \"digits_svm.yml\"..." << endl;
    svm->save("digits_svm.yml");
    svm.release();

    waitKey();

    return 0;
}