/*
 * w4735 assignment 3 step 3 part 1
 * by anders pearson <anders@columbia.edu>
 * 2001-04-19
 */

#include "defHelp.h"
#include <stdlib.h>

#define MAXCOL 160

main(int argc, char *args[]){

  int S = 0;
  struct image img[9]; 
  int i, j, k,r,c,h;
  int w = 0;
  int row, col, sum, samples;
  float speckles[9];
  float diffs[9][9];
  float min = 1000.0;
  float max = 0.0;
  int mini,minj,maxi,maxj;

  float total,trans,mean,stddev;
  float prev;
  char *images[9] = {"iris1.pgm","iris2.pgm","iris3.pgm","iris4.pgm",
		    "iris5.pgm","iris6.pgm","iris7.pgm","iris8.pgm",
		    "iris9.pgm"};
  
  /* read in images */
  for(i = 0; i < 9; i++){
    if(readImage(&img[i], images[i]) != 0) {
      exit(1);
    }
  }
  
  for (S = 3; S <= 15; S += 2) { 
    for(i = 0; i < 9; i++) {
      h = getNRows(&img[i]);
      for(col = (S/2); col < (MAXCOL - (S/2)); col++) {
	for(row = (S/2); row < (h - (S/2)); row++){
	  /* calculate mean for area */
	  total = 0.0;
	  for(c = (col - (S/2)); c < (col + (S/2)); c++) {
	    for(r = (row - (S/2)); r < (row + (S/2)); r++) {
	      total += getPixel(&img[i],r,c);
	    }
	  }
	  mean = total/(S*S);
	  total = 0.0;
	  /* calculate stddev for area */
	  for(c = (col - (S/2)); c < (col + (S/2)); c++) {
	    for(r = (row - (S/2)); r < (row + (S/2)); r++) {
	      total += (getPixel(&img[i],r,c) - mean) * (getPixel(&img[i],r,c) - mean);
	    }
	  }
	  stddev = total/(S*S);
	  stddev = sqrt(stddev);
	  /* decide if it's a speckle */
	  if(getPixel(&img[i],row,col) > (mean + stddev)){
	    speckles[i]++;
	  }
	}
      }
      /* normalize against number of pixels sampled */
      speckles[i] /= ((MAXCOL - S) * (h - S));
    }
    /* generate difference matrix */
    for(i = 0; i < 9; i++) {
      for(j = i+1; j < 9; j++) {
	diffs[i][j] = speckles[i] - speckles[j];
	diffs[i][j] = (diffs[i][j] > 0) ? diffs[i][j] : -diffs[i][j];
      }
    }

    /* find minimum difference */
    min = diffs[0][1];
    max = diffs[0][1];
    mini = 0;
    minj = 1;
    maxi = 0;
    maxj = 1;
    for(i = 0; i < 9; i++) {
      for(j = i+1; j < 9; j++) {
	if(diffs[i][j] < min){
	  min = diffs[i][j];
	  mini = i;
	  minj = j;
	}
	if(diffs[i][j] > max) {
	  max = diffs[i][j];
	  maxi = i;
	  maxj = j;
	}
      }
    }

    if(min > 0.0 && max > 0.0) {
      printf("S: %d, sim (%s,%s,%f), ",S,images[mini],images[minj],min);
      printf("diff (%s,%s,%f)\n",images[maxi],images[maxj],max);
    }
  }
  return 0;
}