#include <bits/stdc++.h>
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/video.hpp"
#include "opencv2/objdetect.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/ml.hpp"
#define inf 2333333333333333
#define N 100010
#define p(a) putchar(a)
#define For(i,a,b) for(int i=a;i<=b;++i)
//by war
//2020.11.19
using namespace std;
using namespace cv;
Mat image0,image1,Ix,Iy,Ix_2,Iy_2,IxIy,out,dst_norm, dst_norm_scaled;
double alpha = 0.05;
void in(int &x){
    int y=1;char c=getchar();x=0;
    while(c<'0'||c>'9'){if(c=='-')y=-1;c=getchar();}
    while(c<='9'&&c>='0'){ x=(x<<1)+(x<<3)+c-'0';c=getchar();}
    x*=y;
}
void o(int x){
    if(x<0){p('-');x=-x;}
    if(x>9)o(x/10);
    p(x%10+'0');
}

signed main(){
    int ksize = 5;
    int scale = 1;
    int delta = 0;
    int ddepth = CV_16S;
    image0 = imread("/Users/war/Downloads/nk.jpg",0);
    image1 = imread("/Users/war/Downloads/nk.jpg");
    Sobel(image0, Ix, ddepth, 1, 0, ksize, scale, delta, BORDER_DEFAULT);
    Sobel(image0, Iy, ddepth, 1, 0, ksize, scale, delta, BORDER_DEFAULT);
    IxIy = Mat::zeros(Ix.rows, Ix.cols, Ix.type());
    For(i,0,Ix.rows-1)
        For(j,0,Iy.cols-1)
            For(k,0,2){
                IxIy.at<Vec3b>(i,j)[k] = Ix.at<Vec3b>(i,j)[k]*Iy.at<Vec3b>(i,j)[k];
            }
    
    Ix_2 = Mat::zeros(Ix.rows, Ix.cols, Ix.type());
    For(i,0,Ix.rows-1)
        For(j,0,Iy.cols-1)
            For(k,0,2){
                Ix_2.at<Vec3b>(i,j)[k] = Ix.at<Vec3b>(i,j)[k]*Ix.at<Vec3b>(i,j)[k];
            }

    Iy_2 = Mat::zeros(Ix.rows, Ix.cols, Ix.type());
    For(i,0,Ix.rows-1)
        For(j,0,Iy.cols-1)
            For(k,0,2){
                Iy_2.at<Vec3b>(i,j)[k] = Iy.at<Vec3b>(i,j)[k]*Iy.at<Vec3b>(i,j)[k];
            }
    GaussianBlur(Ix_2,Ix_2,Size(5,5),1);
    GaussianBlur(Iy_2,Iy_2,Size(5,5),1);
    GaussianBlur(IxIy,IxIy,Size(5,5),1);
    out = Mat::zeros(Ix.rows, Ix.cols, Ix.type());
    For(i,0,Ix.rows-1)
        For(j,0,Iy.cols-1)
            For(k,0,2){
                out.at<Vec3b>(i,j)[k] = Ix_2.at<Vec3b>(i,j)[k]*Iy_2.at<Vec3b>(i,j)[k]-IxIy.at<Vec3b>(i,j)[k]*IxIy.at<Vec3b>(i,j)[k]-
                alpha*((Ix_2.at<Vec3b>(i,j)[k]+Iy_2.at<Vec3b>(i,j)[k])*(Ix_2.at<Vec3b>(i,j)[k]+Iy_2.at<Vec3b>(i,j)[k]));
            }
    threshold(out, out, 200, 255, THRESH_BINARY);
    imshow("out1",out);
//    For(i,0,Ix.rows-1)
//        For(j,0,Iy.cols-1)
//            if((int) out.at<Vec3b>(i,j)[0] > 240 && (int) out.at<Vec3b>(i,j)[1] > 240 && (int) out.at<Vec3b>(i,j)[2] > 240){
//                circle( image1, Point(j,i),0.0001, Scalar(0,0,255), 3, LINE_AA );
//            }
//    imshow("out",image1);
    waitKey();
    return 0;
}