/*	fuzzy_c_means.c
**      This programme is used to for clustering images (unsupervised classifciation)
**      and output three files:
**      (1)A classified image.
**      (2)A fuzzy memership matrix called U_C matrix records fuzzy
**         membership for each pixel(see below); and
**      (2)A Mean_record stores the mean and pixel number for each cluster.                            pixels in each cluster!
**
**********************************************************************
**Fuzzy Membership Matrix U_C:
**===================================================================
**         cluster1 (cover1)|cluster2 (cover2)|cluster3 (cover3)|...
**pixel1  |    0.898        |  0.003          |  0.112  ...     |
**pixel2  |    0.100        |  0.947          |  0.825  ...     |
**pixel3  |    0.002        |  0.050          |  0.063  ...     |
**   .                      |            .    |            .    |
**   .                      |            .    |            .    |
**   .                      |            .    |            .    |
**===================================================================
*/

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


main()
{
  FILE 		*fin[21], *fout[10];
  unsigned int num_of_img, img_row;
  unsigned int img_col, total_pix;
  unsigned int user_num, pix_num, num_of_cycle;
  unsigned int step, pix, cluster_num, num_in_cover[40];
  unsigned int count1, count2, count3, count4;
  unsigned char i[21][5];
  unsigned char in_file[50], out_file[50];
  double out1, m;
  double v1[40][21],v2[40][21],dik[40],po[40],pow_num1,pow_num2;
  double sum_of_uik, po_sum;
  double weight,sum_dijk,class_type[40], max;


      printf("\n====Images clustering using fuzzy c means algorithm====\n");
	  printf("Please input the number of total images (maximum 20):");
      scanf("%d",&num_of_img);
	  for(count1=1; count1<=num_of_img; count1++)
		{
         printf("Input image file %d:", count1);
	     scanf("%s",in_file);
      	 if((fin[count1] = fopen(in_file,"rb")) == NULL)
			{
	         printf("Failed to open %s\n",in_file);
	         exit(EXIT_FAILURE);
			}
		}

	  printf("Please input the number for image rows:");
      scanf("%d",&img_row);
	  printf("Please input the number for image columns:");
      scanf("%d",&img_col);
	  total_pix = img_row*img_col;
	  printf("Input files completed!\n");

	  printf("Ouput clustered image file:");
	  scanf("%s",out_file);
	  if((fout[1] = fopen(out_file,"wb")) == NULL)
	   {
		printf("Failed to open %s\n",out_file);
		exit(EXIT_FAILURE);
		}
      if((fout[2] = fopen("U_C","w+")) == NULL)     /*Membership matrix*/
       {
        printf("Failed to open %s\n","U_C");
		exit(EXIT_FAILURE);
       }
      if((fout[3] = fopen("Mean_record","w")) == NULL)       /*Mean fo each cluster*/
       {
         fprintf(stderr,"Failed to open %s\n","Mean_record");
         exit(EXIT_FAILURE);
       }

      printf("Please input fuzziness m value (suggested range 1.5 to 2.5):");
      scanf("%lf",&m);
      printf("Please input random number seed:");
      scanf("%d",&user_num);
      printf("Please input the number of cluster to be formed:");
      scanf("%d",&cluster_num);
      printf("Please input the number for clustering iteration: ");
      scanf("%d",&num_of_cycle);

      printf("\n***Initializing means for each cluster! Please wait!***\n");

      for(count1 = 1;count1 <= cluster_num;count1++)
	  {
        pix_num = rand()+user_num*count1;
        if(pix_num > total_pix)
           pix_num = pix_num%total_pix;
        printf("\nRandomly located initial pixel position %d for cluster number %d",pix_num, count1);
        for(count2 = 1; count2 <= num_of_img; count2++)
			{
		  	  fseek(fin[count2],pix_num,0);
              fread(i[count2],1,1,fin[count2]);
              v1[count1][count2] =(double) i[count2][0] +(double) count1/8;
             }
       }

      for(count2 = 1; count2 <= num_of_img; count2++)
	  	  fseek(fin[count2],0,0);
      printf("\nInitialization finished! Start calculating!\n");

      step = cluster_num - 1;

      for(count1=1; count1 <= cluster_num;count1++)
		{
		printf("\nThe mean for cluster %d:",count1);
		for(count2=1; count2<= num_of_img; count2++)
		   printf("\n%lf",v1[count1][count2]);
		}

      pow_num1 = 2;
      pow_num2 = 0.5;

      for(count1 = 0;count1 < num_of_cycle;count1++)
		 {
		printf("\nCluster Cycle %d\n",count1);
		for(count2 = 1;count2 <= total_pix;count2++)
			{
		   for(count3 = 1; count3 <= cluster_num; count3++)
				{
			  po_sum = 0;
			  for(count4 = 1;count4 <= num_of_img;count4++)
					{
				     if(count3 == 1)
				        fread(i[count4],1,1,fin[count4]);
					 po[count4]= pow(i[count4][0]-v1[count3][count4],pow_num1);
					 po_sum = po_sum + po[count4];
				    }
			   dik[count3] =(double) pow(po_sum,pow_num2);
				}

           for(count3 = 1;count3 <= cluster_num;count3++)
				{
				sum_dijk = 0;
			    for(count4 = 1;count4 <= cluster_num;count4++)
					{
					sum_dijk = sum_dijk + pow(dik[count3]/dik[count4],pow_num1/(m-1));
					}
				weight = (double) 1/sum_dijk;
				fprintf(fout[2],"%lf ",weight);
				}
			}//end of count2 (total_pix)//

        for(count2 = 1;count2 <= num_of_img;count2++)
               fseek(fin[count2],0,0);
        fseek(fout[2],0,0);

		for(count2 = 1;count2 <= cluster_num;count2++)
			{
			for(count3 = 0;count3 < count2 - 1; count3++)
				fscanf(fout[2],"%lf ",&out1);
			for(count3 = 1; count3<= num_of_img; count3++)
				{
		         v2[count2][count3] =  v1[count2][count3];
				 v1[count2][count3] = 0;
				}
			sum_of_uik = 0;
			for(count3 = 1;count3 <= total_pix;count3++)
				{
				for(count4 = 1; count4<=num_of_img; count4++)
					fread(i[count4],1,1,fin[count4]);
				fscanf(fout[2],"%lf ",&out1);
          	    sum_of_uik = sum_of_uik + pow(out1,m);
                for(count4 = 1; count4<=num_of_img; count4++)
					v1[count2][count4] = v1[count2][count4]+pow(out1,m)*i[count4][0];
				if(count3 != total_pix)
				for(count4 = 0;count4<step;count4++)
					fscanf(fout[2],"%lf ",&out1);
	         	}

			printf("Sum of cluster %d (uik^m): %lf\n",count2,sum_of_uik);
			printf("mean of cluster %d :\n",count2);
			for(count3 = 1; count3<=num_of_img; count3++)
				{
				v1[count2][count3] = (double) v1[count2][count3]/sum_of_uik;
				printf("%lf\n", v1[count2][count3]);
				fseek(fin[count3],0,0);
				}
			fseek(fout[2],0,0);
		}//end of count2 (num_of_cluster)
      }//end of count1 (num_of_cycle)

      pix = 0;
	  max = 0;
      for(count1 = 1; count1 <= cluster_num; count1++)
        num_in_cover[count1] = 0;

      for(count1 = 0;count1 < total_pix;count1++)
		{
         for(count2 = 1;count2 <= cluster_num;count2++)
			{
            fscanf(fout[2],"%lf ",&class_type[count2]);
            if(class_type[count2]>max)
				{
				max = class_type[count2];
				pix = count2;
				}
			}
         fprintf(fout[1],"%c",pix);
         num_in_cover[pix]++;
         max = 0;
		}

      fprintf(fout[3],"m value = %lf\n",m);
      for(count1 = 1;count1 <= cluster_num;count1++)
		{
		fprintf(fout[3],"mean of cluster %d:\n",count1);
		for(count2 = 1; count2<=num_of_img; count2++)
			fprintf(fout[3],"%lf\n",v1[count1][count2]);
        fprintf(fout[3],"total number %d\n",num_in_cover[count1]);
        }

      for(count2 = 1; count2<=num_of_img; count2++)
         fclose(fin[count2]);
      fclose(fout[1]);
	  fclose(fout[2]);
      fclose(fout[3]);
      printf("\n===Running completed!\n");
}