/*********************************************************************/
/* Program samedge1.c    */
/*********************************************************************/
/* This program accepts as input a float BIL HSI image cube.  It then
   calculates a new image cube composed of the magnitude of the
   spectral angle of each original image cube vector (spectrum)
   wrt every other pixel in the cube.  In other words, the first
   "band" of the output cube is composed of the spectral angle
   of each original image cube spectrum wrt the pixel at
   (1,1).  The 2nd "band" contains spectral angle of each image
   spectrum wrt the pixel at (1,2), &etc.

   This verison employs calloc-generated dynamic menory allocation.

   Most current version as of 30 May 2011; supercedes HySPADE.cpp.

   Program written by: Ronald G. Resmini, Ph.D.
                       May 30, 2011
*/
/*********************************************************************/
/* Modification History:                                             */
/*********************************************************************/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>

// Dimensions for the input HSI data cube:
#define ROWFULL  600 // Lines
#define COLFULL  320 // Samples
#define BANDS     30 // Bands
#define HEADER     0 // Header

// Dimensions for window
#define ROW    70
#define COL    70
#define RXC    4900
#define WINDOW_ROW_STEP		68
#define WINDOW_COL_STEP		68
#define EP     20
#define EPDIV   5

/* Subroutine Prototypes:  */
  void inputcube(float ***hycube, FILE *infile);
  void calcspecmag(float ***hycube, float **vectormag);
  void procspecX(float ***hycube, float **vectormag, float ***samcube);
  void procspecY(float ***hycube, float **vectormag, float ***samcube);
  void outputmat(float ***samcube, FILE *outfile1);
  void outputcub(float ***hycube, FILE *outfile3);
  void outputedges(short ***edges, FILE *outfile2);
  void findedgesX(float ***samcube, short ***edges);
  void findedgesY(float ***samcube, short ***edges);
  float calcstdev(float [RXC]);

int main() /* Begin main() program processing. */
{
  FILE  *infile,
//      *outfile1,
        *outfile2;
  char  infile_name[80];
//  char  outfile1_name[80];
  char  outfile2_name[80];
  short ***edges_full,
        ***edges_window;
  float ***hycube,  // full hsi cube
        ***window,  // N x N x(NxN)
        ***samcube,
         **vectormag;
  short i, j, m;

  printf("%s\n\n", "Beginning of program HySPADE");

  /* Build the dynamically created arrays:                     */

  hycube = (float***) calloc(ROWFULL, sizeof(float**));
  for (i = 0; i <= (ROWFULL-1); i++)
    {
     if ((hycube[i] = (float**) calloc(COLFULL, sizeof(float*))) == NULL)
	 {
		 printf("\n*** Error allocating hycube array ***\n");
		 exit(0);
	 }
    }
  for (i = 0; i <= (ROWFULL-1); i++)
    {
     for (j = 0; j <= (COLFULL-1); j++)
       {
        if ((hycube[i][j] = (float*) calloc(BANDS, sizeof(float))) == NULL)
		{
			printf("\n*** Error allocating hycube array ***\n");
			exit(0);
		}
       }
    }

  window = (float***) calloc(ROW, sizeof(float**));
  for (i = 0; i <= (ROW-1); i++)
    {
     if ((window[i] = (float**) calloc(COL, sizeof(float*))) == NULL)
	 {
		 printf("\n*** Error allocating hycube array ***\n");
		 exit(0);
	 }
    }
  for (i = 0; i <= (ROW-1); i++)
    {
     for (j = 0; j <= (COL-1); j++)
       {
        if ((window[i][j] = (float*) calloc(BANDS, sizeof(float))) == NULL)
		{
			printf("\n*** Error allocating hycube array ***\n");
			exit(0);
		}
       }
    }

  edges_window = (short***) calloc(ROW, sizeof(short**));
  for (i = 0; i <= (ROW-1); i++)
    {
     if ((edges_window[i] = (short**) calloc(COL, sizeof(short*))) == NULL)
	 {
		 printf("\n*** Error allocating edges_window array ***\n");
		 exit(0);
	 }
    }
  for (i = 0; i <= (ROW-1); i++)
    {
     for (j = 0; j <= (COL-1); j++)
       {
        if ((edges_window[i][j] = (short*) calloc(EP, sizeof(short))) == NULL)
		{
			printf("\n*** Error allocating edges_window array ***\n");
			exit(0);
		}
       }
    }

  edges_full = (short***) calloc(ROWFULL, sizeof(short**));
  for (i = 0; i <= (ROWFULL-1); i++)
    {
     if ((edges_full[i] = (short**) calloc(COLFULL, sizeof(short*))) == NULL)
	 {
		 printf("\n*** Error allocating edges_full array ***\n");
		 exit(0);
	 }
    }
  for (i = 0; i <= (ROWFULL-1); i++)
    {
     for (j = 0; j <= (COLFULL-1); j++)
       {
        if ((edges_full[i][j] = (short*) calloc(EP, sizeof(short))) == NULL)
		{
			printf("\n*** Error allocating edges_full array ***\n");
			exit(0);
		}
       }
    }

  samcube = (float***) calloc(ROW, sizeof(float**));
  for (i = 0; i <= (ROW-1); i++)
    {
     if ((samcube[i] = (float**) calloc(COL, sizeof(float*))) == NULL)
	 {
		 printf("\n*** Error allocating samcube array ***\n");
		 exit(0);
	 }
    }
  for (i = 0; i <= (ROW-1); i++)
    {
     for (j = 0; j <= (COL-1); j++)
       {
        if ((samcube[i][j] = (float*) calloc(RXC, sizeof(float))) == NULL)
		{
			printf("\n*** Error allocating samcube array ***\n");
			exit(0);
		}
       }
    }

  vectormag = (float**) calloc(ROW, sizeof(float*));
  for (i = 0; i <= (ROW-1); i++)
    {
     vectormag[i] = (float*) calloc(COL, sizeof(float));
    }

  /* Poll user for input information; e.g., file names, etc.:  */

  printf("%s", "Enter the input BIL cube file_name: ");
  scanf("%s", infile_name);
  printf("\n");
//  printf("%s", "Enter the output samedge cube file #1 name: ");
//  scanf("%s", outfile1_name);
//  printf("\n");
  printf("%s", "Enter the output edges file #2 name: ");
  scanf("%s", outfile2_name);
  printf("\n");

  /* Open the input and output files: */
  infile    = fopen(infile_name, "rb");
//  outfile1  = fopen(outfile1_name, "wb");
  outfile2  = fopen(outfile2_name, "wb");

  /* Call subroutine to input the BIL cube: */
  inputcube(hycube, infile);
  /*Close the input file: */
  fclose(infile);

unsigned lineStep = WINDOW_ROW_STEP;
unsigned sampStep = WINDOW_COL_STEP;
for(unsigned ulWindowLine = 0; ulWindowLine < ROWFULL-ROW+1; ulWindowLine += lineStep) {
	for(unsigned ulWindowSamp = 0; ulWindowSamp < COLFULL-COL+1; ulWindowSamp += sampStep) {
	
	printf("----> ulWindowLine=%d, ulWindowSamp=%d <-----\n", ulWindowLine, ulWindowSamp);

	// copy spectra into window:
	for(unsigned line = 0; line < ROW; ++line)
        {
	   for(unsigned samp = 0; samp < COL; ++samp)
           {
		 for(unsigned band = 0; band < BANDS; ++band)
               {
		     window[line][samp][band] = hycube[ulWindowLine+line][ulWindowSamp+samp][band];
		    }
		}
	  }

  /* To save on computation, the vector magnitudes are calculated...*/
  /* Calculate image cube vector (spectra) magnitudes: */
  calcspecmag(window, vectormag);

// Initialize edges_window[][][]:
for (i=0; i<=(ROW-1); i++)
   {
    for (j=0; j<=(COL-1); j++)
       {
        for (m=0; m<=(EP-1); m++)
           {
             edges_window[i][j][m] = 0;
           }
       }
    }

  /* Perform main spectral processing X direction: */
  procspecX(window, vectormag, samcube);

  /* Call subroutine to find the edges in the X direction: */
  findedgesX(samcube, edges_window);

  /* Write the x-direction samedge cube.           */
  /* Call subroutine to write the BSQ output file: */
//  outputmat(samcube, outfile1);
//  fclose(outfile1);
  
  /* Perform main spectral processing Y direction: */
  procspecY(window, vectormag, samcube);
  
  /* Call subroutine to find the edges in the Y direction: */
  findedgesY(samcube, edges_window);


	// copy window edges onto entire edge plane:
	// (starting at 1 and ending at -1 so as to avoid adding edge artifacts from the window outline!)
    for (m=0; m<=(EP-1); m++)
     {
	for(unsigned line = 1; line < ROW-1; ++line) {
		for(unsigned samp = 1; samp < COL-1; ++samp) {
			edges_full[ulWindowLine+line][ulWindowSamp+samp][m] += edges_window[line][samp][m];
		}
	}
    }

	}
}

  /* Call subroutine to write the output file: */
  outputedges(edges_full, outfile2);
  fclose(outfile2);

  printf("%s\n", "End of program samedge1.c");

  return 0;

} /* End of main() program. */

/*********************************************************************/
/* Function declarations: */
/*********************************************************************/
void inputcube(float ***hycube, FILE *infile)
 {
    int i, j, k;
    char junk;

    printf("%s\n\n", "Begin procedure inputcube. ");

/* Read in the header data: */
     for (i=1; i<=HEADER; i++)
          fread(&junk, sizeof(char), 1, infile);

/* Read in the BIL cube: */
  for (i=0; i<=(ROWFULL-1); i++)
    {
     for (j=0; j<=(BANDS-1); j++)
       {
        for (k=0; k<=(COLFULL-1); k++)
	  {
           fread(&hycube[i][k][j], sizeof(float), 1, infile);
          } /* End of for k. */
        } /* End for j. */
     if (i%2 == 0) printf("%s%d%s\n", "Line: ", (i+1), " read in.");
    } /* End for i. */

  printf("%s\n\n", "End of procedure inputcube. ");

  } /* End of subroutine inputcube. */
/*********************************************************************/
void calcspecmag(float ***hycube, float **vectormag)
  {
   int i, j, k;

   printf("%s\n\n", "Beginning of procedure calcspecmag. ");

   for (i=0; i<=(ROW-1); i++)
     {
      for (j=0; j<=(COL-1); j++)
        {
		 vectormag[i][j] = 0.0;
         for (k=0; k<=(BANDS-1); k++)
	   {
	    vectormag[i][j] += ((float) hycube[i][j][k] * (float) hycube[i][j][k]);
           } /* End of for k. */
/*	 printf("i=%d, j=%d, vectormag=%f\n", i, j, vectormag[i][j]);    */
	 vectormag[i][j] = sqrt(vectormag[i][j]);
        } /* End for j. */
     } /* End for i. */

   printf("%s\n\n", "End of procedure calcspecmag. ");

  } /* End of subroutine calcspecmag. */
/*********************************************************************/
void procspecX(float ***hycube, float **vectormag, float ***samcube)

  {
    int i, j, k, l, m, samband;
    float numerator;

    printf("%s\n\n", "Begin procedure procspecX. ");

    samband = 0;
    for (i=0; i<=(ROW-1); i++)
      {
       for (j=0; j<=(COL-1); j++)
         {
          for (k=0; k<=(ROW-1); k++)
	      {
             for (l=0; l<=(COL-1); l++)
               {
                numerator = 0.0;
	          for (m=0; m<=(BANDS-1); m++)
	            {
	             numerator += ((float) hycube[i][j][m] * (float) hycube[k][l][m]);
	            } /* End of for m */
if ((numerator / (vectormag[i][j] * vectormag[k][l])) > 1.000)
   {
     samcube[k][l][samband] = 0.00;
   }
else
   {
    samcube[k][l][samband] = acos(((numerator / (vectormag[i][j] * vectormag[k][l]))));
   }
               } /* End of for l */
            } /* End of for k */
/*		  printf("%s%d\n", "samband: ", samband);  */
          samband = samband + 1;
         } /* End of for j */
      } /* End of for i */

   printf("%s\n\n", "End of procedure procspecX. ");

  } /* End of subroutine procspecX. */
/*********************************************************************/
void procspecY(float ***hycube, float **vectormag, float ***samcube)

  {
    int i, j, k, l, m, samband;
    float numerator;

    printf("%s\n\n", "Begin procedure procspecY. ");

    samband = 0;
    for (i=0; i<=(COL-1); i++)
      {
       for (j=0; j<=(ROW-1); j++)
         {
          for (k=0; k<=(COL-1); k++)
	      {
             for (l=0; l<=(ROW-1); l++)
               {
                numerator = 0.0;
	          for (m=0; m<=(BANDS-1); m++)
	            {
	             numerator += ((float) hycube[j][i][m] * (float) hycube[l][k][m]);
	            } /* End of for m */
if ((numerator / (vectormag[j][i] * vectormag[l][k])) > 1.000)
   {
     samcube[l][k][samband] = 0.00;
   }
else
   {
    samcube[l][k][samband] = acos(((numerator / (vectormag[j][i] * vectormag[l][k]))));
   }
               } /* End of for l */
            } /* End of for k */
          samband = samband + 1;
         } /* End of for j */
      } /* End of for i */

   printf("%s\n\n", "End of procedure procspecY. ");

  } /* End of subroutine procspecY. */
/*********************************************************************/
void outputmat(float ***samcube, FILE *outfile)
{
   int i, j, k;

   printf("%s\n\n", "Begin procedure outputmat. ");

   /* Write out the BSQ cube: */
   for (i=0; i<=(RXC-1); i++)
     {
      for (j=0; j<=(ROW-1); j++)
        {
         for (k=0; k<=(COL-1); k++)
	   {
            fwrite(&samcube[j][k][i], sizeof(float), 1, outfile);
           } /* End of for k. */
        } /* End for j. */
 /*     printf ("Wrote row #%d/%d\n",i,ROW);    */
     } /* End for i. */

   printf("%s\n\n", "End of procedure outputmat. ");

  } /* End of subroutine outputmat. */

/*********************************************************************/
void outputedges(short ***edges, FILE *outfile)
{
   int i, j, m;

   printf("%s\n\n", "Begin procedure outputedges. ");

   /* Write out the edge-detected planes BSQ cube: */
   for (m=0; m<=(EP-1); m++)
     {
      for (i=0; i<=(ROWFULL-1); i++)
        {
         for (j=0; j<=(COLFULL-1); j++)
           {
            fwrite(&edges[i][j][m], sizeof(short), 1, outfile);
           } /* End for j. */
        } /* End for i. */
     } /* End of for m. */

   printf("%s\n\n", "End of procedure outputmat. ");

  } /* End of subroutine outputedges. */
/*********************************************************************/
void findedgesX(float ***samcube, short ***edges)
 {
  short i, j, k, m;
  short **tally, edgeX, edgeY;
  float *diff, tstdev, stdev;
  
  printf ("%s\n\n", "Begin procedure findedgesX. ");

  diff  = (float*) calloc(RXC, sizeof(float));
  tally = (short**) calloc(RXC, sizeof(short*));  //Now a 2D array
  for (i=0; i<=(RXC-1); i++)
    {
     tally[i] = (short*) calloc(EP, sizeof(short));
    }

  /* Initialize the tally and difference arrays to zero: */
  for (i=0; i<=(RXC-1); i++)
    {
     diff[i] = 0.0;
    } /* End of for i */

  for (i=0; i<=(RXC-1); i++)
    {
      for (m=0; m<=(EP-1); m++)
        {
         tally[i][m] = 0;
        }
    } /* End of for i */

  for (i=0; i<=(ROW-1); i++)
    {
     for (j=0; j<=(COL-1); j++)
       {
        diff[0] = samcube[i][j][1] - samcube[i][j][0];
        for (k=1; k<=(RXC-1); k++)
          {
           diff[k] = samcube[i][j][k] - samcube[i][j][k-1];
          } /* End of for k */
           stdev = calcstdev(diff);
           for (m=0; m<=(EP-1); m++)
             {
              tstdev = stdev * ((float)(m+1)) / EPDIV;
              for (k=0; k<=(RXC-1); k++)
                {
                 if ((diff[k] >= tstdev) || (diff[k] <= -tstdev))
                   {
                    tally[k][m] = tally[k][m] + 1;
                   }
                } /* End of for k */
             } /* End of for m */
        } /* End of for j */
       } /* End of for i */

  for (m=0; m<=(EP-1); m++)
    {
     for (k=0; k<=(RXC-1); k++)
       {
        if (((k+1) % COL) == 0)
         {
          edgeY = COL;
         }
        else
         {
          edgeY = (k+1) % COL;
         } /* End if */
        edgeX = (((k+1) + (COL - edgeY)) / COL);
        edges[edgeX-1][edgeY-1][m] += tally[k][m];
       } /* End of for k */
    } /* End of for m */

  printf ("%s\n\n", "End procedure findedgesX. ");

 } /* End of subroutine findedgesX. */
/*********************************************************************/
void findedgesY(float ***samcube, short ***edges)
 {
  short i, j, k, m;
  short **tally, edgeX, edgeY;
  float *diff, tstdev, stdev;
  
  printf ("%s\n\n", "Begin procedure findedgesY. ");
     diff  = (float*) calloc(RXC, sizeof(float));
     tally = (short**) calloc(RXC, sizeof(short*));  //Now a 2D array
     for (i=0; i<=(RXC-1); i++)
       {
        tally[i] = (short*) calloc(EP, sizeof(short));
       }

  /* Initialize the tally and difference arrays to zero: */
  for (i=0; i<=(RXC-1); i++)
    {
     diff[i] = 0.0;
    } /* End of for i */

     for (i=0; i<=(RXC-1); i++)
       {
         for (m=0; m<=(EP-1); m++)
           {
            tally[i][m] = 0;
           }
       } /* End of for i */

  for (i=0; i<=(ROW-1); i++)
    {
     for (j=0; j<=(COL-1); j++)
       {
        diff[0] = samcube[i][j][1] - samcube[i][j][0];
        for (k=1; k<=(RXC-1); k++)
          {
           diff[k] = samcube[i][j][k] - samcube[i][j][k-1];
          } /* End of for k */
           stdev = calcstdev(diff);
           for(m=0; m<=(EP-1); m++)
              {
               tstdev = stdev * ((float)(m+1)) / EPDIV;
               for (k=0; k<=(RXC-1); k++)
                 {
                  if ((diff[k] >= tstdev) || (diff[k] <= -tstdev))
                    {
                     tally[k][m] = tally[k][m] + 1;
                    }
                 } /* End of for k */
              } /* End of for m */
     } /* End of for j */
    } /* End of for i */

  for(m=0; m<=(EP-1); m++)
    {
     for (k=0; k<=(RXC-1); k++)
       {
        if (((k+1) % ROW) == 0)
         {
          edgeX = ROW;
         }
        else
         {
          edgeX = (k+1) % ROW;
         } /* End if */
        edgeY = (((k+1) + (ROW - edgeX)) / ROW);
        edges[edgeX-1][edgeY-1][m] += tally[k][m];
       } /* End of for k */
    } /* End of for m */

  printf ("%s\n\n", "End procedure findedgesY. ");

 } /* End of subroutine findedgesY. */
/*********************************************************************/
float calcstdev(float diff[RXC])
 {
  int i;
  float sum, sumSq, stdev;

  //printf ("%s\n\n", "Begin function calc_stdev. ");

  sum = 0.0;
  sumSq = 0.0;

  for (i=0; i<=(RXC-1); i++)
    {
     sum = sum + diff[i];
     sumSq = sumSq + (diff[i] * diff[i]);
    } /* End of for i */

  stdev = sqrt( (sumSq - ((sum * sum) / RXC)) / (RXC-1) );

  //printf ("%s\n\n", "End function calc_stdev. ");

  return stdev;

 } /* End of function calc_stdev. */

/*********************************************************************/
void outputcub(float ***hycube, FILE *outfile)
{
   int i, j, k;

   printf("%s\n\n", "Begin procedure outputcub. ");

   /* Write out the BIL cube: */
   for (i=0; i<=(ROW-1); i++)
     {
      for (j=0; j<=(BANDS-1); j++)
        {
         for (k=0; k<=(COL-1); k++)
	   {
            fwrite(&hycube[i][k][j], sizeof(float), 1, outfile);
           } /* End of for k. */
        } /* End for j. */
     } /* End for i. */

   printf("%s\n\n", "End of procedure outputcub. ");

  } /* End of subroutine outputcub. */
/*********************************************************************/
/* End of function declarations.  End of program listing.  */
/*********************************************************************/