math_mod_lab2

// TEST_OPENCV.cpp: определяет точку входа для консольного приложения.
//

#include "stdafx.h"
#include <bits\stdc++.h>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <conio.h>
#include <windows.h>

using namespace cv;
using namespace std;

vector<vector<double>> Intensity(Mat & img){
    vector<vector<double>> a(img.rows , vector<double>(img.cols));
    for(int i = 0 ; i < img.rows ; ++i){
        for(int j = 0 ; j < img.cols ; ++j){
            a[i][j] = (int)img.at<uchar>(i,j);
        }
    }
    return a;
}

vector<vector<double>> Transpose(vector<vector<double>> & m){
    assert(m.size() > 0);
    vector<vector<double>> a(m[0].size() , vector<double>(m.size()));
    for(int i = 0 ; i < m.size() ; ++i){
        for(int j = 0 ; j < m[0].size() ; ++j){
            a[j][i] = m[i][j];
        }
    }
    return a;
}

vector<vector<double>> Multiply(vector<vector<double>> & a , vector<vector<double>> & b){
    assert(a.size() && b.size());
    assert(a.size() == b[0].size());
    vector<vector<double>> m(a.size() , vector<double>(b[0].size()));
    for(int i = 0 ; i < a.size() ; ++i){
        for(int j = 0 ; j < b[0].size() ; ++j){
            for(int k = 0 ; k < b.size() ; ++k){
                m[i][j] += a[i][k] * b[k][j];
            }
        }
    }
    return m;
}

vector<vector<double>> Add(vector<vector<double>> & a , vector<vector<double>> & b , double k = 1){// A + kB
    assert(a.size() && a.size() == b.size());
    assert(a[0].size() == b[0].size());
    vector<vector<double>> m(a.size() , vector<double>(a[0].size()));
    for(int i = 0 ; i < a.size() ; ++i){
        for(int j = 0 ; j < a[i].size() ; ++j){
            m[i][j] += a[i][j] + k * b[i][j];
        }
    }
    return m;
}

vector<vector<double>> MultK(vector<vector<double>> & a , double k){
    assert(a.size());
    assert(a[0].size());
    vector<vector<double>> m(a.size() , vector<double>(a[0].size()));
    for(int i = 0 ; i < a.size() ; ++i){
        for(int j = 0 ; j < a[i].size() ; ++j){
            m[i][j] = a[i][j] / k;
        }
    }
    return m;
}

vector<vector<double>> Inverse(vector<vector<double>> & a){
    assert(a.size());
    assert(a.size() == a[0].size());
    int n = a.size();
    vector<vector<double>> m(n , vector<double>(2 * n , 0.0));
    for(int i = 0 ; i < n ; ++i){
        for(int j = 0 ; j < n ; ++j){
            m[i][j] = a[i][j];
        }
        m[i][i + n] = 1.0;
    }
    for(int i = 0 ; i < n ; ++i){
        double t = m[i][i];
        if(fabs(t) < 1e-7)return vector<vector<double>>(0);
        for(int j = i ; j < 2 * n ; ++j){
            m[i][j] = m[i][j] / t;
        }
        for(int j = 0 ; j < n ; ++j){
            if(i != j){
                t = m[j][i];
                for(int k = 0 ; k < 2 * n ; ++k){
                    m[j][k] = m[j][k] - t * m[i][k];
                }
            }
        }
    }
    vector<vector<double>> res(n , vector<double>(n , 0.0));
    for(int i = 0 ; i < n ; ++i){
        for(int j = 0 ; j < n ; ++j){
            res[i][j] = m[i][j + n];
        }
    }
    return res;
}

vector<vector<double>> Moore_Penrose(vector<vector<double>> & a , double eps){
    vector<vector<double>> xt = Transpose(a);
    vector<vector<double>> xxt = Multiply(xt , a);
    vector<vector<double>> e(xxt.size() , vector<double>(xxt.size() , 0.0));
    for(int i = 0 ; i < xxt.size() ; ++i){
        e[i][i] = 1.0;
    }
    vector<vector<double>> ad = Add(xxt , e , eps);
    vector<vector<double>> iad = Inverse(ad);
    vector<vector<double>> res = Multiply(iad , xt);
    return res;
}

double norm(vector<vector<double>> & a){
    double sum = 0.0;
    for(int i = 0 ; i < a.size(); ++i){
        for(int j = 0 ; j < a[i].size() ; ++j){
            sum += pow(a[i][j] , 2);
        }
    }
    return sum;
}

vector<vector<double>> Solve_Moore_Penrose(vector<vector<double>> & a , vector<vector<double>> & b){//vector<vector<double>>
    double best_norm = 1e9 , best_eps = 1e9;
    vector<vector<double>> px;
    for(int i = 10 ; i <= 100 ; i += 10){
        double eps = (double)i / 100.0;
        vector<vector<double>> aa = Moore_Penrose(a , eps);
        vector<vector<double>> t = Multiply(b , aa);
        vector<vector<double>> ta = Multiply(t , a);
        vector<vector<double>> ck = Add(ta , b , -1.0);
        double cur = norm(ck);
        if(cur < best_norm){
            best_norm = cur;
            best_eps = eps;
            px = aa;
        }
    }
    vector<vector<double>> A = Multiply(b , px);
    return Multiply(A , a);
}


vector<vector<double>> Greville(vector<vector<double>> & a){
    vector<vector<double>> x(1 , vector<double>(a[0].size()));
    x[0] = a[0];

    vector<vector<double>> tr = Transpose(x);
    vector<vector<double>> div = Multiply(x , tr);
    vector<vector<double>> pin = MultK(tr , div[0][0]);

    for(int i = 1 ; i < a.size() ; ++i){
        vector<vector<double>> cur(1 , vector<double>(a[i].size()));
        cur[0] = a[i];
        vector<vector<double>> tcur = Transpose(cur);
        vector<vector<double>> e(a[0].size() , vector<double>(a[0].size() , 0.0));
        for(int j = 0 ; j < a[0].size() ; ++j){
            e[j][j] = 1.0;
        }
        vector<vector<double>> su = Multiply(pin , x);
        vector<vector<double>> z = Add(e , su , -1);
        vector<vector<double>> ck = Multiply(cur , z);
        ck = Multiply(ck , tcur);

        if(ck[0][0] > 1e-6){
            double d = ck[0][0];
            vector<vector<double>> fi = Multiply(z , tcur);
            fi = Multiply(fi , cur);
            fi = Multiply(fi , pin);
            fi = MultK(fi , d);
            fi = Add(pin , fi , -1);

            vector<vector<double>> se = Multiply(z , tcur);
            se = MultK(se , d);

            pin = fi;
            for(int j = 0 ; j < pin.size() ; ++j){
                pin[j].push_back(se[j][0]);
            }
        }else{
            vector<vector<double>> tpin = Transpose(pin);
            vector<vector<double>> r = Multiply(pin , tpin);
            vector<vector<double>> ck2 = Multiply(cur , r);
            ck2 = Multiply(ck , tcur);
            double d = 1 + ck2[0][0];
            vector<vector<double>> fi = Multiply(r , tcur);
            fi = Multiply(fi , cur);
            fi = Multiply(fi , pin);
            fi = MultK(fi , d);
            fi = Add(pin , fi , -1);

            vector<vector<double>> se = Multiply(r , tcur);
            se = MultK(se , d);

            pin = fi;
            for(int j = 0 ; j < pin.size() ; ++j){
                pin[j].push_back(se[j][0]);
            }
        }
        x.push_back(a[i]);
    }
    return pin;
}


int round(double d){
    double d1 = floor(d);
    if(fabs(fabs(d) - fabs(d1)) - 0.5 >= 1e-7){
        return max(0 , min((int)d + 1 , 255));
    }
    return max(0 , min((int)d , 255));
}

int main( int argc, char** argv )
{
//  freopen("output.txt" , "w" , stdout);
    cout.setf(ios::fixed);
    cout.precision(14);
    Mat X , Y;
    X = imread("C:\\Users\\Max\\Downloads\\Lab_mat_model_student\\Lab2\\x1.bmp" , CV_LOAD_IMAGE_ANYDEPTH);
    Y = imread("C:\\Users\\Max\\Downloads\\Lab_mat_model_student\\Lab2\\y8.bmp" , CV_LOAD_IMAGE_ANYDEPTH);
    if(!X.data || !Y.data){
        cout <<  "Could not open or find images" << endl;
        system("pause");
        return -1;
    }

    vector<vector<double>> intensity_x = Intensity(X);
    vector<vector<double>> intensity_y = Intensity(Y);

    intensity_x.push_back(vector<double>(intensity_x[0].size() , 1));

    vector<vector<double>> A = Solve_Moore_Penrose(intensity_x , intensity_y);

    Mat res_Moore_Penrose(A.size() , A[0].size() , CV_8UC1);

    for(int i = 0 ; i < A.size() ; ++i){
        for(int j = 0 ; j < A[i].size() ; ++j){
            res_Moore_Penrose.at<uchar>(i,j) = round(A[i][j]);
        }
    }


    vector<vector<double>> Agr = Greville(intensity_x);
    Agr = Multiply(intensity_y , Agr);
    Agr = Multiply(Agr , intensity_x);

    Mat res_Greville(Agr.size() , Agr[0].size() , CV_8UC1);

    for(int i = 0 ; i < Agr.size() ; ++i){
        for(int j = 0 ; j < Agr[i].size() ; ++j){
            res_Greville.at<uchar>(i,j) = round(Agr[i][j]);
        }
    }
    namedWindow("Moore_Penrose", WINDOW_AUTOSIZE);
    imshow("Moore_Penrose", res_Moore_Penrose);
    waitKey(0);
    namedWindow("Greville", WINDOW_AUTOSIZE);
    imshow("Greville", res_Greville);
    waitKey(0);
    system("pause");
    return 0;
}