Chinaunix首页 | 论坛 | 博客
  • 博客访问: 1218105
  • 博文数量: 950
  • 博客积分: 10000
  • 博客等级: 上将
  • 技术积分: 13070
  • 用 户 组: 普通用户
  • 注册时间: 2008-08-04 09:23
文章分类

全部博文(950)

文章存档

2011年(1)

2008年(949)

我的朋友

分类: C/C++

2008-08-04 09:35:25

下载本文示例代码
下载源代码

一、引言

DCT变换是数字图像处理中重要的变换,很多重要的图像算法、图像应用都是基于DCT变换的,如JPEG图像编码方式。对于大尺寸的二维数值矩阵,倘若采用普通的DCT变换来进行,其所花费的时间将是让人难以忍受甚至无法达到实用。而要克服这一难点,DCT变换的快速算法无非是非常吸引人的。

就目前而言,DCT变换的快速算法无非有以下两种方式:

1.由于FFT算法的普便采用,直接利用FFT来实现DCT变换的快速算法相比来说就相对容易。但是此种方法也有不足:计算过程会涉及到复数的运算。由于DCT变换前后的数据都是实数,计算过程中引入复数,而一对复数的加法相当于两对实数的加法,一对复数的乘法相当于四对实数的乘法和两对实数的加法,显然是增加了运算量,也给硬件存储提出了更高的要求。

2.直接在实数域进行DCT快速变换。显然,这种方法相比于前一种而言,计算量和硬件要求都要优于前者。

鉴于此,本文采用第二种方法来实现DCT变换的快速算法。

二、理论推导

  限于篇幅,在此不能罗列,具体推导过程可参见《DCT快速新算法及滤波器结构研究与子波变换域图像降噪研究》华南理工大学博士论文。

三、程序实现

DCT快速变换
考虑到DCT变换中的系数要重复计算,可使用查找表来提高运行的效率,只要开始DCT变换之前计算一次,DCT变换中就可以只查找而无需计算系数。

void initDCTParam(int deg)

{

      // deg 为DCT变换数据长度的幂



      if(bHasInit)

      {

             return; //不用再计算查找表

      }



      int total, halftotal, i, group, endstart, factor;



      total = 1 << deg;



      if(C != NULL) delete []C;



      C = (double *)new double[total];



      halftotal = total >> 1;



      for(i=0; i < halftotal; i  )

             C[total-i-1]=(double)(2*i 1);



      for(group=0; group < deg-1; group  )

      { 



             endstart=1 << (deg-1-group);



             int len = endstart >> 1;



             factor=1 << (group 1);



             for(int j = 0;j < len; j  )

                    C[endstart-j-1] = factor*C[total-j-1];

      }



      for(i=1; i < total; i  )

             C[i] = 2.0*cos(C[i]*PI/(total << 1)); ///C[0]空着,没使用



      bHasInit=true;

}

DCT变换过程可分为两部分:前向追底和后向回根

前向追底:

void dct_forward(double *f,int deg)

{

      // f中存储DCT数据



      int i_deg, i_halfwing, total, wing, wings, winglen, halfwing;



      double temp1,temp2;



      total = 1 << deg;



      for(i_deg = 0; i_deg < deg; i_deg  )

      {

             wings = 1 << i_deg;

             winglen = total >> i_deg;

             halfwing = winglen >> 1;

             for(wing = 0; wing < wings; wing  )

             {

                    for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing  )

                    {

                           temp1 = f[wing*winglen i_halfwing];

                           temp2 = f[(wing 1)*winglen-1-i_halfwing];

                           if(wing%2)

                                  swap(temp1,temp2); // 交换temp1与temp2的值



                           f[wing*winglen i_halfwing] = temp1 temp2;

                           f[(wing 1)*winglen-1-i_halfwing] = 

                                (temp1-temp2)*C[winglen-1-i_halfwing];

                    }

             }

      }

}

后向回根:
void dct_backward(double *f,int deg)

{

      // f中存储DCT数据



      int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;



      total = 1 << deg;



      for(i_deg = deg-1; i_deg >= 0; i_deg--)

      {

             wings = 1 << i_deg;

             winglen = 1 << (deg-i_deg);



             halfwing = winglen >> 1;



             for(wing = 0; wing < wings; wing  )

             {

                    for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing  )

                    {  

                           //f[i_halfwing wing*winglen] = f[i_halfwing wing*winglen];

                           if(i_halfwing == 0)

                           {

                                    f[halfwing wing*winglen i_halfwing] = 

                                        0.5*f[halfwing wing*winglen i_halfwing];

                            }

                           else

                           {

                                  temp1=bitrev(i_halfwing,deg-i_deg-1);   // bitrev为位反序

                                  temp2=bitrev(i_halfwing-1,deg-i_deg-1); // 第一参数为要变换的数

                     // 第二参数为二进制长度

                                  f[halfwing wing*winglen temp1] =

                                       f[halfwing wing*winglen temp1]-f[halfwing wing*winglen temp2];

                           }     

                    }

             }

      }

}



位反序函数如下:
int bitrev(int bi,int deg)

{    

      int j = 1, temp = 0, degnum, halfnum;



      degnum = deg;



      //if(deg<0) return 0;



      if(deg == 0) return bi;



      halfnum = 1 << (deg-1);



      while(halfnum)

      {            

             if(halfnum&bi)

                    temp  = j;



             j<<=1;



             halfnum >>= 1;

      }



      return temp;

}

完整实现一维DCT变换:
void fdct_1D_no_param(double *f,int deg)

{

      initDCTParam(deg);

      dct_forward(f,deg);

      dct_backward(f,deg);

      fbitrev(f,deg);     // 实现位反序排列

      f[0] = 1/(sqrt(2.0))*f[0];

}



void fdct_1D(double *f,int deg)

{

      fdct_1D_no_param(f,deg);

      int total = 1 << deg;

      double param = sqrt(2.0/total);

      for(int i = 0; i < total; i  )

             f[i] = param*f[i];

}



利用一维DCT变换来实现二维DCT变换:
void fdct_2D(double *f,int deg_row,int deg_col)

{    

      int rows,cols,i_row,i_col;

      double two_div_sqrtcolrow;

      rows=1 << deg_row;

      cols=1 << deg_col;

      init2D_Param(rows,cols);

      two_div_sqrtcolrow = 2.0/(sqrt(double(rows*cols)));  



      for(i_row = 0; i_row < rows; i_row  )

      {

             fdct_1D_no_param(f i_row*cols,deg_col);

      }



      for(i_col = 0; i_col < cols; i_col  )

      {

             for(i_row = 0; i_row < rows; i_row  )

             {

                    temp_2D[i_row] = f[i_row*cols i_col];

             }



             fdct_1D_no_param(temp_2D, deg_row);



             for(i_row = 0; i_row < rows; i_row  )

             {

                    f[i_row*cols i_col] = temp_2D[i_row]*two_div_sqrtcolrow;

             }          

      }



      bHasInit = false;

}



IDCT快速变换
初始化查找表:
void initIDCTParam(int deg)

{

      if(bHasInit)

             return;    //不用再计算查找表



      int total, halftotal, i, group, endstart, factor;

      total = 1 << deg;



      // if(C!=NULL) delete []C;

      // C=(double *)new double[total];



      // 由于正变换已经为C申请了空间,所以逆变换就需再申请空间了!



      halftotal = total >> 1;



      for(i = 0; i < halftotal; i  )

             C[total-i-1] = (double)(2*i 1);



      for(group = 0; group < deg-1; group  )

      { 

             endstart = 1 << (deg-1-group);

             int len = endstart>>1;

             factor = 1 << (group 1);

             for(int j = 0; j < len; j  )

                    C[endstart-j-1] = factor*C[total-j-1];

      }



      for(i = 1; i < total; i  )

             C[i] = 1.0/(2.0*cos(C[i]*PI/(total << 1)));       // C[0]空着没用



      bHasInit=true;

}

IDCT变换过程也可分为两部分:前向追底和后向回根
前向追底
void idct_forward(double *F,int deg)

{

      int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;



      total = 1 << deg;

      for(i_deg = 0; i_deg < deg; i_deg  )

      {

             wings = 1 << i_deg;

             winglen = 1 << (deg-i_deg);

             halfwing = winglen >> 1;

             for(wing = 0; wing < wings; wing  )

             {

                    for(i_halfwing = halfwing-1; i_halfwing >= 0; i_halfwing--)

                    {

                           if(i_halfwing == 0)

                           {

                                  F[halfwing wing*winglen i_halfwing] = 

                                    2.0*F[halfwing wing*winglen i_halfwing];

                            }

                           else

                           { 

                                  temp1 = bitrev(i_halfwing,deg-i_deg-1);

                                  temp2 = bitrev(i_halfwing-1,deg-i_deg-1);

                                  F[halfwing wing*winglen temp1] = F[halfwing wing*winglen temp1]

                                           F[halfwing wing*winglen temp2];

                           }

                    }

             }

      }

}

后向回根
void idct_backward(double *F, int deg)

{

      int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;



      double temp1, temp2;

      total = 1 << deg;

      for(i_deg = deg-1; i_deg >= 0; i_deg--)

      {

             wings = 1 << i_deg;

             winglen = total >> i_deg;

             halfwing = winglen >> 1;

             for(wing = 0; wing < wings; wing  )

             {

                    for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing  )

                    {

                           temp1 = F[wing*winglen i_halfwing];

                           temp2 = F[(wing 1)*winglen-1-i_halfwing]*C[winglen-1-i_halfwing];

                           if(wing % 2)

                           {

                                  F[wing*winglen i_halfwing] = (temp1-temp2)*0.5;

                                  F[(wing 1)*winglen-1-i_halfwing] = (temp1 temp2)*0.5;

                           }

                           else

                           {

                                  F[wing*winglen i_halfwing] = (temp1 temp2)*0.5;

                                  F[(wing 1)*winglen-1-i_halfwing] = (temp1-temp2)*0.5;

                           }

                    }

             }

      }

}

完整实现一维IDCT变换:
void fidct_1D_no_param(double *F, int deg)

{

      initIDCTParam(deg);

      F[0] = F[0]*sqrt(2.0);

      fbitrev(F, deg);

      idct_forward(F, deg);

      idct_backward(F, deg);

}



void fdct_1D(double *f, int deg)

{

      fdct_1D_no_param(f, deg);

      int total = 1 << deg;



      double param = sqrt(2.0/total);

      for(int i = 0; i < total; i  )

             f[i] = param*f[i];

}



利用一维IDCT变换来实现二维IDCT变换:
void fidct_2D(double *F, int deg_row, int deg_col)

{

      int rows,cols,i_row,i_col;



      double     sqrtcolrow_div_two;

      rows = 1 << deg_row;

      cols = 1 << deg_col;

      init2D_Param(rows,cols);

      sqrtcolrow_div_two = (sqrt(double(rows*cols)))/2.0;



      for(i_row = 0; i_row < rows; i_row  )

      {

             fidct_1D_no_param(F i_row*cols,deg_col);

      }



      for(i_col = 0; i_col < cols; i_col  )

      {

             for(i_row = 0; i_row < rows; i_row  )

             {

                    temp_2D[i_row] = F[i_row*cols i_col];

             }



             fidct_1D_no_param(temp_2D, deg_row);

             for(i_row = 0; i_row < rows; i_row  )

             {

                    F[i_row*cols i_col] = temp_2D[i_row]*sqrtcolrow_div_two;

             }

      }



      bHasInit=false;

}



多线程的考量 由于DCT变换要花费一定的时间,特别是在数据矩阵尺寸比较大的时候。此时,如果没有增加一个线程来执行DCT变换,操作界面可能因程序忙于DCT变换的计算而失去响应,所以,增加一个用来进行DCT变换的线程是十分必要的。
首先定义一个结构
typedef struct

{    

      int row;

      int col;

      double *data;

      //double *data2;

      //double *data3; // 在计算彩色图象的数据矩阵时,彩色图象用RGB三个分量



      bool m_bfinished;



      DWORD m_start,m_end; //以毫秒计,用来计算DCT变换所用的时间;

}RUNINFO;
DCT变换进程函数:
UINT ThreadProcfastDct(LPVOID pParam)

{

      RUNINFO *pinfo = (RUNINFO*)pParam;

      pinfo->m_start = ::GetTickCount();

      fdct_2D((double *)pinfo->data, GetTwoIndex(pinfo->row), GetTwoIndex(pinfo->col));

      pinfo->m_end = ::GetTickCount();

      pinfo->m_bfinished = true;



      return 1;

}
IDCT变换进程函数:
UINT ThreadProcfastIDct(LPVOID pParam)

{

      RUNINFO *pinfo = (RUNINFO*)pParam;     

      pinfo->m_start = ::GetTickCount();

      fidct_2D((double *)pinfo->data, GetTwoIndex(pinfo->row), GetTwoIndex(pinfo->col));

      pinfo->m_end = ::GetTickCount();

      pinfo->m_bfinished = true;



      return 1;

}


四、程序运行

图1 普通DCT变换


图2 快速DCT变换


图3 快速IDCT变换

从以上可以看出,采用上述快速DCT变换对一幅256灰度的256*256的图像进行DCT正变换只需94ms,IDCT逆变换也只需94ms,而如果采用普通DCT变换,所需时间要575172ms。由此可见,DCT快速变换的巨大的优势,计算速度快,效率高。 下载本文示例代码
阅读(849) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~