DCT高速変換

12398 ワード
グラフィック画像
DCT高速変換著者:陳祖尚(blog:darnshong)ダウンロードソース1、引用DCT変換はデジタル画像処理において重要な変換であり、多くの重要な画像アルゴリズム、画像応用はJPEG画像符号化方式のようなDCT変換に基づいている.大きなサイズの2次元数値マトリクスでは,通常のDCT変換を用いると,その時間は耐えられず,実用化にも至らない.この難点を克服するには、DCT変換の高速アルゴリズムは非常に魅力的にほかならない.現在のところ、DCT変換の高速アルゴリズムは以下の2つの方式にほかならない:1.FFTアルゴリズムの一般的な採用のため、直接FFTを利用してDCT変換を実現する高速アルゴリズムは比較的に容易である.しかし、この方法には、計算プロセスが複数の演算にかかわるという不足もあります.DCT変換前後のデータはすべて実数であるため、計算中に複数を導入し、1対の複素の加算は2対の実数の加算に相当し、1対の複素の乗算は4対の実数の乗算と2対の実数の加算に相当し、明らかに演算量を増加させ、ハードウェアストレージにもより高い要求を提出した.2.実数領域で直接DCT高速変換を行う.明らかに、この方法は、前の方法に比べて、計算量およびハードウェア要件が前者よりも優れている.これに鑑みて,本論文では,DCT変換の高速アルゴリズムを実現するために第2の方法を採用した.二、理論の導出は紙面に限られ、ここで羅列することはできない.具体的な導出過程は「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変換プロセスは2つの部分に分けることができる:順方向追跡と後方向回帰ルート順方向追跡:
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;
}
完全に1次元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];
}
1次元DCT変換を利用して2次元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変換プロセスも2つの部分に分けることができる:順方向追跡と後方向回帰ルート前追底
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変換を実行するためにスレッドを1つ増やさなければ、動作インタフェースは、プログラムがDCT変換の計算に追われて応答を失う可能性があるため、DCT変換を行うためのスレッドを1つ増やす必要がある.まず、構造
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正変換するには94 ms、IDCT逆変換も94 msしかかからないが、一般DCT変換を用いると、所要時間は575172 msです.このことから,DCT高速変換の大きな利点は,計算速度が速く,効率が高いことである.http://www.vckbase.com/document/viewdoc/?id=1511
Talibで使えるローソク足についての関数まとめ(4)
Restful APIとGraphQLの違いは何ですか?!