/*	ISODATA.c
**	this programm clusters the
**      generic binary image file
**      into user defined cluster
**      number based on hard c mean!
**      The programm is used for 3 band
**      images. If you want to cluster
**      more bands, the programm has to be
**      slightly modified!
*/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <time.h>

#define num_of_band 3
#define cluster_num 4


main(int argc, char **argv)
{
  FILE 	 *fin1,*fin2,*fin3,*fin4,*fin5,*fin6,*fin7,*fin8,*fin9,*fin10;
  FILE 	 
*fout,*fout1,*fout2,*fout3,*fout4,*fout5,*fout6,*fout7,*fout8,*fout9,*fout10;
  unsigned long int img_length, sum_of_uik;
  unsigned long int img_width, total_pix, step;
  unsigned long int count, count1, count2, count3;
  unsigned char 
out1[1],i1[1],i2[1],i3[1],i4[1],i5[1],i6[1],i7[1],i8[1],i9[1],i10[1];
  long float  v1[10][10],v2[10][10],dik[10],po[10],min,pow_num1,pow_num2;
  int weight;

  if(argc == 7)
    {
      img_length = atoi(argv[2]);
      img_width = atoi(argv[3]);

      if((fin1 = fopen(argv[4],"rb")) == NULL)               /*image 1*/
	{
	  fprintf(stderr,"Failed to open %s\n",argv[4]);
	  exit(EXIT_FAILURE);
	}
      if((fin2 = fopen(argv[5],"rb")) == NULL)               /*image 2*/
	{
	  fprintf(stderr,"Failed to open %s\n",argv[5]);
	  exit(EXIT_FAILURE);
	}
      if((fin3 = fopen(argv[6],"rb")) == NULL)               /*image 3*/
	{
	  fprintf(stderr,"Failed to open %s\n",argv[6]);
	  exit(EXIT_FAILURE);
	}
      if((fout = fopen(argv[1],"w+b")) == NULL)              /*output 
image*/
       {
         fprintf(stderr,"Failed to open %s\n",argv[1]);
         exit(EXIT_FAILURE);
       }

      printf("\n");
      printf("***Initializing Classification Matrix U! Please wait!***\n");
      initial_U(fout,img_width,img_length);
      printf("Initialization Finish!\n");
      total_pix = img_length*img_width;

      fseek(fout,0,0);

      for(count1 = 1;count1 <= cluster_num;count1++)/*v[cluster][band]*/
        for(count2 = 1;count2 <= num_of_band;count2++)
           v1[count1][count2] = 0;

      printf("Start Calculating!\n");
      count = 1;

      for(count = 0; count < 20; count++)/*Iterations chosen as 20 is 
sufficient for convergence*/
      {
      for(count1 = 1;count1 <= cluster_num;count1++)/*calculate mean for 
each cluster*/
	{
        sum_of_uik = 0;
        v2[count1][1] =  v1[count1][1];
        v2[count1][2] =  v1[count1][2];
        v2[count1][3] =  v1[count1][3];
        for(count2 = 1;count2 <= num_of_band;count2++)
           v1[count1][count2] = 0;

        for(count2 = 1;count2 <= total_pix;count2++)
	  {
	   fread(out1,1,1,fout);
           fread(i1,1,1,fin1);
           fread(i2,1,1,fin2);
           fread(i3,1,1,fin3);
	   if(out1[0] == count1)
	     {
	      sum_of_uik++;
              v1[count1][1] = v1[count1][1] + i1[0];
              v1[count1][2] = v1[count1][2] + i2[0];
              v1[count1][3] = v1[count1][3] + i3[0];
	     }
	  }
        printf("number of cluster %d : %d\n",count1,sum_of_uik);
        v1[count1][1] = (float) v1[count1][1]/sum_of_uik;
        v1[count1][2] = (float) v1[count1][2]/sum_of_uik;
        v1[count1][3] = (float) v1[count1][3]/sum_of_uik;
        printf("mean of cluster %d : (%lf,%lf)\n",count1, v1[count1][1], 
v1[count1][2], v1[count1][3]);
        fseek(fin1,0,0);
        fseek(fin2,0,0);
        fseek(fin3,0,0);
        fseek(fout,0,0);
	}

      pow_num1 = 2;
      pow_num2 = 0.5;
      for(count1 = 1;count1 <= total_pix;count1++)/*update U matrix based on 
the nearest distance*/
	{
         fread(i1,1,1,fin1);
         fread(i2,1,1,fin2);
         for(count2 = 1;count2 <= cluster_num;count2++)
	   {
           po[1]= powf(i1[0]-v1[count2][1],pow_num1);
           po[2]= powf(i2[0]-v1[count2][2],pow_num1);
           po[3]= powf(i3[0]-v1[count2][3],pow_num1);
           dik[count2] = powf(po[1]+po[2]+po[3],pow_num2); /*distance*/
	   }
         min = dik[1];
         weight = 1;
         for(count2 = 1;count2 <= cluster_num;count2++)
	   {
            if(dik[count2] < min)
               weight = count2;
	   }
         fprintf(fout,"%c",weight);
	}
        fseek(fin1,0,0);
        fseek(fin2,0,0);
        fseek(fin3,0,0);
        fseek(fout,0,0);
      }

      fclose(fin1);
      fclose(fin2);
      fclose(fin3);
      fclose(fout);


    }
  else
      printf("\nUsage:\n\tISODATA 'output image file'  'rows in file' ' 
cols' 'image 1' 'image 2' 'image 3'.\n\n");
}



initial_U(c1,width,length)
FILE *c1;
unsigned long int width,length;
{
unsigned long int total_pix, num_per_c;
unsigned long int count1, count2, count3;
int anxi_x, anxi_y;
int weight;
char a[1];

total_pix = length*width;
num_per_c = total_pix/cluster_num;

for(count1 = 0;count1 < length;count1++)
   for(count2 = 0;count2 < width;count2++)
     fprintf(c1,"%c",1);
printf("***Step 1 Finished!\n");
printf("***Start up Random Clustering\n");
fseek(c1,0,0);
for(count1 = 1;count1 < cluster_num;count1++)
   {
   weight = count1+1;
   for(count2=0;count2 < num_per_c;count2++)
     {
      anxi_x = random()%width;
      anxi_y = random()%length;
      if(anxi_x < 0)
         anxi_x = -anxi_x;
      if(anxi_y < 0)
         anxi_y = -anxi_y;
      for(count3 = 0;count3 < anxi_y;count3++)
         fseek(c1,width,1);
      fseek(c1,anxi_x,1);
      fprintf(c1,"%c",weight);
      fseek(c1,0,0);
     }
   printf("***Step %d Finished!\n",count1+1);
   }

}