#include <iostream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define ACROSS 204
#define DOWN 382

const int ASCPGM = 2,  ASCPPM = 3,   // "MAGIC NUMBERS"
          BINPGM = 5,  BINPPM = 6;

int  getPic(int [DOWN][ACROSS], char[]);       // get either, known dim
void putPic(int [DOWN][ACROSS], char[], int);  // put specified, known dim

int  sobel(int [DOWN][ACROSS], int [DOWN][ACROSS]);
void marrH(int [DOWN][ACROSS], int [DOWN][ACROSS]);
//void canny(int [DOWN][ACROSS], int [DOWN][ACROSS]);


int main(void)
{
   int pic[DOWN][ACROSS], // row-column form, as opposed to x-y
       edg[DOWN][ACROSS],
       mN;

   // input file
   mN = getPic(pic, "sodacan.pgm");

   // sobel edge detection with primitive "on-the-fly" threshhold
// NOTE:  "variable" threshholding POORLY implemented--REU students should try harder...   sobel(pic, edg);
   putPic(edg, "sobel.pgm", mN);

   // marr-hildreth edge-detection, based on a theory of how the eye works
   marrH(pic, edg);
   putPic(edg, "marrh.pgm", mN);

   // canny edge detection
   // REU students to implement

   return 0;
}

int getPic(int pic[DOWN][ACROSS], char fn[])
{
   ifstream from;
   FILE *p5;      // for binarys

   char bil[50], *snafu, c;
   bool hasComment = false;
   int lines = 3, mN, across, down, size;

   from.open(fn);
   if (from.fail()) {cout << "failed to open " << fn << "\n\a"; exit(-1); }

   from >> c >> mN;

   from >> bil;  snafu = bil;

   if (snafu[0] == '#')
   {
	   from.getline(bil, 50);
	   from >> across;
       hasComment = true;
       lines = 4;
   }
   else
   {
	   across = (int) atof(snafu);
   }

   from >> down >> size;

   if (mN == 5)	     // binary
   {
      from.close();

      p5 = fopen(fn, "rb");
      for(int x = 0; x < lines; x++)
         fgets(bil, 50, p5);

      for (int y = 0; y < down; y++)            // each row
         for (int x = 0; x < across; x++)         // across each row
         {
            pic[y][x] = (int) getc(p5); // [row][column]
         }

      fclose(p5);
   }

   else if (mN == 2)   // ASCII
   {
      for (int y = 0; y < down; y++)            // each row
         for (int x = 0; x < across; x++)         // across each row
             from >> pic[y][x];
       from.close();
    }

   else
   {
      from.close();
	  cout << "I can't read this kind of file yet!! \n\n";
	  exit(-1);  system("PAUSE");
   }

   return mN;
} // getPic

void putPic(int pic[DOWN][ACROSS], char fn[], int mN)
{
   ofstream to;
   FILE *p5;

   to.open(fn);
   if (to.fail()) {cout << "failed to open " << fn << "\n\a"; exit(-1); }

   to << 'P' << mN << '\n' << ACROSS << ' ' << DOWN << "\n255\n";

   if (mN == 5)	     // binary pgm
   {
      to.close();
      p5 = fopen(fn, "wb");

      fprintf(p5, "P5\n%d %d\n255\n", ACROSS, DOWN);

      for (int y = 0; y < DOWN; y++)          // each row
         for (int x = 0; x < ACROSS; x++)       // across each row
            fprintf(p5, "%c", pic[y][x]);

      fclose(p5);
   }

   else if (mN == 2)   // ASCII pgm
   {
      for (int y = 0; y < DOWN; y++)          // each row
         for (int x = 0; x < ACROSS; x++)       // across each row
             to << pic[y][x] << ' ';
      to.close();
   }
} // putPic

int sobel(int pic[DOWN][ACROSS], int after[DOWN][ACROSS])
// NOTE:  "variable" threshholding POORLY implemented--REU students should try harder...
{
   int tier[8] = {51, 76, 102, 127, 153, 178, 204, 229}, tiers = 8,    // 50%, 60%, ... 90%
       instances[8], distrib[8], maxDistrib, total, th;  // th = threshhold

   int maskX[3][3] = { {-1,-2,-1}, {0,0,0}, {1,2,1} },
       maskY[3][3] = { {1,0,-1}, {2,0,-2}, {1,0,-1} },
       tempX, tempY, max = 0,
       mr = 1, sum1 = 0, sum2 = 0,  //  mask radius, 2 counters in convolution
       scaled;  // scale the magnitude to fit 0-255

   int x, y, i, j, p, q;  // counting/looping utilities

   for (x = 0; x < ACROSS; x++)
      for (y = 0; y < DOWN; y++)
			after[y][x] = 0;

	for (i = 0; i < 8; i++)
		instances[i] = 0;
	total = 0;  maxDistrib = 0;  th = 0;


   for (i=mr;i<ACROSS-mr;i++)
   { for (j=mr;j<DOWN-mr;j++)
     {
        sum1 = 0;
        sum2 = 0;
        for (p=-mr;p<=mr;p++)
        {
           for (q=-mr;q<=mr;q++)
           {
              sum1 += pic[j+q][i+p] * maskX[q+mr][p+mr];
              sum2 += pic[j+q][i+p] * maskY[q+mr][p+mr];
           }
        }
        tempX = sum1;
        tempY = sum2;

        int root = (int) sqrt(  pow(tempX, 2) + pow(tempY, 2)  );
			 
        if (root > max) max = root;
             
        after[j][i] = root;
     }
   }

   for (x = 0; x < ACROSS; x++)
      for (y = 0; y < DOWN; y++)
      {
         int scaled = (int) (  (255.0/max) * after[y][x]  );

			// tier[8] = {51, 76, 102, 127, 153, 178, 204, 229}, tiers = 8
         if      (scaled > 229)   instances[7]++;    // 90%    // tier[7]
		 else if (scaled > 204)   instances[6]++;    // 80%    // tier[6]
		 else if (scaled > 178)   instances[5]++;    // 70%    // tier[5]
         else if (scaled > 153)   instances[4]++;    // 60%    // tier[4]
         else if (scaled > 127)   instances[3]++;    // 50%    // tier[3]
         else if (scaled > 102)   instances[2]++;    // 40%    // tier[2]
         else if (scaled >  76)   instances[1]++;    // 30%    // tier[1]
         else if (scaled >  51)   instances[0]++;    // 20%    // tier[0]

         after[y][x] = scaled;
      }
   for (i = 0; i < 8; i++)
	{
		distrib[i] = (int) (( (double) instances[i] / (ACROSS * DOWN) ) * 100);
		if (distrib[i] > maxDistrib) maxDistrib = distrib[i];
	}

	i = 7;
	while (th == 0)
	{
		total += distrib[i];
		if   (total > 1.5 * maxDistrib)   th = tier[i];
		else                              i--;
	}

   for (y = 0; y < DOWN; y++)
      for (x = 0; x < ACROSS; x++)
      {
         if( after[y][x] > th ) after[y][x] = 0;
         else                   after[y][x] = 255;
      }

   return th;   // variable threshhold
} // sobel

void marrH(int pic[DOWN][ACROSS], int after[DOWN][ACROSS])
{
   int    temp[DOWN][ACROSS], edgeflag[DOWN][ACROSS],
          i,j,p,q,x,y,
          center = 10, mr;

   double mask[20][20], maskval,
          sigma,
          sum, max, lmax, lmin, scaled,
          TOLERZ, TOLERD;     // how close to 0.0 IS 0.0, how great is the Diff

   sigma = 1;
   mr = (int) (3 * sigma);
   TOLERZ = 6;   TOLERD = 100;
   max = 0; lmax = 0; lmin = 255;

   for(j = -mr; j <= mr; j++)       // j ~corr. to y, row
   {
      for(i = -mr; i <= mr; i++)    // i ~corr. to x, column
      {
         maskval = ((2-(((i*i)+(j*j))/(sigma*sigma)))*
                      (exp(-1.0*(((i*i)+(j*j))/(2*(sigma*sigma))))));
         mask[center + j][center + i] = maskval;
         
      } 
   }

   for (j=mr;j<DOWN-mr;j++)
   {
     for (i=mr;i<ACROSS-mr;i++)
     {
        sum = 0;
        for (q=-mr;q<=mr;q++)
        {
           for (p=-mr;p<=mr;p++)
           {
              sum += pic[j+q][i+p] * mask[q + center][p + center] ;
           }
        }
        temp[j][i] = (int) sum;

        if (fabs(sum) > max)  max = fabs(sum);
     }
   }
   
   for (y = 0; y < DOWN; y++)
      for (x = 0; x < ACROSS; x++)
      {
         scaled = (  (temp[y][x] / max) + 1  ) * 127;
         after[y][x] = (int) scaled;

         if (   temp[y][x] < TOLERZ  &&  temp[y][x] > -TOLERZ   )
            temp[y][x] = 0;
      }

   for (y = 1; y < DOWN-1; y++)
      for (x = 1; x < ACROSS-1; x++)
      {
         if (  temp[y][x] < TOLERZ  &&  temp[y][x] > -TOLERZ  )   // no  == 0.0
         {
            // remember, all zero's must be zero'd before this
            //        one pos, one neg  &&  diff btx > TOLERD
            if (   (  temp[y+1][x-1] * temp[y-1][x+1]  < 0
                         &&   fabs( fabs(temp[y+1][x-1]) - fabs(temp[y-1][x+1]) ) > TOLERD  )   ||
                   (  temp[y][x-1]   * temp[y][x+1]    < 0
                         &&   fabs( fabs(temp[y][x-1]) - fabs(temp[y][x+1]) ) > TOLERD  )       ||
                   (  temp[y-1][x-1] * temp[y+1][x+1]  < 0
                         &&   fabs( fabs(temp[y-1][x-1]) - fabs(temp[y+1][x+1]) ) > TOLERD  )   ||
                   (  temp[y+1][x]   * temp[y-1][x]    < 0
                         &&   fabs( fabs(temp[y+1][x]) - fabs(temp[y-1][x]) ) > TOLERD  )          )
               after[y][x] = 0;
            else  after[y][x] = 255;
         }
         else if (temp[y][x] > TOLERD)
         {
            for ( j = -1; j <= 2; j++)
               for ( i = -1; i <= 2; i++)
               {
                  if (temp[y+j][x+i] < -TOLERD)
                  {  after[y][x] = 0;  edgeflag[y][x] = 1;   }
                  else if (edgeflag[y][x] != 1)
                     after[y][x] = 255;
               }
         }
         else  after[y][x] = 255;
      }  // end inner loop (x)
   // end outer loop (y)   
   
} // marrH