// OpenCVTest.cpp : Defines the entry point for the console application.
//

#define _CRT_SECURE_NO_DEPRECATE
#define _USE_MATH_DEFINES


#include "stdafx.h"
#include <iostream>
#include <opencv2/core.hpp>
#include <stdio.h>
//#include "cxcore.hpp"
#include <opencv\cxcore.h>
#include <opencv2/imgcodecs.hpp>
#include <conio.h>
//#include "imgproc.hpp"
#include <opencv2/imgproc.hpp>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui.hpp>
//#include "highgui.hpp"
#include <stdio.h>

#include <fstream>
#include <math.h>



/// Shortcuts for std and opencv
using namespace std;
using namespace cv;

/// Global variables
#define SegmentSizeX 58//Size of segments in x
#define SegmentSizeY 47//Size of segments in y

float OfficialPercentile; //Percentile size
int i;
int j;

///* Histogram creation *///

void myhist(Mat &med, Mat &decoy) {
    decoy.setTo(0);
    float maxi = 0;

    float *y = decoy.ptr<float>(0);
    for (int i = 0; i < med.rows; i++)
    {
        uchar *x = med.ptr<uchar>(i);
        for (int j = 0; j < med.cols; j++)
        {
            for (int k = 0; k < 256; k++)
            {
                if (k == (int)x[j])
                {
                    y[k] += 1;
                    if (maxi < y[k])
                    {
                        maxi = y[k];
                    }
                    break;
                }
            }
        }
    }
    y[256] = maxi;
}


///* Function that draws the histogram *///
void drawHist(Mat &decoy, char *winName)
{
    Mat histFig(640, 512, CV_8U);
    histFig.setTo(256);
    namedWindow(winName);

    float *y = decoy.ptr<float>(0);
    float max = y[256];
    float factor = ((float)640 / max);

    for (int i = 0; i < 256; i++)
    {
        y[i] *= factor;
        Point pt1, pt2;
        pt1 = Point(2 * i, 640);
        pt2 = Point(2 * i, 640 - (int)y[i]);
        const Scalar col = Scalar(0, 0, 0);
        line(histFig, pt1, pt2, col, 2, CV_AA);

    }
    imshow(winName, histFig);
}


///* Function that writes matrices to a file *///
void writeMatToFile(Mat &m, const char* filename)
{
    ofstream fout(filename);

    if (!fout)
    {
        cout << "File Not Opened" << endl;
        return;
    }

    for (int i = 0; i<m.rows; i++)
    {
        for (int j = 0; j<m.cols; j++)
        {
            fout << m.at<float>(i, j) << "\t";
        }
        fout << endl;
    }

    fout.close();
}

void segmentDivisions(Mat &decoyM, Mat *SegmentArray) {
    int i;
    int j;
    int k = 0;

    //printf("cols: %d \n",decoyM.cols);
    //printf("rows: %d \n", decoyM.rows);
    for (i = 0; i < decoyM.rows; i += SegmentSizeX) {
        for (j = 0; j < decoyM.cols; j += SegmentSizeY) {
            //printf("Test inside void i:%d j:%d k:%d ", i, j, k);
            Rect roiArr(j, i, SegmentSizeY, SegmentSizeX);
            Mat holdingRect = decoyM(roiArr);
            k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
            SegmentArray[k] = holdingRect.clone();

        }
    }

}

void StarMatrix(Mat &decoyM, Mat *SegmentStar) {
    int i;
    int j;
    int k = 0;

    //printf("cols: %d \n",decoyM.cols);
    //printf("rows: %d \n", decoyM.rows);
    for (i = 0; i < decoyM.rows; i += SegmentSizeX) {
        for (j = 0; j < decoyM.cols; j += SegmentSizeY) {
            //printf("Test inside void i:%d j:%d k:%d ", i, j, k);
            Rect roiArr(j, i, SegmentSizeY, SegmentSizeX);
            Mat holdingRect = decoyM(roiArr);
            k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
            SegmentStar[k] = holdingRect.clone();
        }
    }
}


void PercentileCalc(Mat &decoy) {
    int total = 0;
    int total2 = 0;
    int total3 = 0;
    float percentile;
    int number;

    /// Acumulates the intensity (values) of the picture
    for (i = 0; i < 256; i++)
    {
        total += decoy.at<float>(i);
    }

    /// Defining the 40% percentile
    float percenti = total * 0.4; // 0.4 = 40%
                                  //printf("The 40% of the accumuated value is \n %f", percenti);

                                  /// Need this function as we using i-1, to avoid negative values
    for (i = 0; i < 1; i++)
    {
        total2 += decoy.at<float>(i);
    }

    for (i = 1; i < 256; i++)
    {
        /// Checks for what itensity and which color-value is bounded by the accumulated values
        /// of the percentile
        if (percenti > total2)
        {
            /// Acumulates the value as long as its not more than the percenti
            /// This value will eventually become larger than the percenti (upperbounded)
            total2 += decoy.at<float>(i);

            /// Accumulates the value, one position less than the previous one
            /// This function will not surpass the percenti value (lowerbounded)
            total3 += decoy.at<float>(i - 1);

            /// This value will be the color-value for which is correspond to the lower-bounded value
            number = i;
        }
    }

    /// Function to calulate the value between the lower and upper-bounded values
    percentile = (percenti - total3) / (total2 - total3);
    //printf("\n Percenti = %f \n Total 2 = %d \n Total 3 = %d \n", percenti, total2, total3);

    /// This function sums the value found to the lower color-number, thus this will yield the
    /// exact color value for which the percentile value exist
    OfficialPercentile = number + percentile;
    printf("The percentile is found for colorvalue %f", number + percentile);
}

void RemoveBackground(Mat &SegM, float pres) {
    for (int i = 0; i < SegmentSizeX; i++) {
        for (int j = 0; j < SegmentSizeY; j++) {
            if (SegM.at<float>(i, j) <= pres) {
                //SegM->at<float>(i, j) = 0;
            }
        }
    }
    //printf("Value of point is %d \n", data[i]);
    //printf("Value of point is %d \n", data[j]);

    /*
    for (i = 0; i < width_final; i++) {
    for (j = 0; j < height_final; j++) {

    if (data[i] < 15 ) {
    data[i] = 0;
    }
    if (data[j] < 15) {
    data[j] = 0;
    }
    //printf("Value of point is %d \n", data[i]);
    //printf("Value of point is %d \n", data[j]);
    }
    }
    */
}

int main(int argc, char* argv[])
{
    //Mat(MAX_EGDE, MAX_EGDE, CV_64F);

    /// Assigning size of the Matrix-Arrays
    Mat SegmentArray[(752 / SegmentSizeY)*(580 / SegmentSizeX)];
    Mat StarMatrixArray[144];



    /// Uploading binary file and making a matrix out of it
    FILE* f = fopen("C:\\Users\\Runar\\Desktop\\DTU Work\\Image_Ana\\Prosject\\Images\\ICON C1 RS2\\A150722_12132639.unc", "rb");
    unsigned char info[34];
    fread(info, sizeof(unsigned char), 34, f);

    /// Extract image height and width from header
    int width = *(int*)&info[28];
    int height = *(int*)&info[30];

    /// Defining the length of the picture
    unsigned int width_h = info[31];
    unsigned int width_l = info[30];
    unsigned int width_final = (width_h << 8) + width_l;
    //printf("width %d \n", width_final);

    unsigned int height_h = info[29];
    unsigned int height_t = info[28];
    unsigned int height_final = 2 * ((height_h << 8) + height_t); //Multiplied with 2 (see Image format description)
                                                                  //printf("height %d \n", height_final);

                                                                  //To see values in the header file, uncomment this section
                                                                  /*
                                                                  for (i = 0; i < 34; i++) {
                                                                  printf("%d \n", info[i]);
                                                                  }
                                                                  printf("2 \n");
                                                                  */

                                                                  /// Writing of the binary data to a matrix continues
    int size = 1 * width_final * height_final;
    unsigned char* data = new unsigned char[size]; // allocate 3 bytes per pixel
    fread(data, sizeof(unsigned char), size, f); // read the rest of the data at once
    fclose(f);

    Mat image = Mat(width_final*height_final, 1, CV_8U, (unsigned*)data);
    Mat reshaped = image.reshape(0, height_final);
    cvNamedWindow("image", WINDOW_AUTOSIZE);
    imshow("image", reshaped);

    /// Snappshot of a certain area of the picture
    Mat outImg;
    Rect roi(730, 295, 8, 15);
    Mat Rectangle_outImg = reshaped(roi);

    resize(Rectangle_outImg, outImg, cv::Size(reshaped.cols*1.5, reshaped.rows*1.5), 0, 0, CV_INTER_LINEAR);
    cvNamedWindow("resized window", CV_WINDOW_AUTOSIZE);
    imshow("resized window", outImg);

    /// To check the amount of pixels in the new image, a picture is saved in the workspace
    /// and then by checking properties of the picture in the folder, pixel number can be validated

    //imwrite("C:\\Users\\Runar\\Desktop\\DTU Work\\Image_Ana\\Prosject\\Images\\ICON C1 RS2\\EdgdeOfStar",outImg);




    /* Functions should start here */

    /// Calculates the segmentsize, saving it as int (whole number)
    int SegmentSize = (width_final / SegmentSizeY)*(height_final / SegmentSizeX);

    ///* Step 2: Calculations of Segments *///

    /// MERGINGSEGEMENTS:
    // There is several way to do this, but the best way I found is to use cv::hconcat(mat1, mat2, dst)
    // for horizontal merge orcv::vconcat(mat1, mat2, dst) for vertical.

    //printf("SegmentArray: %d", (750 / SegmentSize)*(580 / SegmentSize));

    segmentDivisions(reshaped, SegmentArray);
    //Reading out to a file
    //Mat TestMatr = Mat(100, 100, CV_32F);
    cout << SegmentArray[115] << endl;


    ///* Step 2 Calculating the brigthest star within each segments*///
    /*
    ///Function that defines the max, min and its position
    double Min, Max, MinSeg, MaxSeg;
    Point p_min, p_max;

    //Save min and max as an array?

    Mat StarOutMat;
    Mat mergedMatrix;

    int threshold = 35;

    for (i = 8; i < 9; i++) {

        ///Finding the max, min value and the position
        //cout << SegmentArray[i] << endl;
        minMaxLoc(SegmentArray[i], &Min, &Max, &p_min, &p_max);
        printf("For segment %d : ", i);
        cout << "min: " << Min << " at " << p_min << "max: " << Max << " at " << p_max << endl;

        int MaxValue = Max;

        if (MaxValue > threshold) {

            ///Setting region of interest
            int xUpperBorder = p_max.x + 6;
            int xLowerBorder = p_max.x - 6;
            int yUpperBorder = p_max.y + 6;
            int yLowerBorder = p_max.y - 6;

            /// Have to divide per 16 (y) segment 9 (x)

            /// If the star is at an end, we have to merge two matrices


            if ((xUpperBorder > SegmentSizeY && yUpperBorder > SegmentSizeX) || (xLowerBorder < 0 && yLowerBorder < 0)) {

                /// If a star is in a corner, 4 matrices (or new matrix) has to be made
                cout << "This star is out of reach" << endl;
            }

            if (xUpperBorder > SegmentSizeY) {

                if (i % 47 == 0) {
                    cout << "The region of interest is outside of the picture" << endl;
                }

                else {

                    hconcat(SegmentArray[i], SegmentArray[i + 1], mergedMatrix);
                    cout << "Fails for x1 coordinate" << endl;

                    /// Resizing and interpolating the segment and defining start windowd
                    resize(mergedMatrix, StarOutMat, cv::Size(), 2, 2, INTER_LINEAR);

                    /// Calculates the same procedure for new matrix
                    minMaxLoc(StarOutMat, &Min, &Max, &p_min, &p_max);
                    cout << "min: " << Min << " at " << p_min << "max: " << Max << " at " << p_max << endl;

                    int New_yCoord = p_max.y - 6;
                    int New_xCoord = p_max.x - 6;
                    printf(" \n new ycoord %d and new xcoord %d \n", New_yCoord, New_xCoord);

                    //cout << StarOutMat << endl;

                    /// Making a matrix around the star
                    Rect roiStar(New_xCoord, New_yCoord, 12, 12);
                    Mat HoldingMatrix = StarOutMat(roiStar);
                    Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
                                                            //k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
                                                            //SegmentStar[k] = holdingRect.clone();
                    cout << StarMatrix << endl;
                    StarMatrixArray[i] = StarMatrix;

                }
            }

            if (xLowerBorder < 0) {

                if (i % 46 == 0) {
                    cout << "The region of interest is outside of the picture" << endl;
                }

                else {

                    hconcat(SegmentArray[i - 1], SegmentArray[i], mergedMatrix);
                    cout << "Fails for x_2 coordinate" << endl;

                    /// Resizing and interpolating the segment and defining start windowd
                    resize(mergedMatrix, StarOutMat, cv::Size(), 2, 2, INTER_LINEAR);

                    /// Calculates the same procedure for new matrix
                    minMaxLoc(StarOutMat, &Min, &Max, &p_min, &p_max);
                    cout << "min: " << Min << " at " << p_min << "max: " << Max << " at " << p_max << endl;

                    int New_yCoord = p_max.y - 6;
                    int New_xCoord = p_max.x - 6;
                    printf(" \n new ycoord %d and new xcoord %d \n", New_yCoord, New_xCoord);

                    //cout << StarOutMat << endl;

                    /// Making a matrix around the star
                    Rect roiStar(New_xCoord, New_yCoord, 12, 12);
                    Mat HoldingMatrix = StarOutMat(roiStar);
                    Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
                                                            //k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
                                                            //SegmentStar[k] = holdingRect.clone();
                    cout << StarMatrix << endl;
                    StarMatrixArray[i] = StarMatrix;

                    printf("xLowerBorder is now %d", xLowerBorder);
                }


            }

            if (yUpperBorder > SegmentSizeX) {

                if (i < 16) {
                    cout << "The region of interest is outside of the picture" << endl;
                }

                else {

                    vconcat(SegmentArray[i], SegmentArray[i-16], mergedMatrix);
                    cout << "Fails for yUpperBorder coordinate" << endl;

                    /// Resizing and interpolating the segment and defining start windowd
                    resize(mergedMatrix, StarOutMat, cv::Size(), 2, 2, INTER_LINEAR);

                    /// Calculates the same procedure for new matrix
                    minMaxLoc(StarOutMat, &Min, &Max, &p_min, &p_max);
                    cout << "min: " << Min << " at " << p_min << "max: " << Max << " at " << p_max << endl;

                    int New_yCoord = p_max.y - 6;
                    int New_xCoord = p_max.x - 6;
                    printf(" \n new ycoord %d and new xcoord %d \n", New_yCoord, New_xCoord);

                    //cout << StarOutMat << endl;

                    /// Making a matrix around the star
                    Rect roiStar(New_xCoord, New_yCoord, 12, 12);
                    Mat HoldingMatrix = StarOutMat(roiStar);
                    Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
                                                            //k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
                                                            //SegmentStar[k] = holdingRect.clone();
                    cout << StarMatrix << endl;
                    StarMatrixArray[i] = StarMatrix;

                    printf("Upperborder is now %d", xLowerBorder);

                }
            }

            if (yLowerBorder < 0) {

                if (i > 144 ){
                    cout << "The region of interest is outside of the picture" << endl;
                }

                else {

                    vconcat(SegmentArray[i], SegmentArray[i + 16], mergedMatrix);
                    cout << "Fails for yLowerBorder coordinate" << endl;

                    /// Resizing and interpolating the segment and defining start windowd
                    resize(mergedMatrix, StarOutMat, cv::Size(), 2, 2, INTER_LINEAR);

                    /// Calculates the same procedure for new matrix
                    minMaxLoc(StarOutMat, &Min, &Max, &p_min, &p_max);
                    cout << "min: " << Min << " at " << p_min << "max: " << Max << " at " << p_max << endl;

                    int New_yCoord = p_max.y - 6;
                    int New_xCoord = p_max.x - 6;
                    printf(" \n new ycoord %d and new xcoord %d \n", New_yCoord, New_xCoord);

                    //cout << StarOutMat << endl;

                    /// Making a matrix around the star
                    Rect roiStar(New_xCoord, New_yCoord, 12, 12);
                    Mat HoldingMatrix = StarOutMat(roiStar);
                    Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
                                                            //k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
                                                            //SegmentStar[k] = holdingRect.clone();
                    cout << StarMatrix << endl;
                    StarMatrixArray[i] = StarMatrix;

                    printf("yUpperBorder is now %d", xLowerBorder);

                }
            }

            else {
                cout << "TEST TEST TEST \n\n" << endl;

                /// Resizing and interpolating the segment and defining start windowd
                resize(SegmentArray[i], StarOutMat, cv::Size(), 2, 2, INTER_LINEAR);
                //cout << StarOutMat << endl;

                int yCoord = p_max.y - 6;
                int xCoord = p_max.x - 6;

                /// Making a matrix around the star
                Rect roiStar(xCoord, yCoord, 12, 12);
                Mat HoldingMatrix = StarOutMat(roiStar);
                Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
                //k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
                //SegmentStar[k] = holdingRect.clone();
                StarMatrixArray[i] = StarMatrix;

                //printf("This segment works %d", i);
                //cout << StarMatrixArray[i] << endl;
            }
        }
    }
    */




    /*
    resize(SegmentArray[41], StarOutMat, cv::Size(), 1, 1, INTER_CUBIC);
    minMaxLoc(StarOutMat, &MinSeg, &MaxSeg, &p_min, &p_max);
    int yCoord = p_max.y - 6;
    int xCoord = p_max.x - 6;

    printf("yCoord is %d and xCoord is %d ", yCoord, xCoord);
    printf("\n\n");

    Rect roiStar(xCoord, yCoord, 12, 12);
    Mat HoldingMatrix = StarOutMat(roiStar);
    Mat StarMatrix = HoldingMatrix.clone(); //Clones it in a matrix so original variables dont change
    cout << StarMatrix << endl;
    */






    //myhist(StarMatrix, decoy);
    //drawHist(decoy, "Histogram without interpolation");
    //const char* filename = "StarOut.txt";
    //writeMatToFile(StarMatrix, filename);
    //Rect roiStar2(xCoord, yCoord, 12, 12);
    //Mat StarMatrix2 = StarOutMat(roiStar);

    // Get the values
    //printf("\ntestPoint is : %d\n", testPoint);
    //printf("MaxSeg is this number: %f", p_max); // 7 - 47


    int k = 0;



    //const char* filename = "Output.txt";
    //writeMatToFile(SegmentArray[4], filename);





    ///* Step 3: Calculation of Histogram *///

    /*if (argc < 2)
    {
    //std::cout << "less arguments" << std::end1;
    printf("Fails here");
    return -1;
    }*/

    //Reading out to a file
    //const char* filename = "output3.txt";
    //writeMatToFile(SegmentArray[], filename);

    /// Define the limit of the histogram, calculates and draws it
    Mat decoy(1, 257, CV_32F);
    //for (i = 0; i < 10; i++) {

    printf("\n");
    myhist(SegmentArray[10], decoy);
    drawHist(decoy, "Histogram of section");




    //imwrite("C:\\Users\\Runar\\Desktop\\DTU Work\\Image_Ana\\Prosject\\Images\\WorkInPorgress\\Histogram.png", decoy);
    //}

    ///* Step 4: Calculation of Percentile *///
    PercentileCalc(decoy);
    printf("Percentile %f", OfficialPercentile);
    //count <

    /// Step 5: Calculating the removal of background *///

    Mat Test101 = SegmentArray[22];

    RemoveBackground(Test101, OfficialPercentile); // SegM = Seg(k) and pres = pre
    cout << Test101 << endl;

    //const char* filename = "Output.txt";
    //writeMatToFile(Test101, filename);






    /*Mat holdingRect = decoyM(roiArr);
    k = floor((decoyM.cols / SegmentSizeY)*(i / SegmentSizeX) + j / SegmentSizeY);
    SegmentArray[k] = holdingRect.clone();
    */
    ///* Calculation of the cart-2-polar *///

    /// Average of the first picture
    /// These variables should be changed to the average / or updated per picture
    int averageX = 291; //midle area for j variable
    int averageY = 367; //midle area for i variable

                        /// Variables used to for cart-2-polar transformation
    double magnitude;
    double angle;
    double distanceX;
    double distanceY;

    /// Size of arrays, spesified for area of interest
    float Avalue[15][8];
    float Mvalue[15][8];

    /// A area is manually found, and arrays of their magnitude and angle is created
    for (int i = 295; i < 310; i++)
    {
        /// Is uesed if we want to utilize the intensity (value) of each byte
        unsigned char* pByte = reshaped.ptr<unsigned char>(i);

        for (int j = 730; j < 738; j++) {

            /// We need to check if the area is in the RHP, LFH and what complex plane its in
            if (j > averageX) //Then we need abs sign
            {
                distanceX = abs(averageX - j);
            }
            if (j < averageX)
            {
                distanceX = (averageX - j);
            }
            if (i > averageY) //Then we need abs sign
            {
                distanceY = abs(averageY - i);
            }
            if (i < averageY)
            {
                distanceY = (averageY - i);
            }

            /// Printing out the angle and distance
            angle = atan2(distanceY, distanceX);
            magnitude = sqrt(abs(distanceX*distanceX) + abs(distanceY*distanceY));

            /// Need to substract from the area we are cuurrently looking at
            Avalue[i - 295][j - 730] = angle;
            Mvalue[i - 295][j - 730] = magnitude;


            /// Function if converting to degrees is desired
            /*
            double rad2angle;
            rad2angle = angle*(180 / (M_PI));
            */

            /// If the desired values is liked to be printed out in the console application
            /*
            printf("Mag and angle: %f & %f with value: ", magnitude, angle);
            printf("value %d ", pByte[j]);
            printf("\n");
            */
        }
    }


    /// The folling function stores the value in seperate matrices
    /*
    Mat AngelMatrix = Mat(15,8, CV_32F, Avalue);
    Mat MagnitudeMatrix = Mat(15, 8, CV_32F, Mvalue);
    printf("\n");
    //To print Mat values:
    //cout << MagnitudeMatrix.at<float>(0, 2) << endl;
    cout << MagnitudeMatrix << endl;
    printf("\n");
    //cout << AngelMatrix.at<float>(0, 2) << endl;
    cout << AngelMatrix << endl;
    */







    //imshow("Removed background picture", threshold_Image);
    //cout << threshold_Image << endl;



    waitKey(0);
    return(0);





}