dprdct.c

JPEG Image compression using IJG standards followed

/*
 ********************************************************************
 *
 * File Name:		dprdct.c
 *
 * Designed by:		
 *					Jie Liang
 *
 * Created Date:	June 30, 2000
 *
 * Contents:
 *					Implement the 2D  8x8 inverse lossless binDCT.
 *
 * Usage:	dprdct dctfile recfile hsize vsize
 *		"Arguments:\n"
 *		"	dctfile: 	input DCT-transformed file,\n"
 *		"	recfile:	reconstructed image file,\n"
 *		"	hsize:		Horizontal size of the image,\n"
 *		"	vsize:		Vertical size of the image.\n"
 *		"Examples:\n"
 *		"	./dprdct lena.dct lena.rec 512 512\n"
 *			
 *	
 * Limitations:
 *				
 *			
 * Copyright (C) 2000 Department of Electrical and Computer Engineering 
 * The Johns Hopkins University. All right reserved. 
 *
 *********************************************************************
 */

/*
 *********************************************************************
 *
 *   Modification History:
 *
 *   Date	Programmer	Description
 *-------------------------------------------------------------------
 *
 * $Log$
 *
 *
 *
 *********************************************************************
 */

static const char rcsid[] = 
	"$Id$";

/*
 ********************************************************************
 *
 *	Include Files
 *
 ********************************************************************
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>


/*
 *******************************************************************
 *
 *	Variables
 *
 *******************************************************************
 */

/*
 *******************************************************************
 *
 * 	Local Constants
 *
 *******************************************************************
 */

/*
 *******************************************************************
 *
 *	Local Macros
 *
 *******************************************************************
 */

/*
 *******************************************************************
 *
 *	Local Enumerated Types
 *
 *******************************************************************
 */

/*
 *******************************************************************
 *
 *	Local Data Types
 *
 *******************************************************************
 */
#define DCTSIZE   8
typedef signed short int DCTELEM;

/*
 *******************************************************************
 *
 *	Local Function Prototypes
 *
 *******************************************************************
 */

void jpeg_idct_bin_a1_pr (DCTELEM *coef_block, unsigned char *output_buf) ;

/*
 ***********************************************************************
 *
 * Function Name:		main
 *
 * Descriptions:
 *				open input image, perform DCT transformation
 *		
 *
 * Input Parameters:
 *				argc: argument count;
 *				argv: array of input arguments.
 * Output Parameters:		0:    successful,
 *				-1:   error.
 * Limitations:
 *
 ************************************************************************
 */
int main(int argc, char *argv[])
{
  FILE 				*dctfile, *recfile;
  unsigned char		row, col, rblks, cblks, rbidx, cbidx;	/* horizontal, vertical blocks and indexes */
  int				i, j, hsize, vsize, offset, qtzer, qtzrate;	/* source symbols */ 

  unsigned char		*imgptr_chr, outbuf[64], *tmpoutptr;
  DCTELEM			*imgptr, *hdptr, *tmptr, dct_1d_buf[64];
  clock_t			timer0, timer1;


  printf("\nInverse lossless binDCT:\n");

  if (argc != 5) {
	printf(
		"\nUsage:\n dprdct dctfile dctfile hsize vsize\n"
		"Arguments:\n"
		"	dctfile: 	DCT-transformed file,\n"
 		"	recfile:	reconstructed image file to be generated,\n"
		"	hsize:		Horizontal size of the image,\n"
		"	vsize:		Vertical size of the image.\n"
		"Examples:\n"
		"	./dprdct lena.dct lena_rec.raw 512 512\n"
	);	
	return(-1);
  }
 

  /* Open Input files */
  dctfile = fopen(argv[1], "r"); 
  if (dctfile == NULL)  {
	fprintf(stderr, "Input file not found.\n");
	return(-1);
  }
  
  /*create output file */
  recfile = fopen(argv[2], "w");
  if( recfile == NULL ) {
  	fprintf(stderr, "Error in generating output file.\n");
	return(-1);
  }

  /* read image size */
  hsize = atoi(argv[3]);
  vsize = atoi(argv[4]);

  if ( (hsize & 7) || (vsize & 7) ) {
	fprintf(stderr, "Image size paramemters should be multiples of 8.\n");
	return(-1);
  }
	
  /* block indexes in two directions */ 
  cblks = hsize / 8 ;
  rblks = vsize / 8 ;

  /* allocate memory for the img array */
  if( !(imgptr_chr = (unsigned char *) malloc(hsize * vsize)) ||
	  !(imgptr = (DCTELEM *) malloc(hsize * vsize * sizeof(DCTELEM) )) ) {
	printf("Memory allocation error.\n");
	return(-1);
  }  

  /* read image data */
  fread(imgptr, sizeof(DCTELEM), hsize * vsize, dctfile);
  fclose(dctfile);

  for (i=0; i<64;i++)
	printf("%5d", *(imgptr+i));

  /*
   *************************************************
   *
   * Inverse DCT
   *
   *************************************************
   */

  timer0 = clock();

  for(rbidx = 0; rbidx < rblks; rbidx++) {
	for (cbidx = 0; cbidx < cblks; cbidx++) {
		
	  // pointer of input block
		hdptr = imgptr +  rbidx * hsize * 8 + cbidx * 8 ;	

		//copy block to a 8x8 array
		for (i = 0; i < 8; i++) {
		  tmptr = hdptr + i * hsize;
		  for (j = 0; j < 8; j++) {
			dct_1d_buf[8*i + j] = *(tmptr + j);
		  }
		}

  for (i=0; i<64;i++)
	printf("%5d", dct_1d_buf[i]);
  printf("\n\n");

		jpeg_idct_bin_a1_pr((DCTELEM *)&dct_1d_buf, (unsigned char *)&outbuf);

		//copy output in a 8x8 array to output image block.
		for (i = 0; i < 8; i++) {
		  tmpoutptr = (imgptr_chr  +  rbidx * hsize * 8 + cbidx * 8) + i * hsize;
		  for (j = 0; j < 8; j++) {
			*(tmpoutptr + j) = outbuf[8*i + j];
			printf("%5d",*(tmpoutptr + j) );
		  }
		}

	} 
  }	

  timer1 = clock();
  printf("Time used by Inverse lossless binDCT is: %10.6f sec.\n", ((float) (timer1 - timer0)) / (float )CLOCKS_PER_SEC);

  fwrite(imgptr_chr, sizeof(unsigned char), hsize * vsize, recfile);
  fclose(recfile);

  printf("\nReconstructed file generated.\n");

   return(0);
}

/*
 ***********************************************************************
 *
 * Function Name:		inv_fastdct
 *
 * Descriptions:
 *				perform 8-point DCT for one row data of unsigned char type.
 *				A(i, :) * C_transpose			
 *		
 * Input Parameters:
 *				data:  pointer to the given data.
 * Return Parameters:	none.
 *
 * Notice:
 *				Because butfly * butfly = 2I, the fast DCT will scale
 *				the imgae by a factor of 2.
 *
 ************************************************************************
 */

void jpeg_idct_bin_a1_pr (DCTELEM *coef_block, unsigned char *output_buf) 
{
  DCTELEM tmp0, tmp1, tmp2, tmp3,tmp4,tmp5,tmp6,tmp7;
  DCTELEM tmp10, tmp11, tmp12, tmp13;

  DCTELEM *inptr, *wsptr;
  int ctr;
  unsigned char *outptr;	

  /* Pass 1: process columns from input, store into work array. */
  /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  /* furthermore, we scale the results by 2**PASS1_BITS. */

  inptr = coef_block;
  wsptr = coef_block;

  for (ctr = DCTSIZE; ctr > 0; ctr--) {

	if (inptr[DCTSIZE*1] == 0 && inptr[DCTSIZE*2] == 0 &&
	inptr[DCTSIZE*3] == 0 && inptr[DCTSIZE*4] == 0 &&
	inptr[DCTSIZE*5] == 0 && inptr[DCTSIZE*6] == 0 &&
	inptr[DCTSIZE*7] == 0) {
      int dcval = inptr[DCTSIZE*0] >> 1;
      
      wsptr[DCTSIZE*0] = dcval;
      wsptr[DCTSIZE*1] = dcval;
      wsptr[DCTSIZE*2] = dcval;
      wsptr[DCTSIZE*3] = dcval;
      wsptr[DCTSIZE*4] = dcval;
      wsptr[DCTSIZE*5] = dcval;
      wsptr[DCTSIZE*6] = dcval;
      wsptr[DCTSIZE*7] = dcval;
      
      inptr++;			
      continue;
    }

    tmp0 = inptr[DCTSIZE*0];
    tmp1 = inptr[DCTSIZE*4];
	tmp2 = inptr[DCTSIZE*6];
    tmp3 = inptr[DCTSIZE*2];
    tmp4 = inptr[DCTSIZE*7];
    tmp5 = inptr[DCTSIZE*5];
    tmp6 = inptr[DCTSIZE*3];
    tmp7 = inptr[DCTSIZE*1];

	/*
    fprintf(stderr, "%10d", tmp0);
    fprintf(stderr, "%10d", tmp7);
    fprintf(stderr, "%10d", tmp3);
    fprintf(stderr, "%10d", tmp6);
    fprintf(stderr, "%10d", tmp1);
    fprintf(stderr, "%10d", tmp5);
    fprintf(stderr, "%10d", tmp2);
    fprintf(stderr, "%10d", tmp4);
	fprintf(stderr, "\n");	
	*/

	/* X[0] and X[4] */
	tmp11 = ((tmp0 + 1) >> 1) - tmp1;
	tmp10 = tmp0 - tmp11;
	
	/* X[6] and X[2] */
	tmp13 = tmp3 + (((tmp2 << 1) + tmp2 + 4) >> 3);
	tmp12 = ((tmp13 + 1) >> 1) - ((tmp13 + 8) >> 4) - tmp2;

	tmp0 = tmp10 + tmp13;
	tmp3 = tmp10 - tmp13;
	tmp1 = tmp11 + tmp12;
	tmp2 = tmp11 - tmp12;

	/* X[7] and X[1]: */
	/* 7pi/16 = 3/16d 3/16u */
    tmp13 = tmp7 + ( ((tmp4 << 1) + tmp4 + 8) >> 4 );
	tmp10 = ( ((tmp13 << 1) + tmp13 + 8) >> 4 ) - tmp4;

	/* X[5] and X[3] */
	/* 3pi/16 = 1/2d -7/8u */
	/*    tmp12 = tmp6 + ((tmp5 + 1) >> 1);
		  tmp11 = tmp5 - tmp12 + ((tmp12 + 4) >> 3);*/

	/* new 7/16 and -5/8*/
	tmp12 = tmp6 + ((tmp5 + 1) >> 1) - ((tmp5 + 8) >> 4);
    tmp11 = tmp5 - ((tmp12 + 1) >> 1) - ((tmp12 + 4) >> 3);

	/* Butterfly */
	tmp4 = tmp10 + tmp11;
	tmp5 = tmp10 - tmp11;
	tmp6 = tmp13 - tmp12;
	tmp7 = tmp13 + tmp12;

	/* pi/4 = -3/8u -11/16d 7/16u */
	tmp5 = (((tmp6 << 1) + tmp6 + 4) >> 3) - tmp5;
	tmp6 = tmp6 - tmp5 + ((tmp5 + 2) >> 2) + ((tmp5 + 8) >> 4);
	tmp5 = tmp5 + ((tmp6 + 1) >> 1) - ((tmp6 + 8) >> 4);

	/* last stage: butterfly */
	wsptr[DCTSIZE*0] = (tmp0 + tmp7);
    wsptr[DCTSIZE*7] = (tmp0 - tmp7);
    wsptr[DCTSIZE*1] = (tmp1 + tmp6);
    wsptr[DCTSIZE*6] = (tmp1 - tmp6);
    wsptr[DCTSIZE*2] = (tmp2 + tmp5);
    wsptr[DCTSIZE*5] = (tmp2 - tmp5);
    wsptr[DCTSIZE*3] = (tmp3 + tmp4);
    wsptr[DCTSIZE*4] = (tmp3 - tmp4);
    
    inptr++;			/* advance pointers to next column */

  }
  
  /* Pass 2: process rows from work array, store into output array. */
  /* Note that we must descale the results by a factor of 8 == 2**3, */
  /* and also undo the PASS1_BITS scaling. */

  //fprintf(stderr, "\nAfter inverse DCT:\n");

  wsptr = coef_block;
  outptr = output_buf;

  for (ctr = 0; ctr < DCTSIZE; ctr++) {
	outptr = outptr + ctr*8;

    /* Rows of zeroes can be exploited in the same way as we did with columns.
     * However, the column calculation has created many nonzero AC terms, so
     * the simplification applies less often (typically 5% to 10% of the time).
     * On machines with very fast multiplication, it's possible that the
     * test takes more time than it's worth.  In that case this section
     * may be commented out.
     */
    
#ifndef NO_ZERO_ROW_TEST
    if (wsptr[1] == 0 && wsptr[2] == 0 && wsptr[3] == 0 && wsptr[4] == 0 &&
	wsptr[5] == 0 && wsptr[6] == 0 && wsptr[7] == 0) {

	  /* if all AC are 0, the IDCT will all equal to 1/2 DC, so downscale by 2,
		 After that, apply the downscale of 16 caused by butterflies, so total downscale = 32.*/
	  unsigned char dcval = wsptr[0] >> 5;
      outptr[0] = dcval;
      outptr[1] = dcval;
      outptr[2] = dcval;
      outptr[3] = dcval;
      outptr[4] = dcval;
      outptr[5] = dcval;
      outptr[6] = dcval;
      outptr[7] = dcval;

      wsptr += DCTSIZE;		
      continue;
    }
#endif
    
    /* Even part: reverse the even part of the forward DCT. */
    /* The rotator is sqrt(2)*c(-6). */

    /* Even part */
/**********************/
/* not necessary ??? */
/************************/

/********
    tmp0 = (INT32) wsptr[0];
    tmp1 = (INT32) wsptr[4];
    tmp2 = (INT32) wsptr[6];
    tmp3 = (INT32) wsptr[2];
    tmp4 = (INT32) wsptr[7];
    tmp5 = (INT32) wsptr[5];
    tmp6 = (INT32) wsptr[3];
    tmp7 = (INT32) wsptr[1];
*********/

	/* X[0] and X[4] */
	tmp11 = ((wsptr[0] + 1) >> 1) - wsptr[4];
	tmp10 = wsptr[0] - tmp11;
	
	/* X[6] and X[2] */
	tmp13 = wsptr[2] + (((wsptr[6] << 1) + wsptr[6] + 4) >> 3);
	tmp12 = ((tmp13 + 1) >> 1) - ((tmp13 + 8) >> 4) - wsptr[6];

	tmp0 = tmp10 + tmp13;
	tmp3 = tmp10 - tmp13;
	tmp1 = tmp11 + tmp12;
	tmp2 = tmp11 - tmp12;

	/* 7pi/16 = -3/16d 3/16u */
    tmp13 = wsptr[1] +( ((wsptr[7] << 1) + wsptr[7] + 8) >> 4 );
	tmp10 = ( ((tmp13 << 1) + tmp13 + 8) >> 4 ) - wsptr[7];

	/* 3pi/16 = 1/2d -7/8u */
	/*    tmp12 = wsptr[3] + ((wsptr[5] + 1) >> 1);
		  tmp11 = wsptr[5] - tmp12 + ((tmp12 + 4) >> 3); */

	/* new 7/16 and -5/8*/
	tmp12 = wsptr[3] + ((wsptr[5] + 1) >> 1) - ((wsptr[5] + 8) >> 4);
    tmp11 = wsptr[5] - ((tmp12 + 1) >> 1) - ((tmp12 + 4) >> 3);

	tmp4 = tmp10 + tmp11;
	tmp5 = tmp10 - tmp11;
	tmp6 = tmp13 - tmp12;
	tmp7 = tmp13 + tmp12;

	/* pi/4 = -3/8u -11/16d 7/16u */
	tmp5 = (((tmp6 << 1) + tmp6 + 4) >> 3) - tmp5;
	tmp6 = tmp6 - tmp5 + ((tmp5 + 2) >> 2) + ((tmp5 + 8) >> 4);
	tmp5 = tmp5 + ((tmp6 + 1) >> 1) - ((tmp6 + 8) >> 4);

	/* last stage: butterfly */

    /* Final output stage: scale down by a factor of 8 and range-limit */
    outptr[0] = (tmp0 + tmp7) >> 4;
    outptr[7] = (tmp0 - tmp7) >> 4;

    outptr[1]=(tmp1 + tmp6) >> 4;
    outptr[6]=(tmp1 - tmp6) >> 4;

    outptr[2]=(tmp2 + tmp5) >> 4;
    outptr[5]=(tmp2 - tmp5) >> 4;

    outptr[3]=(tmp3 + tmp4) >> 4;
    outptr[4]=(tmp3 - tmp4) >> 4;
     
    wsptr += DCTSIZE;		/* advance pointer to next row */

/*******************/
/* Jie: test code */
	/*	for (tmp0 = 0; tmp0 < 8; tmp0 ++) {
	  fprintf(stderr, "%10d", outptr[tmp0]);
	}
	fprintf(stderr, "\n");
	*/
  }
}




/*
 *******************************************************
 *
 * End of $Source$
 *
 *******************************************************
 */