小波变换和恢复1（输入数组的形式）

抽样差值的小波变换和恢复（对一层小波效果好，多层小波恢复的效果较差）

#include <stdlib.h>
#include <cv.h>
#include <highgui.h>

// 二维离散小波变换（矩阵数组）
void DWT(float *pMatrix, int height, int width, int nLayer);
// 二维离散小波恢复（矩阵数组）
void IDWT(float *pMatrix, int height, int width, int nLayer);

int main()
{
       // 小波变换层数
       int nLayer = 1;
       // 输入彩色图像
       IplImage *pSrc = cvLoadImage("12.bmp", CV_LOAD_IMAGE_GRAYSCALE);
       //cvSaveImage("222.bmp",pSrc);
       // 计算小波图象大小
       int height = pSrc->height;
       int width = pSrc->width;
       // 创建小波图象
       cvNamedWindow("1",1);
       cvShowImage("1",pSrc);
       //cvWaitKey(0);
       IplImage *pWavelet = cvCreateImage(cvGetSize(pSrc), IPL_DEPTH_32F, 1);
       IplImage *image = cvCreateImage(cvGetSize(pSrc), IPL_DEPTH_8U, 1);
       if (pWavelet)
       {

#if 0   //分块小波变换 
       int nwidth = 8;
        float *pMatrix1 = new float[height*(nwidth+2)];
        for(int i = 0; i < height*(nwidth+2); i++)
        {
               pMatrix1[i] = 255.0f;
        }
        int time = (int)(width / nwidth);
        for(int n = 0; n < time; n++)
         {
               float *pMatrix2 = new float[height*2];
               int temp = 0;
               int temp1 = 0;
               for(int i = 0; i < height; i++)
               {
                      for(int j = n*nwidth; j < (n+1)*nwidth; j++)
                       {
                             pMatrix1[i*nwidth + j%nwidth] = 0.0;
                             temp1 = (unsigned char)(*(pSrc->imageData+i*pSrc->widthStep+j));
                             pMatrix1[i*nwidth + j%nwidth] = (float)(temp1);
                      }
               }
               for(int i = 0; i < height; i++)
               {
                      for(int j =(n+1)*nwidth-2; j < (n+1)*nwidth; j++)
                      {
                             pMatrix1[i*nwidth + j%nwidth] = 0.0;
                             temp1 = (unsigned char)(*(pSrc->imageData+i*pSrc->widthStep+j));
                             pMatrix1[i*nwidth + j%nwidth] = (float)(temp1);
                      }
               }

               {
                      // 二维离散小波变换
                      DWT(pMatrix1,height,nwidth, nLayer);
                      // 二维离散小波恢复
                      IDWT(pMatrix1,height,nwidth, nLayer);
               }

               for(int i = 0; i < height; i++)
               {
                      for(int j = n*nwidth; j < (n+1)*nwidth; j++)
                      {
                             temp = (int)(pMatrix1[i*nwidth + j%nwidth]);
                             *(image->imageData+i*image->widthStep+j) = (unsigned char)temp;
                      } 
               }
               delete pMatrix2;
        }
        delete pMatrix1;

#else  //整幅图的小波变换
        float *pMatrix = new float[height*width];
        int temp = 0;
        int temp1 = 0;
        for(int i = 0; i < height; i++)
        {
               for(int j = 0; j < width; j++)
               {
                      pMatrix[i*width + j] = 0.0;
                      temp1 = (unsigned char)(*(pSrc->imageData+i*pSrc->widthStep+j));
                      pMatrix[i*width + j] = (float)(temp1);
               }
        }

        {
                二维离散小波变换
               //DWT(pMatrix,height,width, nLayer);
                二维离散小波恢复
               //IDWT(pMatrix,height,width, nLayer);

               // 二维离散小波变换
               dwt(pMatrix,height,width, nLayer);
               // 二维离散小波恢复
               idwt(pMatrix,height,width, nLayer);
        }
  
         for(int i = 0; i < height; i++)
        {
               for(int j = 0; j < width; j++)
               {
                      temp = (int)(pMatrix[i*width + j]);
                      if(temp < 0)
                             temp = 0;
                      if(temp > 255)
                             temp = 255;
                      *(image->imageData+i*image->widthStep+j) = (unsigned char)temp;
               } 
        }
        delete pMatrix;
#endif
        cvNamedWindow("123",1);
        cvShowImage("123",image);
        cvWaitKey(0);
        cvReleaseImage(&pWavelet);

  
       }
       // 显示图像pSrc
       // ...
       cvReleaseImage(&pSrc);

       return 0;
}

// 二维离散小波变换（单通道浮点图像）
void DWT(float *pMatrix, int height, int width, int nLayer)
{
 int     i, x, y, n;
 float   fValue   = 0;
 float   fRadius  = sqrt(2.0f);
 int     nWidth   = width;
 int     nHeight  = height;
 int     nHalfW   = nWidth / 2;
 int     nHalfH   = nHeight / 2;
 float **pData    = new float*[height];
 float  *pRow     = new float[width];
 float  *pColumn  = new float[height];
 for (i = 0; i < height; i++)
 {
  pData[i] = &pMatrix[i*width];
 }
 // 多层小波变换
 for (n = 0; n < nLayer; n++, nWidth /= 2, nHeight /= 2, nHalfW /= 2, nHalfH /= 2)
 {
  // 水平变换
  for (y = 0; y < nHeight; y++)
  {
   // 奇偶分离
   memcpy(pRow, pData[y], sizeof(float) * nWidth);
   for (i = 0; i < nHalfW; i++)
   {
    x = i * 2;
    pData[y][i] = pRow[x];
    pData[y][nHalfW + i] = pRow[x + 1];
   }
   //差
   for (i = 0; i < nHalfW - 1; i++)
   {
    fValue = (pData[y][i] + pData[y][i + 1]) / 2;
    pData[y][nHalfW + i] -= fValue;
   }
   fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
   pData[y][nWidth - 1] -= fValue;
   fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
   pData[y][0] += fValue;
   for (i = 1; i < nHalfW; i++)
   {
    fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
    pData[y][i] += fValue;
   }
   // 频带系数
   for (i = 0; i < nHalfW; i++)
   {
    pData[y][i] *= fRadius;
    pData[y][nHalfW + i] /= fRadius;
   }
  }
  // 垂直变换
  for (x = 0; x < nWidth; x++)
  {
   // 奇偶分离
   for (i = 0; i < nHalfH; i++)
   {
    y = i * 2;
    pColumn[i] = pData[y][x];
    pColumn[nHalfH + i] = pData[y + 1][x];
   }
   for (i = 0; i < nHeight; i++)
   {
    pData[i][x] = pColumn[i];
   }
   //差
   for (i = 0; i < nHalfH - 1; i++)
   {
    fValue = (pData[i][x] + pData[i + 1][x]) / 2;
    pData[nHalfH + i][x] -= fValue;
   }
   fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
   pData[nHeight - 1][x] -= fValue;
   fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
   pData[0][x] += fValue;
   for (i = 1; i < nHalfH; i++)
   {
    fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
    pData[i][x] += fValue;
   }
   // 频带系数
   for (i = 0; i < nHalfH; i++)
   {
    pData[i][x] *= fRadius;
    pData[nHalfH + i][x] /= fRadius;
   }
  }
 }
 delete[] pData;
 delete[] pRow;
 delete[] pColumn;
}

// 二维离散小波恢复（单通道浮点图像）
void IDWT(float *pMatrix, int height, int width, int nLayer)
{
 int     i, x, y, n;
 float   fValue   = 0;
 float   fRadius  = sqrt(2.0f);
 int     nWidth   = width >> (nLayer - 1);
 int     nHeight  = height >> (nLayer - 1);
 int     nHalfW   = nWidth / 2;
 int     nHalfH   = nHeight / 2;
 float **pData    = new float*[height];
 float  *pRow     = new float[width];
 float  *pColumn  = new float[height];
 for (i = 0; i < height; i++)
 {
  pData[i] = &pMatrix[i*width];
 }
 // 多层小波恢复
 for (n = 0; n < nLayer; n++, nWidth *= 2, nHeight *= 2, nHalfW *= 2, nHalfH *= 2)
 {
  // 垂直恢复
  for (x = 0; x < nWidth; x++)
  {
   // 频带系数
   for (i = 0; i < nHalfH; i++)
   {
    pData[i][x] /= fRadius;
    pData[nHalfH + i][x] *= fRadius;
   }
   // 提升小波恢复
   fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
   pData[0][x] -= fValue;
   for (i = 1; i < nHalfH; i++)
   {
    fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
    pData[i][x] -= fValue;
   }
   for (i = 0; i < nHalfH - 1; i++)
   {
    fValue = (pData[i][x] + pData[i + 1][x]) / 2;
    pData[nHalfH + i][x] += fValue;
   }
   fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
   pData[nHeight - 1][x] += fValue;
   // 奇偶合并
   for (i = 0; i < nHalfH; i++)
   {
    y = i * 2;
    pColumn[y] = pData[i][x];
    pColumn[y + 1] = pData[nHalfH + i][x];
   }
   for (i = 0; i < nHeight; i++)
   {
    pData[i][x] = pColumn[i];
   }
  }
  // 水平恢复
  for (y = 0; y < nHeight; y++)
  {
   // 频带系数
   for (i = 0; i < nHalfW; i++)
   {
    pData[y][i] /= fRadius;
    pData[y][nHalfW + i] *= fRadius;
   }
    //提升小波恢复
   fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
   pData[y][0] -= fValue;
   for (i = 1; i < nHalfW; i++)
   {
    fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
    pData[y][i] -= fValue;
   }
   for (i = 0; i < nHalfW - 1; i++)
   {
    fValue = (pData[y][i] + pData[y][i + 1]) / 2;
    pData[y][nHalfW + i] += fValue;
   }
   fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
   pData[y][nWidth - 1] += fValue;
   // 奇偶合并
   for (i = 0; i < nHalfW; i++)
   {
    x = i * 2;
    pRow[x] = pData[y][i];
    pRow[x + 1] = pData[y][nHalfW + i];
   }
   memcpy(pData[y], pRow, sizeof(float) * nWidth);
  }
 }
 delete[] pData;
 delete[] pRow;
 delete[] pColumn;
}