/*
* THIS ROUTINE IS PART OF THE CLEMENTINE PDS FILE READER PROGRAM.
* IT WAS WRITTEN BY ACT CORP. IN DIRECT SUPPORT TO THE
* CLEMENTINE (DSPSE) PROGRAM.
*
* IF YOU FIND A PROBLEM OR MAKE ANY CHANGES TO THIS CODE PLEASE CONTACT
* Dr. Erick Malaret at ACT Corp.
* tel: (703) 742-0294
* (703) 683-7431
* email: nrlvax.nrl.navy.mil
*
*
*/
/*
Program to calculate the Discrete Cosine Transform of a 8x8 block of data
*/
#include
#ifdef __BORLANDC__
#include
#else
#include
#endif
#include
#include
#include "jpeg_c.h"
#define TRUE 1
#define FALSE 0
extern int OPTIMIZE;
extern int SLOWDECOMP;
long *DCTHist[64];
float *Rn[64];
float Q[64];
float C[64], Ct[64];
float qtable[64];
float U[64];
short ULS[64];
int zzseq[] = {0,1,8,16,9,2,3,10,17,24,32,25,18,11,4,5,12,19,26,33,40,
48,41,34,27,20,13,6,7,14,21,28,35,42,49,56,57,50,43,36,29,
22,15,23,30,37,44,51,58,59,52,45,38,31,39,46,53,60,61,54,
47,55,62,63};
void core(void);
void getDCTHist(BitStream *bs, long rows, long cols);
void getRn(void);
void initC(float *m, short N);
void transp(float *out, float *in, short N);
void matmult1(float *out, float *in1, float *in2, short N);
/****************************************************************************/
void decomp(BitStream *bs, CHARH *Image, long rows, long cols)
{
short i,j,k,indexi,indexip8,indexj,indexjp8;
float temp;
float V1[64];
long rowsleft, colsleft, icols;
short Done, Diff, Pred;
/*********************** Image Decompression *********************/
Done = FALSE;
rowsleft = rows;
colsleft = cols;
indexj = 0;
indexi = 0;
Pred = 0;
while ( !Done ) {
/* Decode block data */
decode(ULS,bs);
Diff = ULS[0] + Pred;
Pred = Diff;
ULS[0] = Diff;
/* Calculate Inverse Discrete Cosine Transform */
/*dequantize(U,ULS1,pt_qtable); 64 multiplications */
for (i=0; i<64; i++) {
if ( OPTIMIZE ) {
temp = Rn[i][ULS[i]];
if ( SLOWDECOMP )
U[zzseq[i]] = temp * Q[i];
else
U[zzseq[i]] = temp * qtable[i];
} else {
if ( SLOWDECOMP )
U[zzseq[i]] = ULS[i] * Q[i];
else
U[zzseq[i]] = ULS[i] * qtable[i];
}
}
if ( SLOWDECOMP ) {
matmult1(V1,Ct,U,8);
matmult1(U,V1,C,8);
} else {
core();
}
/*ilevelshift(U);*/
for (i=0; i<64; i++) {
U[i] += 128.0;
U[i] = floor( U[i] + 0.5 );
if ( U[i] > 255.0 ) U[i] = 255.0;
else if ( U[i] < 0.0 ) U[i] = 0.0;
else;
}
indexip8 = indexi + 8;
indexjp8 = indexj + 8;
if ( (rowsleft > 8) && (colsleft > 8) ) {
for (i=indexi, k=0; i < indexip8; i++)
for (j=indexj, icols=i*cols; j < indexjp8; j++, k++)
Image[icols+j] = (char)U[k];
indexj += 8;
colsleft -= 8;
}
else if ( (rowsleft > 8) && (colsleft <= 8) ) {
for (i=indexi, k=0; i < indexip8; i++) {
for (j=indexj, icols=i*cols; j < cols; j++, k++)
Image[icols+j] = (char)U[k];
k += (8 - colsleft);
}
indexj = 0;
indexi += 8;
rowsleft -= 8;
colsleft = cols;
}
else if ( (rowsleft <= 8) && (colsleft > 8) ) {
for (i=indexi, k=0; i < rows; i++)
for (j=indexj, icols=i*cols; j < indexjp8; j++, k++)
Image[icols+j] = (char)U[k];
indexj += 8;
colsleft -= 8;
}
else {
for (i=indexi, k=0; i < rows; i++) {
for (j=indexj, icols=i*cols; j < cols; j++, k++)
Image[icols+j] = (char)U[k];
k += (8 - colsleft);
}
Done = TRUE;
}
}
}
void getDCTHist(BitStream *bs, long rows, long cols)
{
short Pred, i;
long nblocks;
nblocks = (rows * cols) / 64;
Pred = 0;
while ( nblocks ) {
/* Decode block data */
decode(ULS,bs);
ULS[0] += Pred;
Pred = ULS[0];
for (i=0; i<64; i++) DCTHist[i][ULS[i]]++;
nblocks--;
}
}
void getRn()
{
short i, j;
float lb, ub, a, b, num, den, phi;
float m, pa, pap, pb, pbp, pm, pmp;
for (i=0; i<64; i++) {
Rn[i][-256] = (float)-256;
Rn[i][256] = (float)256;
/*
for (j=-255; j<256; j++) {
if ( DCTHist[i][j] > 0 ) {
m = ((float)j)*Q[i];
pm = DCTHist[i][j];
pmp = m*pm;
a = ((float)j - 0.5)*Q[i];
phi = 0.5;
pa = phi*DCTHist[i][j] + (1.0-phi)*DCTHist[i][j-1];
pap = ((float)j*pa - (1.0-phi)*DCTHist[i][j-1])*Q[i];
b = ((float)j + 0.5)*Q[i];
phi = 0.5;
pb = phi*DCTHist[i][j+1] + (1.0-phi)*DCTHist[i][j];
pbp = ((float)j*pb + phi*DCTHist[i][j+1])*Q[i];
num = (m-a)*(pap+pmp)+(b-m)*(pbp+pmp);
den = (m-a)*(pa+pm)+(b-m)*(pb+pm);
Rn[i][j] = (num/den)/Q[i];
} else {
Rn[i][j] = (float)j;
}
}
*/
for (j=-255; j<256; j++) {
if ( DCTHist[i][j] > 0 ) {
lb = ((float)j - 0.5)*Q[i];
ub = ((float)j + 0.5)*Q[i];
m = ((float)j)*Q[i];
pm = DCTHist[i][j];
pmp = m*pm;
a = ceil(lb);
phi = a/Q[i] - (float)(j-1);
pa = phi*DCTHist[i][j] + (1.0-phi)*DCTHist[i][j-1];
pap = ((float)j*pa - (1.0-phi)*DCTHist[i][j-1])*Q[i];
b = ceil(ub);
phi = b/Q[i] - (float)j;
pb = phi*DCTHist[i][j+1] + (1.0-phi)*DCTHist[i][j];
pbp = ((float)j*pb + phi*DCTHist[i][j+1])*Q[i];
num = (m-a)*(pap+pmp)+(b-m)*(pbp+pmp);
den = (m-a)*(pa+pm)+(b-m)*(pb+pm);
Rn[i][j] = (num/den)/Q[i];
} else {
Rn[i][j] = (float)j;
}
}
}
}
#define PI 3.141592653589793
void initC(float *m, short N)
{
short i, j, k;
for (i=0, k=0; i
http://clementine.cnes.fr/software/decompression_new/decomp1.c
(possibly inaccurate URL)
07/1998
