求离散余弦变换算法代码，原理（图象处理用）

[zhonghui]

DCT和IDCT算法代码。

[Kingore]

<VISUAL c++ 数字图象处理》一书有讲。

[Kingore]

/*************************************************************************
*
* º¯ÊýÃû³Æ£º
* DCT()
*
* ²ÎÊý:
* double * f - Ö¸ÏòʱÓòÖµµÄÖ¸Õë
* double * F - Ö¸ÏòƵÓòÖµµÄÖ¸Õë
* r £­2µÄÃÝÊý
*
* ·µ»ØÖµ:
* ÎÞ¡£
*
* ˵Ã÷:
* ¸Ãº¯ÊýÓÃÀ´ÊµÏÖ¿ìËÙÀëÉ¢ÓàÏұ任¡£¸Ãº¯ÊýÀûÓÃ2NµãµÄ¿ìËÙ¸¶Á¢Ò¶±ä»»
* À´ÊµÏÖÀëÉ¢ÓàÏұ任¡£
*
************************************************************************/
VOID WINAPI DCT(double *f, double *F, int r)
{
// ÀëÉ¢ÓàÏұ任µãÊý
LONG count;

// Ñ­»·±äÁ¿
int i;

// Öмä±äÁ¿
double dTemp;

complex<double> *X;

// ¼ÆËãÀëÉ¢ÓàÏұ任µãÊý
count = 1<<r;

// ·ÖÅäÄÚ´æ
X = new complex<double>[count*2];

// ¸³³õֵΪ0
memset(X, 0, sizeof(complex<double>) * count * 2);

// ½«Ê±ÓòµãдÈëÊý×éX
for(i=0;i<count;i++)
{
X[i] = complex<double> (f[i], 0);
}

// µ÷ÓÿìËÙ¸¶Á¢Ò¶±ä»»
FFT(X,X,r+1);

// µ÷ÕûϵÊý
dTemp = 1/sqrt(count);

// ÇóF[0]
F[0] = X[0].real() * dTemp;

dTemp *= sqrt(2);

// ÇóF[u]
for(i = 1; i < count; i++)
{
F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;
}

// ÊÍ·ÅÄÚ´æ
delete X;
}

/*************************************************************************
*
* º¯ÊýÃû³Æ£º
* IDCT()
*
* ²ÎÊý:
* double * F - Ö¸ÏòƵÓòÖµµÄÖ¸Õë
* double * f - Ö¸ÏòʱÓòÖµµÄÖ¸Õë
* r £­2µÄÃÝÊý
*
* ·µ»ØÖµ:
* ÎÞ¡£
*
* ˵Ã÷:
* ¸Ãº¯ÊýÓÃÀ´ÊµÏÖ¿ìËÙÀëÉ¢ÓàÏÒ·´±ä»»¡£¸Ãº¯ÊýÒ²ÀûÓÃ2NµãµÄ¿ìËÙ¸¶Á¢Ò¶±ä»»
* À´ÊµÏÖÀëÉ¢ÓàÏÒ·´±ä»»¡£
*
************************************************************************/
VOID WINAPI IDCT(double *F, double *f, int r)
{
// ÀëÉ¢ÓàÏÒ·´±ä»»µãÊý
LONG count;

// Ñ­»·±äÁ¿
int i;

// Öмä±äÁ¿
double dTemp, d0;

complex<double> *X;

// ¼ÆËãÀëÉ¢ÓàÏұ任µãÊý
count = 1<<r;

// ·ÖÅäÄÚ´æ
X = new complex<double>[count*2];

// ¸³³õֵΪ0
memset(X, 0, sizeof(complex<double>) * count * 2);

// ½«ÆµÓò±ä»»ºóµãдÈëÊý×éX
for(i=0;i<count;i++)
{
X[i] = complex<double> (F[i] * cos(i*PI/(count*2)), F[i] * sin(i*PI/(count*2)));
}

// µ÷ÓÿìËÙ¸¶Á¢Ò¶·´±ä»»
IFFT(X,X,r+1);

// µ÷ÕûϵÊý
dTemp = sqrt(2.0/count);
d0 = (sqrt(1.0/count) - dTemp) * F[0];

// ¼ÆËãf(x)
for(i = 0; i < count; i++)
{
f[i] = d0 + X[i].real()* dTemp * 2 * count;
}

// ÊÍ·ÅÄÚ´æ
delete X;
}

[Kingore]

/*************************************************************************
*
* º¯ÊýÃû³Æ£º
* DIBDct()
*
* ²ÎÊý:
* LPSTR lpDIBBits - Ö¸ÏòÔ´DIBͼÏñÖ¸Õë
* LONG lWidth - ԴͼÏñ¿í¶È£¨ÏóËØÊý£©
* LONG lHeight - ԴͼÏñ¸ß¶È£¨ÏóËØÊý£©
*
* ·µ»ØÖµ:
* BOOL - ³É¹¦·µ»ØTRUE£¬·ñÔò·µ»ØFALSE¡£
*
* ˵Ã÷:
* ¸Ãº¯ÊýÓÃÀ´¶ÔͼÏñ½øÐÐÀëÉ¢ÓàÏұ任¡£
*
************************************************************************/
BOOL WINAPI DIBDct(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{

// Ö¸ÏòԴͼÏñµÄÖ¸Õë
unsigned char* lpSrc;

// Ñ­»·±äÁ¿
LONG i;
LONG j;

// ½øÐи¶Á¢Ò¶±ä»»µÄ¿í¶ÈºÍ¸ß¶È£¨2µÄÕûÊý´Î·½£©
LONG w;
LONG h;

// Öмä±äÁ¿
double dTemp;

int wp;
int hp;

// ͼÏñÿÐеÄ×Ö½ÚÊý
LONG lLineBytes;

// ¼ÆËãͼÏñÿÐеÄ×Ö½ÚÊý
lLineBytes = WIDTHBYTES(lWidth * 8);

// ¸³³õÖµ
w = 1;
h = 1;
wp = 0;
hp = 0;

// ¼ÆËã½øÐÐÀëÉ¢ÓàÏұ任µÄ¿í¶ÈºÍ¸ß¶È£¨2µÄÕûÊý´Î·½£©
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

// ·ÖÅäÄÚ´æ
double *f = new double[w * h];
double *F = new double[w * h];

// ÐÐ
for(i = 0; i < h; i++)
{
// ÁÐ
for(j = 0; j < w; j++)
{
// Ö¸ÏòDIBµÚiÐУ¬µÚj¸öÏóËصÄÖ¸Õë
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// ¸øʱÓò¸³Öµ
f[j + i * w] = *(lpSrc);
}
}

for(i = 0; i < h; i++)
{
// ¶Ôy·½Ïò½øÐÐÀëÉ¢ÓàÏұ任
DCT(&f[w * i], &F[w * i], wp);
}

// ±£´æ¼ÆËã½á¹û
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
f[j * h + i] = F[j + w * i];
}
}

for(j = 0; j < w; j++)
{
// ¶Ôx·½Ïò½øÐÐÀëÉ¢ÓàÏұ任
DCT(&f[j * h], &F[j * h], hp);
}

// ÐÐ
for(i = 0; i < h; i++)
{
// ÁÐ
for(j = 0; j < w; j++)
{
// ¼ÆËãƵÆ×
dTemp = fabs(F[j*h+i]);

// ÅжÏÊÇ·ñ³¬¹ý255
if (dTemp > 255)
{
// ¶ÔÓÚ³¬¹ýµÄ£¬Ö±½ÓÉèÖÃΪ255
dTemp = 255;
}

// Ö¸ÏòDIBµÚyÐУ¬µÚx¸öÏóËصÄÖ¸Õë
// lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);

// ¸üÐÂԴͼÏñ
* (lpSrc) = (BYTE)(dTemp);
}
}

// ÊÍ·ÅÄÚ´æ
delete f;
delete F;

// ·µ»Ø
return TRUE;
}

[zhonghui]

这么长？：）

[hahu]

我没完成的Jpeg编码
里面有FDCT,IDCT,Z行程编码，逆
void CJpegImage::FDCT(short pSrcDct[][8], double pDstDct[][8])
{
TRACE("FDCT:\n");
int x, y, u, v;
for(x=0;x<8;x++)
{
for(y=0;y<8;y++)
{

m_pTemDct[x][y]=pSrcDct[x][y]-128;
// TRACE("%f ",m_pTemDct[x][y]);
}
// TRACE("\n");
}
double persqrt2=1.0/sqrt(2.0);
TRACE("persqrt2=%f\n",persqrt2);
double cu,cv;
double pi = 3.14159265;
for( u = 0; u < 8; u++ )
{
for( v = 0; v < 8; v++ )
{
pDstDct[u][v] = 0;
for( x = 0; x < 8; x++ )
{
for( y = 0; y < 8; y++ )
{
pDstDct[u][v] += m_pTemDct[x][y] * cos((2.0*x+1.0)*u*pi/16.0)
*cos((2.0*y+1.0)*v*pi/16.0);
}
}
cu=(u==0)?persqrt2:1;
cv=(v==0)?persqrt2:1;

pDstDct[u][v] = pDstDct[u][v]*cu*cv/4.0;
}
}
}


void CJpegImage::FQuantify(double pSrc[][8],short int pDst[8][8])
{
TRACE("FQuantify:\n");
int u, v;
for( u = 0; u < 8; u++ )
{
for( v = 0; v < 8; v++ )
{
pDst[u][v] = (short int)ROUND(pSrc[u][v]/QU_LI[u][v]);
}
}
}

int CJpegImage::DToZ(short int pSrc[][8], short int pDstZ[])
{
TRACE("DTOZ:\n");
short int nDirect=1;//1左下角往右上角 -1右上角往左下角 0临界点
bool bPreDirect=true;//上一次的方向true左下角往右上角 false右上角往左下角
double pDst[8][8];
int i,j,x,y;
int nValidate=0;
x=0;y=0;
for(i=0;i<8;i++)
{
for(j=0;j<8;j++)
{
pDst[i][j]=pSrc[x][y];
// TRACE("%f ",pDst[i][j]);
if((x==0&&nDirect==1)||(x==7&&nDirect==-1))
{
if(y<7)
y++;
bPreDirect=(nDirect==1);
nDirect=0;
continue;
}
if((y==0&&nDirect==-1)||(y==7&&nDirect==1))
{
if(x<7)
x++;
bPreDirect=(nDirect==1);
nDirect=0;
continue;
}
if(nDirect==0&&bPreDirect)
nDirect=-1;
else if(nDirect==0&&!bPreDirect)
nDirect=1;
if(nDirect==1)
{
x--;
y++;
}
else if(nDirect=-1)
{
x++;
y--;
}
}
// TRACE("\n");
}
// memset(pBits,0,64);
for(i=7;i>=0;i--)
{
for(j=7;j>=0;j--)
{
if(nValidate==0)
{
if(pDst[i][j]!=0)
{
pDstZ[i*8+j]=(short)pDst[i][j];
// memset(pBits+i*8+j,(int)pDst[i][j],sizeof(unsigned char));
nValidate=1;
}
else
{
pDstZ[i*8+j]=0;
}
}
else
{
pDstZ[i*8+j]=(short)pDst[i][j];
// memset(pBits+i*8+j,(int)pDst[i][j],sizeof(unsigned char));
nValidate++;
}
}
}
return nValidate;
}

void CJpegImage::ZTOD(short int pSrcZ[], short pDst[][8])
{
TRACE("ZTOD:\n");
short int nDirect=1;//1左下角往右上角 -1右上角往左下角 0临界点
bool bPreDirect=true;//上一次的方向true左下角往右上角 false右上角往左下角
// double pDst[8][8];
int i,j,x,y;
int nValidate=0;
x=0;y=0;
for(i=0;i<8;i++)
{
for(j=0;j<8;j++)
{
pDst[x][y]=pSrcZ[i*8+j];
// TRACE("%f ",pDst[i][j]);
if((x==0&&nDirect==1)||(x==7&&nDirect==-1))
{
if(y<7)
y++;
bPreDirect=(nDirect==1);
nDirect=0;
continue;
}
if((y==0&&nDirect==-1)||(y==7&&nDirect==1))
{
if(x<7)
x++;
bPreDirect=(nDirect==1);
nDirect=0;
continue;
}
if(nDirect==0&&bPreDirect)
nDirect=-1;
else if(nDirect==0&&!bPreDirect)
nDirect=1;
if(nDirect==1)
{
x--;
y++;
}
else if(nDirect=-1)
{
x++;
y--;
}
}
// TRACE("\n");
}
}

void CJpegImage::IQuantify(short pSrc[][8], double pDst[][8])
{
TRACE("IQuantify:\n");
int u, v;
for( u = 0; u < 8; u++ )
{
for( v = 0; v < 8; v++ )
{
pDst[u][v] = (double)pSrc[u][v]*QU_LI[u][v];
}
}
}


void CJpegImage::IDCT(double pSrcDct[][8], short pDstDct[][8])
{
TRACE("IDCT:\n");
int x, y, u, v;
double persqrt2=1.0/sqrt(2.0);
TRACE("persqrt2=%f\n",persqrt2);
double cu,cv;
double pi = 3.14159265;
/* double tmp[8][8];

for( x = 0; x < 8; x++ )
{
for( y = 0; y < 8; y++ )
{
tmp[x][y]=pSrcDct[x][y];
}
}*/

for( x = 0; x < 8; x++ )
{
for( y = 0; y < 8; y++ )
{
m_pTemDct[x][y] = 0;


for( u = 0; u < 8; u++ )
{
for( v = 0; v < 8; v++ )
{
cu=(u==0)?persqrt2:1;
cv=(v==0)?persqrt2:1;
m_pTemDct[x][y] += cu*cv*pSrcDct[u][v] * cos((2.0*x+1.0)*u*pi/16.0)
*cos((2.0*y+1.0)*v*pi/16.0);
}
}

m_pTemDct[x][y] = m_pTemDct[x][y]/4.0;
}
}
for( x = 0; x < 8; x++ )
{
for( y = 0; y < 8; y++ )
{
pDstDct[x][y]=(short)m_pTemDct[x][y]+128;
}
}
}

DataBlock* CJpegImage::DataToBlock(DataBlock *pBlock)
{
return pBlock;
}

BOOL CJpegImage::Coding(CString strSource,CString strDest)
{
int i,j;
if(strDest=="")
strDest=strSource+".img";
CFile file;
CFileException fe;

DCTFileHeader dctHdr;
BITMAPFILEHEADER bmpHdr;
// LPBITMAPINFOHEADER *infoHdr;
TRACE("here\n");
if(!file.Open(strSource,CFile::modeRead|CFile::shareDenyWrite,&fe))
{
return FALSE;
}
CFile file2;
CFileException fe2;
if(!(file2.Open(strDest,CFile::modeCreate|
CFile::modeReadWrite|CFile::shareExclusive,&fe2)))
{
return FALSE;
}
if(file.Read((LPSTR)&bmpHdr,sizeof(BITMAPFILEHEADER))!=sizeof(BITMAPFILEHEADER))
{
return FALSE;
}
if(bmpHdr.bfType!=DIB_HEADER_MARKER)
return FALSE;
DWORD dwBitsSize=file.GetLength();
DWORD dwBSize=dwBitsSize-sizeof(BITMAPFILEHEADER);
HGLOBAL hDIB=::GlobalAlloc(GMEM_MOVEABLE|GMEM_ZEROINIT,dwBSize);
HGLOBAL hDCT=::GlobalAlloc(GMEM_MOVEABLE|GMEM_ZEROINIT,dwBSize);
if(hDIB==NULL||hDCT==NULL)
return FALSE;
LPSTR lpDIB=(LPSTR)::GlobalLock((HGLOBAL)hDIB);
LPSTR lpDCT=(LPSTR)::GlobalLock((HGLOBAL)hDIB);
if(file.ReadHuge(lpDIB,dwBSize)!=dwBSize)
{
::GlobalUnlock((HGLOBAL)hDIB);
::GlobalFree((HGLOBAL)hDIB);
return FALSE;
}
LPBITMAPINFOHEADER lpInfoHdr=(LPBITMAPINFOHEADER)lpDIB;
lHeight=lpInfoHdr->biHeight;
lWidth=lpInfoHdr->biWidth;
TRACE("Width=%d,Height=%d\n",lWidth,lHeight);
int xBlockCount=(int)((double)lWidth/8.0+0.5);
int yBlockCount=(int)((double)lHeight/8.0+0.5);
dctHdr.xBlockCount=xBlockCount;
dctHdr.yBlockCount=yBlockCount;
dctHdr.dtHdrSize=sizeof(dctHdr);
dctHdr.dtType=DCT_HEADER_MARKER;
dctHdr.lHeight=lHeight;
dctHdr.lWidth=lWidth;
file2.Write((LPSTR)&dctHdr,sizeof(dctHdr));
for(i=0;i<xBlockCount;i++)
{
for(j=0;j<yBlockCount;j++)
{
}
}
::GlobalUnlock((HGLOBAL)hDIB);
file.Close();
return true;
}

BOOL CJpegImage::Decode(CString strSource, CString strDest)
{
return true;
}