//**************************************************************************
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include "matrix.h"
#include "edgedetect.h"
//**************************************************************************
#define MAX_ITER 0
//**************************************************************************
void DetectEdge(Matrix *Edge, Matrix &Fr, Matrix &Fi)
{
cout << "Edge Detection begins...." << endl;
cout.flush();

int height = Fr.nx;
int width = Fr.ny;

Matrix Fn_r(height, width), Fn_i(height, width);
Matrix Fnn_r(height, width), Fnn_i(height, width);

Fn_r = Fr;
Fn_i = Fi;

int n = 0;
while (n < MAX_ITER)
        {
         ++n;
         Fnn_r.Zero();
         Fnn_i.Zero();
         int count = 0;
         for (int i = 1; i < height-1; i++)
                for (int j = 1; j < width-1; j++)
                        {
                         double angle = atan2(Fn_i.p[i][j], Fn_r.p[i][j]);
                         if (angle < 0.0) angle += 2.0*Pi;
                         int quadrant = ((int) rint(angle/(Pi/4.0)))%8;

                         int i_incr = quadrant%4;
                         if (i_incr && quadrant < 4)
                                i_incr = +1;
                         else if (i_incr && quadrant > 4)
                                i_incr = -1;

                         quadrant = (quadrant+2)%8;
                         int j_incr = quadrant%4;
                         if (j_incr && quadrant < 4)
                                j_incr = +1;
                         else if (j_incr && quadrant > 4)
                                j_incr = -1;
                         
                         int in = i + i_incr;
                         int jn = j + j_incr;
                         double dotp = Fn_r.p[i][j]*Fn_r.p[in][jn]+Fn_i.p[i][j]*Fn_i.p[in][jn];
                         double f1 = sqrt(Fn_r.p[i][j]*Fn_r.p[i][j]+Fn_i.p[i][j]*Fn_i.p[i][j]);
                         double f2 = sqrt(Fn_r.p[in][jn]*Fn_r.p[in][jn]+Fn_i.p[in][jn]*Fn_i.p[in][jn]);
                         dotp /= f1*f2;
                         double angle2 = atan2(Fn_i.p[i][j]-Fn_i.p[in][jn],Fn_r.p[i][j]-Fn_r.p[in][jn]);
                         if (dotp > 0.0)
                                {
                                 ++count;
                                 Fnn_r.p[in][jn] += Fn_r.p[i][j];
                                 Fnn_i.p[in][jn] += Fn_i.p[i][j];
                                }
                         else
                                {
                                 Fnn_r.p[i][j] += Fn_r.p[i][j];
                                 Fnn_i.p[i][j] += Fn_i.p[i][j];
                                }
                        }
         Fn_r = Fnn_r;
         Fn_i = Fnn_i;

         cout << "  Number of pixels updated during iteration " << n << " = " 
                << count << "/" << height*width << endl;

         if (count == 0)
                break;
        }

//Write2DTecplotFile(&Fn_r, &Fn_i, "out.tec");

double threshold = 0.0;

Edge->Zero();
for (int i = 1; i < height; i++)
        for (int j = 1; j < width; j++)
                {
                 float energy = 0.0;
                 if (Fn_r.p[i-1][j] > threshold && Fn_r.p[i][j] < -threshold)
                        energy += Fn_r.p[i-1][j] - Fn_r.p[i][j];
                 if (Fn_i.p[i][j-1] > threshold && Fn_i.p[i][j] < -threshold)
                        energy += Fn_i.p[i][j-1] - Fn_i.p[i][j];
                 Edge->p[i][j] = energy/2.0;
                }

cout << "Edge Detection ends!" << endl;
cout.flush();
}
//**************************************************************************
void Combine2Images(GrayImage &Orig, GrayImage &Edge, GrayImage *Final)
{
/*
GrayImage Edge2;
Edge2 = Edge;
Edge2.DilationNxN(3);
*/
for (int i = 0; i < Orig.height; i++)
        for (int j = 0; j < Orig.width; j++)
                {
                 if (Edge.p[i][j] == 255)
                        Final->p[i][j] = Edge.p[i][j];
                 else
                        Final->p[i][j] = Orig.p[i][j];
                }
}
//**************************************************************************