/***************************************************************************
 * Func: homog_op                                                          *
 *                                                                         *
 * Desc: detects edges in an image using the homogeneity operator          *
 *                                                                         *
 * Params: source - pointer to image in memory                             *
 *         cols - number of columns in image                               *
 *         rows - number of rows in image                                  *
 *         filename - name of output file                                  *
 ***************************************************************************/

void homog_op(image_ptr source, int cols, int rows, char *filename)
    {
    int x, y, i;             /* image loop variables */
    int index;               /* image index */
    int winx, winy;          /* indices for filling sample window */
    int xextra, yextra;      /* size of boundary */
    int conv_line;           /* size of output line */
    unsigned long destadd;   /* destination image address */
    unsigned long sourceadd, sourcebase; /* source addressing */
    char dest[1024];         /* destination image line */
    FILE *fp;                /* output file pointer */
    unsigned char window[9]; /* storage for window samples */
    unsigned char left[25];  /* storage for left edge pixel duplication */
    unsigned char right[25]; /* storage for right edge pixel duplication */
    int xpad, ypad;          /* number of pixels to duplicate */
    int last_line;           /* last line number to process */
    int new_pixel;           /* newly computed pixel value */

    yextra = 2;
    ypad = 1;
    xextra = 2;
    xpad = 1;
    conv_line = cols - xextra;

    if((fp=fopen(filename, "wb")) == NULL)
	{
	printf("Unable to open %s for output\n",filename);
	exit(1);
	}
    fprintf(fp, "P5\n%d %d\n255\n", cols, rows); /* print out header */
    last_line = rows - yextra;
    for(y=0; y<last_line; y++)
	{
	sourcebase=(unsigned long) cols * y;
	destadd=0;
	for(x=xpad; x<(cols - xpad); x++)
	   {
	   index=0;
	   for(winy=0; winy<3; winy++)
	       for(winx=0; winx<3; winx++)
		  {
		  sourceadd=sourcebase+winx+winy*cols;
		  window[index++] = source[sourceadd];
		  }
	   new_pixel = homogeneity(window) + 64;
	   CLIP(new_pixel,0.0,255.0);

	   dest[destadd++]=(unsigned char) new_pixel;
	   sourcebase++;
	   } /* for x */
	 for(i=0; i<xpad; i++)
	   left[i]=dest[0];
	 for(i=0; i<xpad; i++)
	   right[i]=dest[conv_line-1];
	 if(y==0)
	     for(i=0; i<ypad; i++)
	     {
	     fwrite(left, 1, xpad, fp);
	     fwrite(dest, 1, conv_line, fp);
	     fwrite(right, 1, xpad, fp);
	     }
	 fwrite(left, 1, xpad, fp);
	 fwrite(dest, 1, conv_line, fp);
	 fwrite(right, 1, xpad, fp);
	 if(y==(last_line-1))
	     for(i=0; i<ypad; i++)
		 {
		 fwrite(left, 1, xpad, fp);
		 fwrite(dest, 1, conv_line, fp);
		 fwrite(right, 1, xpad, fp);
		 }
    } /* for y */

    }

/****************************************************************************
 * Func: homogeneity                                                        *
 *                                                                          *
 * Desc: finds edges in an image window                                     *
 *                                                                          *
 * Param: window - 3x3 window of source image   |0 1 2|                     *
 *                                              |3 4 5|                     *
 *                                              |6 7 8|                     *
 *                                                                          *
 * Returns: the maximum value of the difference between the center pixel and*
 *          all adjacent pixels                                             *
 ****************************************************************************/

int homogeneity(unsigned char *window)
    {
    int max;                 /* current maximum difference */
    int diff;                /* current difference */

    diff = abs(window[4] - window[0]);  /* upper left pixel */
    max = diff;

    diff = abs(window[4] - window[1]);  /* upper pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[2]);  /* upper right pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[5]);  /* right pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[8]);  /* lower left pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[7]);  /* lower pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[6]);  /* bottom left pixel */
    if(diff > max)
	max = diff;

    diff = abs(window[4] - window[3]);  /* left pixel */
    if(diff > max)
	max = diff;

    return max;
    }