dct.c

DCT DCT by using c

/* ========================================================================= */
/* =========           Image transformation using DCT           ============ */
/* =========              for 512 * 512 images                  ============ */
/* =========                                                    ============ */
/* =========            Input: 512*512 image                    ============ */
/* =========            Output: transformed image               ============ */
/* =========                                                    ============ */
/* =========            Last updated: 10 May, 2006              ============ */
/* ========================================================================= */


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

#define dim         8
#define imagesize   512
#define vectorsize  dim*dim
#define vectornum   (imagesize*imagesize)/(vectorsize)
#define pi          3.14159


double        C[dim];               /* store the coefficient of DCT function */
double        dest[dim][dim];       /* store the transformed coefficients in a block*/
double        source[dim][dim];     /* store the spatial pixels in a block */

unsigned char image[imagesize][imagesize];           /* store the original image */
unsigned char Cimage[imagesize][imagesize];          /* store the processed image */
unsigned char vectors[vectornum][vectorsize];        /* store the original vectors */
float         Tvectors[vectornum][vectorsize];       /* store the transformed values */
unsigned char Rvectors[vectornum][vectorsize];       /* store the recovered values */
unsigned long squar[256];                            /* store the squared values */

FILE *imagefile,*result;


main(argc,argv)
int argc;
char **argv;
{
  int i,j;
  int bno;

  if (argc != 3)
     {
       printf("\nError : The number of parameters isn't correct ! \n");
       printf("\nFormat: Command input_grayscale_image output_image\n");
       exit(1);
     }
     
  if ((imagefile=fopen(argv[1],"rb"))==NULL)
     {
       printf("Unable to open input image %s.\n",argv[1]);
       exit(1);
     }
     
  if ((result=fopen(argv[2],"wb"))==NULL)
     {
       printf("Unable to open output image %s.\n",argv[2]);
       exit(1);
     }

  /* get the original image */
  for(i=0;i<imagesize;i++)
     {
       for(j=0;j<imagesize;j++)
          {
            image[i][j]=fgetc(imagefile);
          }
     }   

  /* transfor image to vector form */
  for(i=0;i<imagesize;i++)
     {
       for(j=0;j<imagesize;j++)
          {
            vectors[j/8+(i/8)*(64)][(j%8)+(i%8)*8]=image[i][j];
          }
     }
     
  for(i=0;i<256;i++)
     {
       squar[i]=(unsigned long)(i*i);
     }
       
  for (i=0; i<dim;i++)
      {
        if (i==0)
           { 
             C[i]= 1.0/sqrt(2.0); 
           } 
        else
           { 
             C[i]= 1;
           }
      }
         
           
  FDCT_transform();

  IDCT_transform();
 
  /* ========= transform vectors back to image =========== */
  for(i=0; i<vectornum; i++)
     {
       for(j=0; j< vectorsize; j++)
          {
            Cimage[j/8+(i/64)*8][j%8+(i%64)*8]=(unsigned char)Rvectors[i][j];
          }
     }

  fwrite(Cimage,sizeof(unsigned char),imagesize*imagesize, result);

  Find_PSNR();

  printf("\nPlease input the number of image blocks to be shown\n");
  scanf("%d",&bno);
  
  ShowBlock(bno);
   
  fclose(imagefile);
  fclose(result);

 }

 



FDCT_transform()
{
  int    i,j,k,m,n;
  float  temp1,temp2,temp3;

  printf("\nEnter forward DCT transformation\n\t ");
  for(i=0; i<vectornum; i++)
     {
       /* transform the input vector back to block of (dim,dim) */
       for(j=0; j<vectorsize; j++)
          {
            source[j/dim][j%dim]=vectors[i][j];
          } 
      
       for(m=0; m<dim; m++)
          { 
            for(n=0; n<dim; n++)
               {
                 temp3 = 0;
                 dest[m][n] = 0;
                 for(j=0; j<dim; j++)
                    { 
                      for(k=0; k<dim; k++)
                         {
                           temp1 = cos( ((2*j+1)*m*pi) /(2*dim) );
                           temp2 = cos( ((2*k+1)*n*pi) /(2*dim) );
                           temp3 += source[j][k] * temp1 * temp2 ;
                         }
                    }
                 dest[m][n] = (2.0/(double)dim) * C[m] * C[n] * temp3;
               }
          }

       /* store the result in Tvectors matrix */
       for(j=0; j<vectorsize; j++)
          { 
            Tvectors[i][j] = (int) dest[j/dim][j%dim];
          }
     } 

}






IDCT_transform()
{
  int    i,j,k,m,n;
  float  temp1,temp2,temp3;
 
  printf("\nEnter inverse DCT transformation\n\t");
  for(i=0; i<vectornum; i++)
     {
       /* transform the input vector back to block of (dim,dim) */
       for(j=0; j<vectorsize; j++)
          {
            source[j/dim][j%dim] = Tvectors[i][j];
          } 
      
       for(j=0; j<dim; j++)
          { 
            for(k=0; k<dim; k++)
               {
            temp3 = 0;
            dest[j][k] = 0;
            for(m=0; m<dim; m++)
               { 
                 for(n=0; n<dim; n++)
                    {
                      temp1 = cos( ((2*j+1)*m*pi) /(2*dim) );
                      temp2 = cos( ((2*k+1)*n*pi) /(2*dim) );
                      temp3 += C[m]*C[n]*source[m][n]*temp1*temp2;
                    }
               }
            dest[j][k] = (2.0/(double)dim) * temp3;
          }
    }

    /* store the result in QTTvectors matrix */
  for(j=0; j<vectorsize; j++)
     { 
       Rvectors[i][j] = (int) dest[j/dim][j%dim];
     }
    }
}






ShowBlock(bno)
int bno;
{
 int        i,j;
 
 printf("\nThe pixels of block %3d are listed in the following",bno);
 for(i=0;i<vectorsize;i++)
    {
      if (i%8==0) printf("\n\t");
      printf("%4d ",vectors[bno][i]);
    }  
 

 printf("\n\nTransformed coefficients of block %3d are listed in the following",bno);
 for(i=0;i<vectorsize;i++)
    {
      if (i%8==0) printf("\n\t");
      printf("%6.1f",Tvectors[bno][i]);
    }  



 printf("\n\nRecovered pixels of block %3d are listed in the following",bno);
 for(i=0;i<vectorsize;i++)
    {
      if (i%8==0) printf("\n\t");
      printf("%4d ",Rvectors[bno][i]);
    }  

  
}






Find_PSNR()
{
 int        i,j,k;
 long int   dist;	 
 float      ave,psnr;

 dist=0L;
 for(i=0;i<imagesize;i++)
    {
      for(j=0;j<imagesize;j++)
         {
           dist+=(long)squar[abs(image[i][j]-Cimage[i][j])];
         }
    }
   
 ave=(float)dist/(imagesize*imagesize);
   
 psnr = 10*log10((double)255*255/ave);
 
 printf("\ntotal_dist=%10ld avg_dist=%8.3f  PSNR=%8.3f\n",dist,ave,psnr);

}