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

#include "imgproc.h"
#define MaxScale 10

unsigned char **gDilation(unsigned char **image, short row, short col,
                   unsigned char *element[3], short num)
{
   /* f+g = max{ f(i+n,j+m)+g(n,m) } */
	short i,j,k;
	unsigned char **pimage;		/* pointer to the processed image */
    short max;
    short temp;
    short tmp_row, tmp_col;

    pimage = (unsigned char **)fspace_2d(row,col,sizeof(unsigned char));
    if(pimage == NULL) return NULL;

    for(i=0; i<row; i++) {
       for(j=0; j<col; j++) {
          max = 0;
          for(k=0; k<num; k++) {
             tmp_row = i+element[k][0];
             if(tmp_row<0) {
                j = col+1;
                break;
                }
             tmp_col = j+element[k][1];
             if(tmp_col<0) break;
             temp = image[tmp_row][tmp_col]+element[k][2];
             max = temp>max ? temp : max;
          }
          pimage[i][j] = max>255 ? 255 : (unsigned char)max;
       }
    }
    return pimage;
}

unsigned char **gErosion(unsigned char **image, short row, short col,
                   unsigned char *element[3], short num)
{
   /* f-g = min{ f(i+n,j+m)-g(n,m) } */
	short i,j,k;
	unsigned char **pimage;		/* pointer to the processed image */
        short min;
        short temp;
        short tmp_row, tmp_col;

        pimage = (unsigned char **)fspace_2d(row,col,sizeof(unsigned char));
        if(pimage == NULL) return NULL;

        for(i=0; i<row; i++) {
           for(j=0; j<col; j++) {
              min = 255;
              for(k=0; k<num; k++) {
                 tmp_row = i+element[k][0];
                 if(tmp_row<0) {
                    j = col+1;
                    break;
                    }
                 tmp_col = j+element[k][1];
                 if(tmp_col<0) break;

                 temp = image[tmp_row][tmp_col]+element[k][2];
                 min = temp<min ? temp : min;
              }
              pimage[i][j] = min<0 ? 0 : (unsigned char)min;
           }
        }
        return pimage;
}


void linefit(double *x, double *y, short n, double *k, double *b)
{
   short i;
   double a11,a12,a13,a21,a22,a23;
   double tmp;

   a11 = 0;
   a12 = 0;
   a13 = 0;
   a21 = 0;
   a22 = 0;
   a23 = 0;
   for(i=0; i<n; i++) {
      a11 += x[i]*x[i];
	  a12 += x[i];
	  a13 += -x[i]*y[i];
	  a21 += x[i];
	  a22 += 1;
	  a23 += -y[i];
   }
   
   tmp = a21*a12-a22*a11;
   if(fabs(tmp)>0.01) {
      *b = (a21*a13-a11*a23)/tmp;
   }

   *k = (a13-a12*(*b))/a11;
}

double **calfractal(unsigned char **image, short row, short col, 
           short minscale, short maxscale, unsigned char *element[3], short num)
{
	short          i,j,k,n,m;
	double         **area[MaxScale];          //pointer to saving area at all scale
	unsigned char  **temp;   //temporatory pointer
	unsigned char  **upper;  //
	unsigned char  **lower;  //
	short          winsize;  //the half size of window for calculating area; 
	short          current_scale;
	//double         
	short          tmp_row,tmp_col;
	double         *temp_x,*temp_y;
	double         temp_b,temp_k;
	
	
	double **result;       //save the fractal result 
	
	
	result = (double **)fspace_2d(row,col,sizeof(double));
	if(result == NULL) return NULL;
	
	//processing dilation at scale (minscale)
	upper = gDilation(image,row,col,element,num);
	if(upper == NULL) return NULL;
	for(i=1; i<minscale; i++){
	   temp = gDilation(upper,row,col,element,num);
	   if(temp == NULL) return NULL;
	   memcpy(upper,temp,row*col*sizeof(unsigned char));
	   ffree_2d(temp,row);    
	}
	
	//processing erosion at scale (minscale)
	lower = gErosion(image,row,col,element,num);
	if(lower == NULL) return NULL;
	for(i=1; i<minscale; i++){
	   temp = gErosion(image,row,col,element,num);
	   if(temp == NULL) return NULL;
	   memcpy(lower,temp,row*col*sizeof(unsigned char));
	   ffree_2d(temp,row);
	}
	
	//calculating the area
	// A(i,j,scale) = (
	current_scale = minscale;
    winsize = maxscale;
	for(k=0; k<maxscale+minscale+1; k++){
	   area[k] = (double **)fspace_2d(row,col,sizeof(double));
	   if(area[k] == NULL) return NULL;
	   
	   for(i=0; i<row; i++) {
	      for(j=0; j<col; j++) {
	         area[k][i][j] = 0.0;
	         for(n=-winsize; n<winsize; n++){
	            for(m=-winsize; m<winsize; m++) {
	               tmp_row = i+n;
	               if(tmp_row<0) {
	               	   n = winsize +1;
	               	   j = col + 1;
	               	   break;
	               	}
	               
	               tmp_col = j+m;
	               if(tmp_col<0){
	               	   n = winsize + 1;
	               	   break;
	               	}
				   area[k][i][j] += upper[tmp_row][tmp_col] - lower[tmp_row][tmp_col];
	               } // end of m
	            } //end of n
			 area[k][i][j] /= 2*current_scale;
	         } //end of j
	      } //end of i
	   current_scale++;
	   } //end of k

	//calculating the fractal at any pixels
	temp_x = (double *)calloc(maxscale-minscale+1,sizeof(double));
	temp_y = (double *)calloc(maxscale-minscale+1,sizeof(double));
	for(i=0; i<row; i++) {
		for(j=0; j<col; j++) {
			for(k=0; k<maxscale-minscale+1; k++) {
			   temp_x[k] = minscale+k;
			   temp_y[k] = area[k][i][j];
			}
			linefit(temp_x,temp_y,maxscale-minscale+1,&temp_k,&temp_b);
			result[i][j] = temp_k;
		}
	}
	return result;
	}


int main(int argc,char *argv[])
{
   unsigned char **image;
   char          file[256];
   double        **result;
//   FILE          *fp;

   printf("Enter the input file:\n");
   scanf("%s",file);
   
   image = InputImageWithName(&row,&col,file);
   if(image == NULL) {
      printf("Error in reading image file!\n");
	  return -1;
   }

   result == calfractal(image,row,col,minscale,maxscale,element,num);


}