Set 4 8 16 32 sightings of fast DCT algorithm...

#include "stdafx.h"

#include "GlobalApi.h"

#include "cdib.h"

#include "math.h"

#include <direct.h>

#include <complex>

#include <afx.h>

//#include "Afxtempl.h"

using namespace std;

/*************************************************************************

 *

 * 函数名称：

 *   MakeGauss()

 *

 * 输入参数:

 *   double sigma         - 高斯函数的标准差

 *   double **pdKernel - 指向高斯数据数组的指针

 *   int *pnWindowSize - 数据的长度

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   这个函数可以生成一个一维的高斯函数的数字数据，理论上高斯数据的长度应

 *   该是无限长的，但是为了计算的简单和速度，实际的高斯数据只能是有限长的

 *   pnWindowSize就是数据长度

 *   

 *************************************************************************

 */

void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)

{

// 循环控制变量

int i   ;

// 数组的中心点

int nCenter;

// 数组的某一点到中心点的距离

double  dDis  ; 

double PI = 3.14159;

// 中间变量

double  dValue; 

double  dSum  ;

dSum = 0 ; 

// 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。

// 这些数据会覆盖绝大部分的滤波系数

*pnWindowSize = 1 + 2 * ceil(3 * sigma);

// 中心

nCenter = (*pnWindowSize) / 2;

// 分配内存

*pdKernel = new double[*pnWindowSize] ;

for(i=0; i< (*pnWindowSize); i++)

{

dDis = (double)(i - nCenter);

dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma)) / (sqrt(2 * PI) * sigma );

(*pdKernel)[i] = dValue ;

dSum += dValue;

}

// 归一化

for(i=0; i<(*pnWindowSize) ; i++)

{

(*pdKernel)[i] /= dSum;

}

}

/*************************************************************************

 *

 * 函数名称：

 *   GaussianSmooth()

 *

 * 输入参数:

 *   unsigned char * pUnchImg - 指向图象数据的指针

 *   int nWidth - 图象数据宽度

 *   int nHeight - 图象数据高度

 *   double dSigma - 高斯函数的标准差

 *   unsigned char * pUnchSmthdImg - 指向经过平滑之后的图象数据

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   为了抑止噪声，采用高斯滤波对图象进行滤波，滤波先对x方向进行，然后对

 *   y方向进行。

 *   

 *************************************************************************

 */

void GaussianSmooth(unsigned char *pUnchImg, int nWidth, int nHeight, 

double sigma, unsigned char * pUnchSmthdImg)

{

// 循环控制变量

    int y;

int x;

int i;

// 高斯滤波器的数组长度

int nWindowSize;

//  窗口长度的1/2

int nHalfLen;   

// 一维高斯数据滤波器

double *pdKernel ;

// 高斯系数与图象数据的点乘

double  dDotMul     ;

// 高斯滤波系数的总和

double  dWeightSum     ;          

  

// 中间变量

double * pdTmp ;

// 分配内存

pdTmp = new double[nWidth*nHeight];

// 产生一维高斯数据滤波器

// MakeGauss(sigma, &dKernel, &nWindowSize);

MakeGauss(sigma, &pdKernel, &nWindowSize) ;

// MakeGauss返回窗口的长度，利用此变量计算窗口的半长

nHalfLen = nWindowSize / 2;

  // x方向进行滤波

for(y=0; y<nHeight; y++)

{

for(x=0; x<nWidth; x++)

{

dDotMul = 0;

dWeightSum = 0;

for(i=(-nHalfLen); i<=nHalfLen; i++)

{

// 判断是否在图象内部

if( (i+x) >= 0  && (i+x) < nWidth )

{

dDotMul += (double)pUnchImg[y*nWidth + (i+x)] * pdKernel[nHalfLen+i];

dWeightSum += pdKernel[nHalfLen+i];

}

}

pdTmp[y*nWidth + x] = dDotMul/dWeightSum ;

}

}

// y方向进行滤波

for(x=0; x<nWidth; x++)

{

for(y=0; y<nHeight; y++)

{

dDotMul = 0;

dWeightSum = 0;

for(i=(-nHalfLen); i<=nHalfLen; i++)

{

// 判断是否在图象内部

if( (i+y) >= 0  && (i+y) < nHeight )

{

dDotMul += (double)pdTmp[(y+i)*nWidth + x] * pdKernel[nHalfLen+i];

dWeightSum += pdKernel[nHalfLen+i];

}

}

pUnchSmthdImg[y*nWidth + x] = (unsigned char)(int)dDotMul/dWeightSum ;

}

}

// 释放内存

delete []pdKernel;

pdKernel = NULL ;

delete []pdTmp;

pdTmp = NULL;

}

/*************************************************************************

 *

 * 函数名称：

 *   DirGrad()

 *

 * 输入参数:

 *   unsigned char *pUnchSmthdImg - 经过高斯滤波后的图象

 *   int nWidht - 图象宽度

 *   int nHeight       - 图象高度

 *   int *pnGradX                         - x方向的方向导数

 *   int *pnGradY                         - y方向的方向导数

 * 返回值:

 *   无

 *

 * 说明:

 *   这个函数计算方向倒数，采用的微分算子是(-1 0 1) 和 (-1 0 1)'(转置)

 *   计算的时候对边界象素采用了特殊处理

 *   

 *   

 *************************************************************************

 */

void DirGrad(unsigned char *pUnchSmthdImg, int nWidth, int nHeight,

 int *pnGradX , int *pnGradY)

{

// 循环控制变量

int y ;

int x ;

// 计算x方向的方向导数，在边界出进行了处理，防止要访问的象素出界

for(y=0; y<nHeight; y++)

{

for(x=0; x<nWidth; x++)

{

pnGradX[y*nWidth+x] = (int) ( pUnchSmthdImg[y*nWidth+min(nWidth-1,x+1)]  

          - pUnchSmthdImg[y*nWidth+max(0,x-1)]     );

}

}

// 计算y方向的方向导数，在边界出进行了处理，防止要访问的象素出界

for(x=0; x<nWidth; x++)

{

for(y=0; y<nHeight; y++)

{

pnGradY[y*nWidth+x] = (int) ( pUnchSmthdImg[min(nHeight-1,y+1)*nWidth + x]  

- pUnchSmthdImg[max(0,y-1)*nWidth+ x ]     );

}

}

}

/*************************************************************************

 *

 * 函数名称：

 *   GradMagnitude()

 *

 * 输入参数:

 *   int *pnGradX                         - x方向的方向导数

 *   int *pnGradY                         - y方向的方向导数

 *   int nWidht - 图象宽度

 *   int nHeight       - 图象高度

 *   int *pnMag                           - 梯度幅度   

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   这个函数利用方向倒数计算梯度幅度，方向倒数是DirGrad函数计算的结果

 *   

 *************************************************************************

 */

void GradMagnitude(int *pnGradX, int *pnGradY, int nWidth, int nHeight, int *pnMag)

{

// 循环控制变量

int y ;

int x ;

// 中间变量

double dSqtOne;

double dSqtTwo;

for(y=0; y<nHeight; y++)

{

for(x=0; x<nWidth; x++)

{

dSqtOne = pnGradX[y*nWidth + x] * pnGradX[y*nWidth + x];

dSqtTwo = pnGradY[y*nWidth + x] * pnGradY[y*nWidth + x];

pnMag[y*nWidth + x] = (int)(sqrt(dSqtOne + dSqtTwo) + 0.5);

}

}

}

/*************************************************************************

 *

 * 函数名称：

 *   NonmaxSuppress()

 *

 * 输入参数:

 *   int *pnMag                - 梯度图

 *   int *pnGradX  - x方向的方向导数

 *   int *pnGradY              - y方向的方向导数

 *   int nWidth                - 图象数据宽度

 *   int nHeight               - 图象数据高度

 *   unsigned char *pUnchRst   - 经过NonmaxSuppress处理后的结果

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   抑止梯度图中非局部极值点的象素。

 *   

 *************************************************************************

 */

void NonmaxSuppress(int *pnMag, int *pnGradX, int *pnGradY, int nWidth, 

int nHeight, unsigned char *pUnchRst)

{

// 循环控制变量

int y ;

int x ;

int nPos;

// x方向梯度分量

int gx  ;

int gy  ;

// 临时变量

int g1, g2, g3, g4 ;

double weight  ;

double dTmp1   ;

double dTmp2   ;

double dTmp    ;

// 设置图象边缘部分为不可能的边界点

for(x=0; x<nWidth; x++)

{

pUnchRst[x] = 0 ;

pUnchRst[nHeight-1+x] = 0;

}

for(y=0; y<nHeight; y++)

{

pUnchRst[y*nWidth] = 0 ;

pUnchRst[y*nWidth + nWidth-1] = 0;

}

for(y=1; y<nHeight-1; y++)

{

for(x=1; x<nWidth-1; x++)

{

nPos = y*nWidth + x ;

// 如果当前象素的梯度幅度为0，则不是边界点

if(pnMag[nPos] == 0 )

{

pUnchRst[nPos] = 0 ;

}

else

{

// 当前象素的梯度幅度

dTmp = pnMag[nPos] ;

// x，y方向导数

gx = pnGradX[nPos]  ;

gy = pnGradY[nPos]  ;

// 如果方向导数y分量比x分量大，说明导数的方向更加“趋向”于y分量。

if (abs(gy) > abs(gx)) 

{

// 计算插值的比例

weight = fabs(gx)/fabs(gy); 

g2 = pnMag[nPos-nWidth] ; 

g4 = pnMag[nPos+nWidth] ;

// 如果x，y两个方向的方向导数的符号相同

// C是当前象素，与g1-g4的位置关系为：

// g1 g2 

//  C         

//  g4 g3 

if (gx*gy > 0) 

{ 

g1 = pnMag[nPos-nWidth-1] ;

g3 = pnMag[nPos+nWidth+1] ;

} 

// 如果x，y两个方向的方向导数的符号相反

// C是当前象素，与g1-g4的位置关系为：

//    g2 g1

//  C         

// g3 g4  

else 

{ 

g1 = pnMag[nPos-nWidth+1] ;

g3 = pnMag[nPos+nWidth-1] ;

} 

}

// 如果方向导数x分量比y分量大，说明导数的方向更加“趋向”于x分量

// 这个判断语句包含了x分量和y分量相等的情况

else

{

// 计算插值的比例

weight = fabs(gy)/fabs(gx); 

g2 = pnMag[nPos+1] ; 

g4 = pnMag[nPos-1] ;

// 如果x，y两个方向的方向导数的符号相同

// C是当前象素，与g1-g4的位置关系为：

// g3   

// g4 C g2       

//       g1

if (gx*gy > 0) 

{

g1 = pnMag[nPos+nWidth+1] ;

g3 = pnMag[nPos-nWidth-1] ;

} 

// 如果x，y两个方向的方向导数的符号相反

// C是当前象素，与g1-g4的位置关系为：

//      g1

// g4 C g2       

//  g3     

else 

{ 

g1 = pnMag[nPos-nWidth+1] ;

g3 = pnMag[nPos+nWidth-1] ;

}

}

// 下面利用g1-g4对梯度进行插值

{

dTmp1 = weight*g1 + (1-weight)*g2 ;

dTmp2 = weight*g3 + (1-weight)*g4 ;

// 当前象素的梯度是局部的最大值

// 该点可能是个边界点

if(dTmp>=dTmp1 && dTmp>=dTmp2)

{

pUnchRst[nPos] = 128 ;

}

else

{

// 不可能是边界点

pUnchRst[nPos] = 0 ;

}

}

} //else

} // for

}

} 

/*************************************************************************

 *

 * 函数名称：

 *   TraceEdge()

 *

 * 输入参数:

 *   int    x - 跟踪起点的x坐标 

 *   int    y - 跟踪起点的y坐标

 *   int nLowThd - 判断一个点是否为边界点的低阈值

 *   unsigned char *pUnchEdge - 记录边界点的缓冲区

 *   int *pnMag               - 梯度幅度图

 *   int nWidth               - 图象数据宽度

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   递归调用  

 *   从(x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pUnchEdge中没有处理并且

 *   可能是边界点的那些象素(=128),象素值为0表明该点不可能是边界点，象素值

 *   为255表明该点已经被设置为边界点，不必再考虑

 *   

 *   

 *************************************************************************

 */

void TraceEdge (int y, int x, int nLowThd, unsigned char *pUnchEdge, int *pnMag, int nWidth) 

{ 

// 对8邻域象素进行查询

int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;

int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;

int yy ;

int xx ;

int k  ;

for(k=0; k<8; k++)

{

yy = y + yNb[k] ;

xx = x + xNb[k] ;

// 如果该象素为可能的边界点，又没有处理过

// 并且梯度大于阈值

if(pUnchEdge[yy*nWidth+xx] == 128  && pnMag[yy*nWidth+xx]>=nLowThd)

{

// 把该点设置成为边界点

pUnchEdge[yy*nWidth+xx] = 255 ;

// 以该点为中心进行跟踪

TraceEdge(yy, xx, nLowThd, pUnchEdge, pnMag, nWidth);

}

}

} 

/*************************************************************************

 *

 * 函数名称：

 *   EstimateThreshold()

 *

 * 输入参数:

 *   int *pnMag               - 梯度幅度图

 *  int nWidth               - 图象数据宽度

 *  int nHeight              - 图象数据高度

 *   int *pnThdHigh           - 高阈值

 *   int *pnThdLow            - 低阈值

 *  double dRatioLow         - 低阈值和高阈值之间的比例

 *  double dRatioHigh        - 高阈值占图象象素总数的比例

 *   unsigned char *pUnchEdge - 经过non-maximum处理后的数据

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   经过non-maximum处理后的数据pUnchEdge，统计pnMag的直方图，确定阈值。

 *   本函数中只是统计pUnchEdge中可能为边界点的那些象素。然后利用直方图，

 *   根据dRatioHigh设置高阈值，存储到pnThdHigh。利用dRationLow和高阈值，

 *   设置低阈值，存储到*pnThdLow。dRatioHigh是一种比例：表明梯度小于

 *   *pnThdHigh的象素数目占象素总数目的比例。dRationLow表明*pnThdHigh

 *   和*pnThdLow的比例，这个比例在canny算法的原文里，作者给出了一个区间。

 *

 *************************************************************************

 */

void EstimateThreshold(int *pnMag, int nWidth, int nHeight, int *pnThdHigh,int *pnThdLow, 

 unsigned char * pUnchEdge, double dRatioHigh, double dRationLow) 

{ 

// 循环控制变量

int y;

int x;

int k;

// 该数组的大小和梯度值的范围有关，如果采用本程序的算法，那么梯度的范围不会超过pow(2,10)

int nHist[1024] ;

// 可能的边界数目

int nEdgeNb     ;

// 最大梯度值

int nMaxMag     ;

int nHighCount  ;

nMaxMag = 0     ; 

// 初始化

for(k=0; k<1024; k++) 

{

nHist[k] = 0; 

}

// 统计直方图，然后利用直方图计算阈值

for(y=0; y<nHeight; y++)

{

for(x=0; x<nWidth; x++)

{

// 只是统计那些可能是边界点，并且还没有处理过的象素

if(pUnchEdge[y*nWidth+x]==128)

{

nHist[ pnMag[y*nWidth+x] ]++;

}

}

}

nEdgeNb = nHist[0]  ;

nMaxMag = 0         ;

// 统计经过“非最大值抑止(non-maximum suppression)”后有多少象素

for(k=1; k<1024; k++)

{

if(nHist[k] != 0)

{

// 最大梯度值

nMaxMag = k;

}

// 梯度为0的点是不可能为边界点的

// 经过non-maximum suppression后有多少象素

nEdgeNb += nHist[k];

}

// 梯度比高阈值*pnThdHigh小的象素点总数目

nHighCount = (int)(dRatioHigh * nEdgeNb +0.5);

k = 1;

nEdgeNb = nHist[1];

// 计算高阈值

while( (k<(nMaxMag-1)) && (nEdgeNb < nHighCount) )

{

k++;

nEdgeNb += nHist[k];

}

// 设置高阈值

*pnThdHigh = k ;

// 设置低阈值

*pnThdLow  = (int)((*pnThdHigh) * dRationLow+ 0.5);

}

/*************************************************************************

 *

 * 函数名称：

 *   Hysteresis()

 *

 * 输入参数:

 *   int *pnMag               - 梯度幅度图

 *  int nWidth               - 图象数据宽度

 *  int nHeight              - 图象数据高度

 *  double dRatioLow         - 低阈值和高阈值之间的比例

 *  double dRatioHigh        - 高阈值占图象象素总数的比例

 *   unsigned char *pUnchEdge - 记录边界点的缓冲区

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   本函数实现类似“磁滞现象”的一个功能，也就是，先调用EstimateThreshold

 *   函数对经过non-maximum处理后的数据pUnchSpr估计一个高阈值，然后判断

 *   pUnchSpr中可能的边界象素(=128)的梯度是不是大于高阈值nThdHigh，如果比

 *   该阈值大，该点将作为一个边界的起点，调用TraceEdge函数，把对应该边界

 *   的所有象素找出来。最后，当整个搜索完毕时，如果还有象素没有被标志成

 *   边界点，那么就一定不是边界点。

 *   

 *************************************************************************

 */

void Hysteresis(int *pnMag, int nWidth, int nHeight, double dRatioLow, 

double dRatioHigh, unsigned char *pUnchEdge)

{

// 循环控制变量

int y;

int x;

int nThdHigh ;

int nThdLow  ;

int nPos;

// 估计TraceEdge需要的低阈值，以及Hysteresis函数使用的高阈值

EstimateThreshold(pnMag, nWidth, nHeight, &nThdHigh, 

               &nThdLow, pUnchEdge,dRatioHigh, dRatioLow);

  // 这个循环用来寻找大于nThdHigh的点，这些点被用来当作边界点，然后用

// TraceEdge函数来跟踪该点对应的边界

   for(y=0; y<nHeight; y++)

 {

      for(x=0; x<nWidth; x++)

{

nPos = y*nWidth + x ; 

// 如果该象素是可能的边界点，并且梯度大于高阈值，该象素作为

// 一个边界的起点

if((pUnchEdge[nPos] == 128) && (pnMag[nPos] >= nThdHigh))

{

// 设置该点为边界点

pUnchEdge[nPos] = 255;

TraceEdge(y, x, nThdLow, pUnchEdge, pnMag, nWidth);

}

      }

   }

 // 那些还没有被设置为边界点的象素已经不可能成为边界点

   for(y=0; y<nHeight; y++)

 {

 for(x=0; x<nWidth; x++)

 {

 nPos = y*nWidth + x ; 

 if(pUnchEdge[nPos] != 255)

 {

 // 设置为非边界点

 pUnchEdge[nPos] = 0 ;

 }

 }

 }

}

/*************************************************************************

 *

 * 函数名称：

 *   Canny()

 *

 * 输入参数:

 *   unsigned char *pUnchImage- 图象数据

 *  int nWidth               - 图象数据宽度

 *  int nHeight              - 图象数据高度

 *   double sigma             - 高斯滤波的标准方差

 *  double dRatioLow         - 低阈值和高阈值之间的比例

 *  double dRatioHigh        - 高阈值占图象象素总数的比例

 *   unsigned char *pUnchEdge - canny算子计算后的分割图

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   canny分割算子，计算的结果保存在pUnchEdge中，逻辑1(255)表示该点为

 *   边界点，逻辑0(0)表示该点为非边界点。该函数的参数sigma，dRatioLow

 *   dRatioHigh，是需要指定的。这些参数会影响分割后边界点数目的多少

 *************************************************************************

 */

void Canny(unsigned char *pUnchImage, int nWidth, int nHeight, double sigma,

 double dRatioLow, double dRatioHigh, unsigned char *pUnchEdge)

{

// 经过高斯滤波后的图象数据

unsigned char * pUnchSmooth ;

  

// 指向x方向导数的指针

int * pnGradX ; 

// 指向y方向导数的指针

int * pnGradY ;

// 梯度的幅度

int * pnGradMag ;

pUnchSmooth  = new unsigned char[nWidth*nHeight] ;

pnGradX      = new int [nWidth*nHeight]          ;

pnGradY      = new int [nWidth*nHeight]          ;

pnGradMag    = new int [nWidth*nHeight]          ;

// 对原图象进行滤波

GaussianSmooth(pUnchImage, nWidth, nHeight, sigma, pUnchSmooth) ;

// 计算方向导数

DirGrad(pUnchSmooth, nWidth, nHeight, pnGradX, pnGradY) ;

// 计算梯度的幅度

GradMagnitude(pnGradX, pnGradY, nWidth, nHeight, pnGradMag) ;

// 应用non-maximum 抑制

NonmaxSuppress(pnGradMag, pnGradX, pnGradY, nWidth, nHeight, pUnchEdge) ;

// 应用Hysteresis，找到所有的边界

Hysteresis(pnGradMag, nWidth, nHeight, dRatioLow, dRatioHigh, pUnchEdge);

// 释放内存

delete pnGradX      ;

pnGradX      = NULL ;

delete pnGradY      ;

pnGradY      = NULL ;

delete pnGradMag    ;

pnGradMag    = NULL ;

delete pUnchSmooth ;

pUnchSmooth  = NULL ;

}

void Bianli (unsigned char *pUnchEdge,int x,int y,int nWidth, int nHeight,unsigned short *L,int n) 

{

//往左遍历

if((x-1)>=0)

{

          if(pUnchEdge[y*nWidth+x-1]==255)

  {

             L[y*nWidth+x-1]=n;

 pUnchEdge[y*nWidth+x-1]=0;

             Bianli(pUnchEdge,x-1,y,nWidth,nHeight,L,n);

  }

}

    //往右遍历

if((x+1)<=nWidth)

{

          if(pUnchEdge[y*nWidth+x+1]==255)

  {

             L[y*nWidth+x+1]=n;

 pUnchEdge[y*nWidth+x+1]=0;

             Bianli(pUnchEdge,x+1,y,nWidth,nHeight,L,n);

  }

}

    //往上遍历

if((y-1)>=0)

{

          if(pUnchEdge[(y-1)*nWidth+x]==255)

  {

             L[(y-1)*nWidth+x]=n;

 pUnchEdge[(y-1)*nWidth+x]=0;

             Bianli(pUnchEdge,x,y-1,nWidth,nHeight,L,n);

  }

}

    //往下遍历

if((y+1)>=0)

{

          if(pUnchEdge[(y+1)*nWidth+x]==255)

  {

             L[(y+1)*nWidth+x]=n;

 pUnchEdge[(y+1)*nWidth+x]=0;

             Bianli(pUnchEdge,x,y+1,nWidth,nHeight,L,n);

  }

}

}

double*  BlockDct(double *lpBlockBits, int BlockWidth, int BlockHeight)

{

// 指向源图像的指针

//unsigned char* lpSrc;

// 循环变量

LONG i;

LONG j;

// 进行付立叶变换的宽度和高度（2的整数次方）

LONG w;

LONG h;

//高度和宽度的幂次（2）

int wp;

int hp;

// 图像每行的字节数

//LONG lLineBytes;

// 计算图像每行的字节数

//lLineBytes = WIDTHBYTES(BlockWidth * 8);

// 赋初值

w = 1;

h = 1;

wp = 0;

hp = 0;

// 计算进行离散余弦变换的宽度和高度（2的整数次方）

while(w * 2 <= BlockWidth)

{

w *= 2;

wp++;

}

while(h * 2 <= BlockHeight)

{

h *= 2;

hp++;

}

//分配内存

double *lpdctblock = new double[w*h];

for(i = 0; i < h; i++)

{

// 对y方向进行离散余弦变换

DCT(&lpBlockBits[w * i], &lpdctblock[w * i], wp);

}

// 保存计算结果

for(i = 0; i < h; i++)

{

for(j = 0; j < w; j++)

{

lpBlockBits[j * h + i] = lpdctblock[j + w * i];

}

}

for(j = 0; j < w; j++)

{

// 对x方向进行离散余弦变换

DCT(&lpBlockBits[j * h], &lpdctblock[j * h], hp);

}

/*

// 行

for(i = 0; i < h; i++)

{

// 列

for(j = 0; j < w; j++)

{

// 计算频谱

dTemp = fabs(lpdctblock[j*h+i]);

// 判断是否超过255

if (dTemp > 255)

{

// 对于超过的，直接设置为255

dTemp = 255;

}

// 指向DIB第y行，第x个象素的指针

//lpSrc = (unsigned char*)lpBlockBits + lLineBytes * (BlockHeight - 1 - i) + j;

// 更新源图像

//* (lpSrc) = (BYTE)(dTemp);

lpdctblock[i * w + h] = dTemp ;

}

}

*/

// 返回

return lpdctblock;

}

/*************************************************************************

 *

 * 函数名称：

 *   DIBDct()

 *

 * 参数:

 *   CDib  *pDib       - 指向CDib类的指针

 *

 * 返回值:

 *   BOOL               - 成功返回TRUE，否则返回FALSE。

 *

 * 说明:

 *   该函数用来对图像进行二维离散余弦变换。

 *

 ************************************************************************/

BYTE* DIBDct(BYTE* lpdib ,int lWidth,int lHeight)

{

// 指向源图像的指针

BYTE * lpSrc;

//图象的宽度和高度

// LONG    lWidth;

// LONG    lHeight;

// 循环变量

LONG i;

LONG j;

// 离散余弦变换的宽度和高度，必须为2的整数次方

int lW = 1;

int lH = 1;

int wp = 0;

int hp = 0;

// 中间变量

double dTemp;

// 保证离散余弦变换的宽度和高度为2的整数次方

while(lW * 2 <= lWidth)

{

lW = lW * 2;

wp++;

}

while(lH * 2 <= lHeight)

{

lH = lH * 2;

hp++;

}

// 分配内存

double *dpf = new double[lW * lH];

double *dpF = new double[lW * lH];

// 时域赋值

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

// 指针指向位图i行j列的象素

// lpSrc = lpdib + lWidth * (i) + j;

lpSrc = lpdib + lWidth * (lHeight - 1 - i) + j;

// 将象素值赋给时域数组

dpf[j + i * lW] = *(lpSrc);

}

}

for(i = 0; i < lH; i++)

// y方向进行离散余弦变换

DCT(&dpf[lW * i], &dpF[lW * i], wp);

// 保存计算结果

for(i = 0; i < lH; i++)

for(j = 0; j < lW; j++)

dpf[j * lH + i] = dpF[j + lW * i];

for(j = 0; j < lW; j++)

// x方向进行离散余弦变换

DCT(&dpf[j * lH], &dpF[j * lH], hp);

// 频谱的计算

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

dTemp = fabs(dpF[i*lW + j]);

// 是否超过255

if (dTemp > 255)

// 如果超过，设置为255

dTemp = 255;

// 指针指向位图i行j列的象素

//lpSrc= lpdib + lWidth * i + j;

lpSrc = lpdib + lWidth * (lHeight - 1 - i) + j;

// 更新源图像

* (lpSrc) = (BYTE)(dTemp);

}

}

// 释放内存

delete dpf;

delete dpF;

// 返回

return lpSrc;

}

/*************************************************************************

 *

 * 函数名称：

 *   DCT()

 *

 * 参数:

 *   double * dpf - 指向时域值的指针

 *   double * dpF - 指向频域值的指针

 *   r     - 2的幂数

 *

 * 返回值:

 *   无。

 *

 * 说明:

 *   实现一维快速离散余弦变换。

 *

 ************************************************************************/

VOID WINAPI DCT(double *dpf, double *dpF, int r)

{

double PI = 3.1415926; 

// 离散余弦变换点数

LONG lNum;

// 循环变量

int i;

// 中间变量

double dTemp;

complex<double> *comX;

// 离散余弦变换点数

lNum = 1<<r;

// 分配内存

comX = new complex<double>[lNum*2];

// 赋初值0

memset(comX, 0, sizeof(complex<double>) * lNum * 2);

// 将时域点转化为复数形式

for(i=0;i<lNum;i++)

{

comX[i] = complex<double> (dpf[i], 0);

}

// 调用快速付立叶变换

FFT_1D(comX,comX,r+1);

// 调整系数

dTemp = 1/sqrt(lNum);

// 求dpF[0]

dpF[0] = comX[0].real() * dTemp;

dTemp *= sqrt(2);

// 求dpF[u]

for(i = 1; i < lNum; i++)

{

dpF[i]=(comX[i].real() * cos(i*PI/(lNum*2)) 

+ comX[i].imag() * sin(i*PI/(lNum*2))) * dTemp;

}

// 释放内存

    delete comX;

}

/*************************************************************************

 *

 * 函数名称：

 *   IDCT()

 *

 * 参数:

 *   double * dpF - 指向频域值的指针

 *   double * dpf - 指向时域值的指针

 *   r －2的幂数

 *

 * 返回值:

 *   无。

 *

 * 说明:

 *   实现一维快速离散余弦反变换。

 *

 ************************************************************************/

/*

VOID WINAPI IDCT(double *dpF, double *dpf, int r)

{

double PI = 3.1415926; 

// 离散余弦反变换点数

LONG lNum;

// 计算离散余弦变换点数

lNum = 1<<r;

// 循环变量

int i;

// 中间变量

double dTemp, d0;

complex<double> *comX;

    // 给复数变量分配内存

comX = new complex<double>[lNum*2];

// 赋初值0

memset(comX, 0, sizeof(complex<double>) * lNum * 2);

// 将频域变换后点写入数组comX

for(i=0;i<lNum;i++)

{

comX[i] = complex<double> (cos(i*PI/(lNum*2)) * dpF[i], sin(i*PI/(lNum*2) * dpF[i]));

}

// 调用快速付立叶反变换

IFFT_1D(comX,comX,r+1);

// 调整系数

dTemp = sqrt(2.0/lNum);

d0 = (sqrt(1.0/lNum) - dTemp) * dpF[0];

// 计算dpf(x)

for(i = 0; i < lNum; i++)

{

dpf[i] = d0 + 2 * lNum * comX[i].real()* dTemp ;

}

// 释放内存

delete comX;

}

*/

/*************************************************************************

 *

 * 函数名称：

 *   FFT_1D()

 *

 * 输入参数:

 *   complex<double> * pCTData - 指向时域数据的指针，输入的需要变换的数据

 *   complex<double> * pCFData - 指向频域数据的指针，输出的经过变换的数据

 *   nLevel －傅立叶变换蝶形算法的级数，2的幂数，

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   一维快速傅立叶变换。

 *

 **************************************************************************/

VOID WINAPI FFT_1D(complex<double> * pCTData, complex<double> * pCFData, int nLevel)

{

// 循环控制变量

int i;

int     j;

int     k;

double PI = 3.1415926; 

// 傅立叶变换点数

int nCount =0 ;

// 计算傅立叶变换点数

nCount =(int)pow(2,nLevel) ;

// 某一级的长度

int nBtFlyLen;

nBtFlyLen = 0 ;

// 变换系数的角度 ＝2 * PI * i / nCount

double dAngle;

complex<double> *pCW ;

// 分配内存，存储傅立叶变化需要的系数表

pCW  = new complex<double>[nCount / 2];

    // 计算傅立叶变换的系数

for(i = 0; i < nCount / 2; i++)

{

dAngle = -2 * PI * i / nCount;

pCW[i] = complex<double> ( cos(dAngle), sin(dAngle) );

}

// 变换需要的工作空间

complex<double> *pCWork1,*pCWork2; 

// 分配工作空间

pCWork1 = new complex<double>[nCount];

pCWork2 = new complex<double>[nCount];

// 临时变量

complex<double> *pCTmp;

// 初始化，写入数据

memcpy(pCWork1, pCTData, sizeof(complex<double>) * nCount);

// 临时变量

int nInter; 

nInter = 0;

// 蝶形算法进行快速傅立叶变换

for(k = 0; k < nLevel; k++)

{

for(j = 0; j < (int)pow(2,k); j++)

{

//计算长度

nBtFlyLen = (int)pow( 2,(nLevel-k) );

//倒序重排，加权计算

for(i = 0; i < nBtFlyLen/2; i++)

{

nInter = j * nBtFlyLen;

pCWork2[i + nInter] = 

pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2];

pCWork2[i + nInter + nBtFlyLen / 2] =

(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2]) 

* pCW[(int)(i * pow(2,k))];

}

}

// 交换 pCWork1和pCWork2的数据

pCTmp   = pCWork1 ;

pCWork1 = pCWork2 ;

pCWork2 = pCTmp ;

}

// 重新排序

for(j = 0; j < nCount; j++)

{

nInter = 0;

for(i = 0; i < nLevel; i++)

{

if ( j&(1<<i) )

{

nInter += 1<<(nLevel-i-1);

}

}

pCFData[j]=pCWork1[nInter];

}

// 释放内存空间

delete pCW;

delete pCWork1;

delete pCWork2;

pCW = NULL ;

pCWork1 = NULL ;

pCWork2 = NULL ;

}

/*************************************************************************

 *

 * 函数名称：

 *    IFFT_1D()

 *

 * 输入参数:

 *   complex<double> * pCTData - 指向时域数据的指针，输入的需要反变换的数据

 *   complex<double> * pCFData - 指向频域数据的指针，输出的经过反变换的数据

 *   nLevel －傅立叶变换蝶形算法的级数，2的幂数，

 *

 * 返回值:

 *   无

 *

 * 说明:

 *   一维快速傅立叶反变换。

 *

 *************************************************************************/

/*

VOID WINAPI IFFT_1D(complex<double> * pCFData, complex<double> * pCTData, int nLevel)

{

// 循环控制变量

int i;

// 傅立叶反变换点数

int nCount;

// 计算傅立叶变换点数

nCount = (int)pow(2,nLevel) ;

// 变换需要的工作空间

complex<double> *pCWork;

// 分配工作空间

pCWork = new complex<double>[nCount];

// 将需要反变换的数据写入工作空间pCWork

memcpy(pCWork, pCFData, sizeof(complex<double>) * nCount);

// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭

// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭

for(i = 0; i < nCount; i++)

{

pCWork[i] = complex<double> (pCWork[i].real(), -pCWork[i].imag());

}

// 调用快速傅立叶变换实现反变换，结果存储在pCTData中

FFT_1D(pCWork, pCTData, nLevel);

// 求时域点的共轭，求得最终结果

// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据

// 相差一个系数nCount

for(i = 0; i < nCount; i++)

{

pCTData[i] = 

complex<double> (pCTData[i].real() / nCount, -pCTData[i].imag() / nCount);

}

// 释放内存

delete pCWork;

pCWork = NULL;

}

*/

/*************************************************************************

 *

 * 函数名称：

 *   DIBDct()

 *

 * 参数:

 *   LPBYTE  lpDIBBits - 指向CDib类的指针

 *   int nWidth        - 待变换的图像的宽度

 *   int nHeight       - 待变换的图像的高度 

 *

 * 返回值:

 *   BOOL               - 成功返回TRUE，否则返回FALSE。

 *

 * 说明:

 *   该函数用来对图像进行二维离散余弦变换。

 *

 ************************************************************************/

 /*

BYTE* DIB8Dct(LPBYTE  lpDIBBits ,int nWidth,int nHeight)

{

// 指向源图像的指针

BYTE* lpSrc;

// 循环变量

LONG i;

LONG j;

LONG    w;

LONG    h;

// 离散余弦变换的宽度和高度，必须为2的整数次方

LONG lW = 1;

LONG lH = 1;

int wp = 0;

int hp = 0;

//DCT块大小

const int N=8;

// 中间变量

double dTemp;

// 保证离散余弦变换的宽度和高度为2的整数次方

while(lW * 2 <= nWidth)

{

lW = lW * 2;

wp++;

}

while(lH * 2 <= nHeight)

{

lH = lH * 2;

hp++;

}

// 分配内存

double *dpf = new double[lW * lH];

double *dpF = new double[lW * lH];

double *block = new double[N * N];

double *lpdctblock = new double[N * N]; 

// 时域赋值

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

// 指针指向位图i行j列的象素

// lpSrc = lpDIBBits + nWidth * ( i) + j;

lpSrc = lpDIBBits + nWidth * (nHeight - 1 - i) + j;

// 将象素值赋给时域数组

dpf[j + i * lW] = *(lpSrc);

}

}

   int nw=lW/N;

   int nh=lH/N;

   for (i = 0; i < nh; i++)

   {

   for(j = 0; j < nw; j++)

   {

for (h = 0; h< N; h++)

{

for(w = 0; w< N; w++)

{

block[h * N + w] = dpf[(N * i + h) * nWidth + N * j + w];

}

}

//对指定的块进行DCT变换

lpdctblock = BlockDct(block,N ,N );

//把变换的结果再赋给图像矩阵

for (h = 0; h< N; h++)

{

for(w = 0; w< N; w++)

{

dpF[(N * i + h) * nWidth + N * j + w] = lpdctblock[h * N + w];

}

}

   }

   }

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

// 计算频谱

dTemp = fabs(dpF[i*lW + j]);

// 判断是否超过255

if (dTemp > 255)

{

// 对于超过的，直接设置为255

dTemp = 255;

}

lpSrc = lpDIBBits + nWidth * (nHeight - 1 - i) + j;

// 更新源图像

* (lpSrc) = (BYTE)(dTemp);

}

}

// 释放内存

delete dpf;

delete dpF;

delete block;

// delete lpdctblock;

// delete lpdctblock;

// delete lpDIBBits;

 

// 返回

return lpSrc;

}

*/

/*************************************************************************

 *

 * 函数名称：

 *   DIBfitDct()

 *

 * 参数:

 *   LPBYTE  lpImage   - 指向待处理图像的指针

 *   int nWidth        - 待变换的图像的宽度

 *   int nHeight       - 待变换的图像的高度 

 *   int *Threshold    - 图像分割的阈值 

 *

 * 返回值:

 *   BYTE*             - 自适应图像DCT压缩的数据

 *

 * 说明:

 *   该函数用来对图像根据平滑边缘区自适应进行二维离散余弦变换。

 *

 ************************************************************************/

double * DIB16Dct(double * sourcedata)

{

double s08,s412,s210,s614,s19,s513,s311,s715,d08,d412,d210,d614,d19,d513,d311,d715;

double * data;

double dtemp;

int i;

dtemp = sqrt(2)/4;

   //行DCT变换

data=sourcedata;

for(i=0; i<16; i++)

{

s08 = data[0] + data[8];

s412 = data[4] + data[12];

s210 = data[2] + data[10]; 

s614 = data[6] + data[14];

s19 = data[1] + data[9];

s513 = data[5] + data[13];

s311 = data[3] + data[11];

s715 = data[7] + data[15];

d08 = data[0] - data[8];

d412 = data[4] - data[12];

d210 = data[2] - data[10]; 

d614 = data[6] - data[14];

d19 = data[1] - data[9];

d513 = data[5] - data[13];

d311 = data[3] - data[11];

d715 = data[7] - data[15];

data[0] =0.25* (s08+s412+s210+s614+s19+s311+s715);

data[1] = dtemp*((-0.4716)*d513+0.2904*d311-0.09806*d412+0.9952*d08+0.6344*d210+0.8820*d19-0.9572*d715-0.7732*d614);

data[2] = dtemp*(-0.8316*s311+0.8316*s715-0.1951*s210+0.1951*s614+0.5556*s19-0.5556*s513+0.9808*s08-0.9808*s412);

data[3] = dtemp*(-0.7732*d311+0.9954*d513+0.2904*d412+0.9570*d08-0.8820*d210+0.09800*d19-0.6348*d715+0.4716*d614);

data[4] = dtemp*((-0.3828)*(s19+s513-s311-s715)+0.9240*(s08+s412-s210-s614));

data[5] = dtemp*(-0.6344*d513+0.9954*d311-0.4716*d412+0.8820*d08-0.2904*d210-0.7732*d19-0.09814*d715+0.9572*d614);

data[6] = dtemp*(0.1951*s311-0.1951*s715-0.9810*s19+0.9810*s513+0.5556*s210-0.5556*s614+0.8316*s08-0.8316*s412); 

    data[7] = dtemp*(-0.8824*d311-0.2904*d513+0.6344*d412+0.7730*d08+0.9954*d210-0.9572*d19+0.4716*d715-0.09796*d614);

data[8] = dtemp*0.7071*(s08+s412+s210+s614-s19-s513-s311-s715);

data[9] = dtemp*(0.9572*d513+0.4716*d311-0.7732*d412+0.6344*d08-0.09806*d210-0.2904*d19+0.8820*d715-0.9954*d614);

data[10] = dtemp*(0.9810*s311-0.9810*s715-0.8316*s210+0.8316*s614+0.1951*s19-0.1951*s513+0.5556*s08-0.5556*s412);

data[11] = dtemp*(0.09796*d311-0.7732*d513+0.8820*d412+0.4714*d08-0.9572*d210+0.6344*d19+0.9956*d715-0.2904*d614);

data[12] = dtemp*((-0.9238)*(-s19-s513+s311+s715)+0.3827*(s08+s412-s210-s614));

    data[13] = dtemp*(-0.09814*d513-0.6344*d311-0.9572*d412+0.2904*d08+0.4716*d210+0.9954*d19+0.7734*d715+0.8820*d614);

data[14] = dtemp*(0.5558*s311-0.5558*s715+0.8316*s19-0.8316*s513+0.9810*s210-0.9810*s614+0.1951*s08-0.1951*s412);

data[15] = dtemp*(0.9572*d311+0.8824*d513+0.9954*d412+0.09800*d08+0.7732*d210+0.4716*d19+0.2906*d715+0.6348*d614);

data += 16;

}

//列DCT变换

    data = sourcedata;

    for(i=0; i<16; i++)

{

s08 = data[0*16] + data[8*16];

s412 = data[4*16] + data[12*16];

s210 = data[2*16] + data[10*16]; 

s614 = data[6*16] + data[14*16];

s19 = data[1*16] + data[9*16];

s513 = data[5*16] + data[13*16];

s311 = data[3*16] + data[11*16];

s715 = data[7*16] + data[15*16];

d08 = data[0*16] -data[8*16];

d412 = data[4*16] - data[12*16];

d210 = data[2*16] - data[10*16]; 

d614 = data[6*16] - data[14*16];

d19 = data[1*16] - data[9*16];

d513 = data[5*16] - data[13*16];

d311 = data[3*16] - data[11*16];

    d715 = data[7*16] - data[15*16];

    data[0*16] =0.25* (s08+s412+s210+s614+s19+s311+s715);

data[1*16] = dtemp*((-0.4716)*d513+0.2904*d311-0.09806*d412+0.9952*d08+0.6344*d210+0.8820*d19-0.9572*d715-0.7732*d614);

data[2*16] = dtemp*(-0.8316*s311+0.8316*s715-0.1951*s210+0.1951*s614+0.5556*s19-0.5556*s513+0.9808*s08-0.9808*s412);

data[3*16] = dtemp*(-0.7732*d311+0.9954*d513+0.2904*d412+0.9570*d08-0.8820*d210+0.09800*d19-0.6348*d715+0.4716*d614);

data[4*16] = dtemp*((-0.3828)*(s19+s513-s311-s715)+0.9240*(s08+s412-s210-s614));

data[5*16] = dtemp*(-0.6344*d513+0.9954*d311-0.4716*d412+0.8820*d08-0.2904*d210-0.7732*d19-0.09814*d715+0.9572*d614);

data[6*16] = dtemp*(0.1951*s311-0.1951*s715-0.9810*s19+0.9810*s513+0.5556*s210-0.5556*s614+0.8316*s08-0.8316*s412); 

    data[7*16] = dtemp*(-0.8824*d311-0.2904*d513+0.6344*d412+0.7730*d08+0.9954*d210-0.9572*d19+0.4716*d715-0.09796*d614);

data[8*16] = dtemp*0.7071*(s08+s412+s210+s614-s19-s513-s311-s715);

data[9*16] = dtemp*(0.9572*d513+0.4716*d311-0.7732*d412+0.6344*d08-0.09806*d210-0.2904*d19+0.8820*d715-0.9954*d614);

data[10*16] = dtemp*(0.9810*s311-0.9810*s715-0.8316*s210+0.8316*s614+0.1951*s19-0.1951*s513+0.5556*s08-0.5556*s412);

data[11*16] = dtemp*(0.09796*d311-0.7732*d513+0.8820*d412+0.4714*d08-0.9572*d210+0.6344*d19+0.9956*d715-0.2904*d614);

data[12*16] = dtemp*((-0.9238)*(-s19-s513+s311+s715)+0.3827*(s08+s412-s210-s614));

    data[13*16] = dtemp*(-0.09814*d513-0.6344*d311-0.9572*d412+0.2904*d08+0.4716*d210+0.9954*d19+0.7734*d715+0.8820*d614);

data[14*16] = dtemp*(0.5558*s311-0.5558*s715+0.8316*s19-0.8316*s513+0.9810*s210-0.9810*s614+0.1951*s08-0.1951*s412);

data[15*16] = dtemp*(0.9572*d311+0.8824*d513+0.9954*d412+0.09800*d08+0.7732*d210+0.4716*d19+0.2906*d715+0.6348*d614);

data ++;

}

return sourcedata;

}

double * DIB4Dct(double * sourcedata)

{

double const C1=0.9239;

double const C2=0.7071;

double const C3=0.3827;

double S02,S13,D02,D13,SS0213,SD0213,DSO213;

double * data;

double dtemp;

int i;

data=sourcedata;

//  后面有专门的地方使亮度值减去128

//  行DCT变换

for(i=0; i<4; i++)

{

S02 = data[0] + data[2];

S13 = data[1] + data[3];

D02 = data[0] - data[2];

D13 = data[1] - data[3];

SS0213 = S02 + S13;

SD0213 = D02 + D13;

DSO213 = S02 - S13;

dtemp = sqrt(2)/2; 

data[0] = 0.5 * SS0213;

data[1] = dtemp * C1 * SD0213;

data[2] = dtemp * C2 * DSO213;

data[3] = dtemp * C3 * SD0213;

data += 4;

}

/////////////////////////////////////////////////////

//  列DCT变换

    data=sourcedata;

for(i=0; i<4; i++)

{

S02 = data[0*4] + data[2*4];

S13 = data[1*4] + data[3*4];

D02 = data[0*4] - data[2*4];

D13 = data[1*4] - data[3*4];

SS0213 = S02 + S13;

SD0213 = D02 + D13;

DSO213 = S02 - S13;

dtemp = sqrt(2)/2; 

data[0*4] = 0.5 * SS0213;

data[1*4] = dtemp * C1 * SD0213;

data[2*4] = dtemp * C2 * DSO213;

data[3*4] = dtemp * C3 * SD0213;

data ++;

}

   return sourcedata;

}

double * DIB8Dct(double * sourcedata)

{

double const C1=0.9808;

double const C2=0.9239;

double const C3=0.8315;

double const C4=0.7071;

double const C5=0.5556;

double const C6=0.3827;

double const C7=0.1951;

double S18,S27,S36,S45,S1845,S2736;

double D18,D27,D36,D45,D1845,D2736;

double * data;

int i;

/* for(i=0;i<64;i++)

{

sourcedata[i]=sourcedata[i]-128;//图象数据减去128

}

*/

//  后面有专门的地方使亮度值减去128

//  行DCT变换

data=sourcedata;

   for(i=0;i<8;i++)

   {

    S18=data[0]+data[7];

    S27=data[1]+data[6];

    S36=data[2]+data[5];

    S45=data[3]+data[4];

    S1845=S18+S45;

    S2736=S27+S36;

    D18=data[0]-data[7]; 

    D27=data[1]-data[6];

    D36=data[2]-data[5];

    D45=data[3]-data[4];

    D1845=S18-S45;

    D2736=S27-S36;

    data[0]=0.5*(C4*(S1845+S2736));

    data[1]=0.5*(C1*D18+C3*D27+C5*D36+C7*D45);

    data[2]=0.5*(C2*D1845+C6*D2736);

    data[3]=0.5*(C3*D18-C7*D27-C1*D36-C5*D45);

    data[4]=0.5*(C4*(S1845-S2736));

    data[5]=0.5*(C5*D18-C1*D27+C7*D36+C3*D45);

    data[6]=0.5*(C6*D1845-C2*D2736);

    data[7]=0.5*(C7*D18-C5*D27+C3*D36-C1*D45);

data+=8;

   }

/////////////////////////////////////////////////////

//  列DCT变换

   data=sourcedata;

   for(i=0;i<8;i++)

   {

    S18=data[0*8]+data[7*8];

    S27=data[1*8]+data[6*8];

    S36=data[2*8]+data[5*8];

    S45=data[3*8]+data[4*8];

    S1845=S18+S45;

    S2736=S27+S36;

    D18=data[0*8]-data[7*8]; 

    D27=data[1*8]-data[6*8];

    D36=data[2*8]-data[5*8];

    D45=data[3*8]-data[4*8];

    D1845=S18-S45;

    D2736=S27-S36;

    data[0*8]=0.5*(C4*(S1845+S2736));

    data[1*8]=0.5*(C1*D18+C3*D27+C5*D36+C7*D45);

    data[2*8]=0.5*(C2*D1845+C6*D2736);

    data[3*8]=0.5*(C3*D18-C7*D27-C1*D36-C5*D45);

    data[4*8]=0.5*(C4*(S1845-S2736));

    data[5*8]=0.5*(C5*D18-C1*D27+C7*D36+C3*D45);

    data[6*8]=0.5*(C6*D1845-C2*D2736);

    data[7*8]=0.5*(C7*D18-C5*D27+C3*D36-C1*D45);

data++;

   }

   return sourcedata;

}

BYTE* DIBnDct(LPBYTE  lpDIBBits ,int nWidth,int nHeight,int N)

{

// 指向源图像的指针

BYTE* lpSrc;

// 循环变量

LONG i;

LONG j;

LONG    w;

LONG    h;

// LONG    k;

// 离散余弦变换的宽度和高度，必须为2的整数次方

LONG lW = 1;

LONG lH = 1;

int wp = 0;

int hp = 0;

//DCT块大小

//const int N=8;

// 中间变量

double dTemp;

//量化表

// 保证离散余弦变换的宽度和高度为2的整数次方

while(lW * 2 <= nWidth)

{

lW = lW * 2;

wp++;

}

while(lH * 2 <= nHeight)

{

lH = lH * 2;

hp++;

}

// 分配内存

double *dpf = new double[lW * lH];

double *dpF = new double[lW * lH];

double *block = new double[N * N];

double *lpdctblock = new double[N * N]; 

   

// 时域赋值(t=16ms)

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

// 指针指向位图i行j列的象素

// lpSrc = lpDIBBits + nWidth * ( i) + j;

// lpSrc = lpDIBBits + nWidth * (nHeight - 1 - i) + j;

// 将象素值赋给时域数组

// dpf[j + i * lW] = *(lpSrc);

dpf[j + i * lW] = lpDIBBits[nWidth * (nHeight - 1 - i) + j];

}

}

   int nw=lW/N;

   int nh=lH/N;

   for (i = 0; i < nh; i++)

   {

    

   for(j = 0; j < nw; j++)

   {

for (h = 0; h< N; h++)

{

for(w = 0; w< N; w++)

{

block[h * N + w] = dpf[(N * i + h) * nWidth + N * j + w];

}

}

//对指定的块进行DCT变换

//lpdctblock = BlockDct(block,N ,N );

//lpdctblock = Pre_DCT(block);

 if(N==4)

 lpdctblock = DIB4Dct(block);

 else if(N==8)

lpdctblock = DIB8Dct(block);

 else if (N==16)

lpdctblock = DIB16Dct(block);

 else if (N==32)

lpdctblock = DIB32Dct(block);

 else

lpdctblock = BlockDct(block,N ,N );

 

 

//把变换的结果再赋给图像矩阵

for (h = 0; h< N; h++)

{

for(w = 0; w< N; w++)

{

dpF[(N * i + h) * nWidth + N * j + w] = lpdctblock[h * N + w];

}

}

 

   }

 

   }

for(i = 0; i < lH; i++)

{

for(j = 0; j < lW; j++)

{

// 计算频谱

dTemp = fabs(dpF[i*lW + j]);

// 判断是否超过255

if (dTemp > 255)

{

// 对于超过的，直接设置为255

dTemp = 255;

}

//lpSrc = lpDIBBits + nWidth * i + j;

lpSrc = lpDIBBits + nWidth * (nHeight - 1 - i) + j;

// 更新源图像

* (lpSrc) = (BYTE)(dTemp);

// lpDIBBits[nWidth * (nHeight - 1 - i) + j] = (BYTE)(dTemp)；

}

}

// 释放内存

delete dpf;

delete dpF;

delete block;

// delete lpdctblock;

// delete lpdctblock;

// delete lpDIBBits;

 

// 返回

return lpSrc;

}

double * DIB32Dct(double * sourcedata)

{

double  s016, s824, s420, s1228, s218, s1026, s622, s1430 ;

double  s117, s925, s521, s1329, s319, s1127, s723, s1531;

double  d016, d824, d420, d1228, d218, d1026, d622, d1430 ;

double  d117, d925, d521, d1329, d319, d1127, d723, d1531;

double * data;

int i;

   //行DCT变换

data=sourcedata;

for(i=0; i<32; i++)

{

s016=data[0]+data[31];

s824=data[16]+data[15];

s420=data[8]+data[23];

s1228=data[24]+data[7];

s218=data[4]+data[27];

    s1026=data[20]+data[11];

s622=data[12]+data[19];

    s1430=data[28]+data[3];

s117=data[2]+data[29];

s925=data[18]+data[13];

s521=data[10]+data[21];

s1329=data[26]+data[5];

s319=data[6]+data[25];

s1127=data[22]+data[9];

s723=data[14]+data[17];

s1531=data[30]+data[1];

    

d016=data[0]-data[31];

d824=data[16]-data[15];

d420=data[8]-data[23];

d1228=data[24]-data[7];

d218=data[4]-data[27];

    d1026=data[20]-data[11];

d622=data[12]-data[19];

    d1430=data[28]-data[3];

d117=data[2]-data[29];

d925=data[18]-data[13];

d521=data[10]-data[21];

d1329=data[26]-data[5];

d319=data[6]-data[25];

d1127=data[22]-data[9];

d723=data[14]-data[17];

d1531=data[30]-data[1];

    data[0]=.1768*(s1430+s1026+s824+s117+s925+s521+s1228+s016+s420+s1329+s319+s1127+s218+s723+s1531+s622);

data[1]=.8422e-1*d622-.2473*d1531+.2425*d117+.2008*d319+.3668e-1*d723+.1679*d420-.1227e-1*d824-.2144*d1329+.1285*d521-.1489*d1127-.1069*d1026-.2354*d1430+.2497*d016+.2260*d218-.6074e-1*d925-.1852*d1228;

data[2]=0.25*((.8819)*(s117-s925)+(-.7730)*(s622-s1430)+(-.4714)*(s521-s1329)+(.6344)*(s218-s1026)+(-.9801e-1)*(s420-s1228)+(-.9569)*(s723-s1531)+(.2903)*(s319-s1127)+(.9952)*(s016-s824));

data[3]=-.2144*d622-.2260*d1531+.1852*d117-.8422e-1*d319-.1069*d723-.2008*d420+.3668e-1*d824-.2497*d521+.1227e-1*d1329+.2354*d1127+.2425*d1026-.1285*d1430+.2473*d016+.1679*d925+.6075e-1*d218+.1489*d1228;

data[4]=0.25*((-.8315)*(s319+s1127-s723-s1531)+(-.1951)*(s218+s1026-s622-s1430)+(.5556)*(s117+s925-s521-s1329)+(.9808)*(s016+s824-s420-s1228));

data[5]=.2473*d622-.1852*d1531-.2497*d319+.8422e-1*d117+.1679*d723-.1285*d420-.6074e-1*d824+.2260*d1329+.1069*d521-.2008*d1026+.1227e-1*d1127+.3668e-1*d1430+.2425*d016-.1489*d218-.2354*d925+.2144*d1228;

data[6]=0.25*((.9802e-1)*(s117-s925)+(-.4714)*(-s622+s1430)+(-.8819)*(s218-s1026)+(-.9952)*(-s521+s1329)+(-.7730)*(s319-s1127)+(.6344)*(-s723+s1531)+(.9569)*(s016-s824)+(-.2903)*(-s420+s1228));

data[7]=-.1679*d622-.1285*d1531-.6075e-1*d319-.3668e-1*d117-.2144*d723+.2260*d420+.8422e-1*d824+.2008*d1329+.1489*d521-.2425*d1127+.1227e-1*d1026+.1852*d1430+.2354*d016+.2473*d925-.2497*d218-.1069*d1228;

data[8]=0.25*((-.3827)*(s117+s925+s521+s1329-s319-s1127-s723-s1531)+(.9239)*(s016+s824+s420+s1228-s218-s1026-s622-s1430));

data[9]=.1227e-1*d622-.6074e-1*d1531-.1489*d117+.2144*d319+.2425*d723+.8422e-1*d420-.1069*d824-.2473*d521-.3668e-1*d1329+.1285*d1127+.1852*d1026+.2497*d1430+.2260*d016-.2008*d925-.1679*d218-.2354*d1228;

data[10]=0.25*((-.7730)*(s117-s925)+(-.6344)*(s521-s1329)+(-.9569)*(-s622+s1430)+(-.4714)*(s420-s1228)+(-.9952)*(-s319+s1127)+(.2903)*(-s218+s1026)+(.9801e-1)*(-s723+s1531)+(.8819)*(s016-s824));

data[11]=.1489*d622+.1227e-1*d1531+.1852*d319-.2260*d117-.2497*d723-.2425*d420+.1285*d824-.2354*d1329+.8422e-1*d521+.1679*d1127-.2473*d1026+.2008*d1430+.2144*d016+.1069*d925+.3668e-1*d218+.6075e-1*d1228;

data[12]=0.25*((-.1951)*(-s319-s1127+s723+s1531)+(-.9808)*(s117+s925-s521-s1329)+(-.5556)*(-s218-s1026+s622+s1430)+(.8315)*(s016+s824-s420-s1228));

data[13]=-.2425*d622+.8422e-1*d1531-.2497*d117-.1069*d319+.2354*d723-.3668e-1*d420-.1489*d824+.1679*d521-.1852*d1329+.1285*d1026-.2260*d1127+.6074e-1*d1430+.2008*d016+.1227e-1*d925+.2144*d218+.2473*d1228;

data[14]=0.25*((-.9802e-1)*(s622-s1430)+(.4714)*(s723-s1531)+(.2903)*(-s521+s1329)+(-.9569)*(s117-s925)+(-.6344)*(-s420+s1228)+(.8819)*(-s319+s1127)+(.7730)*(s016-s824)+(-.9952)*(-s218+s1026));

data[15]=.2260*d622+.1489*d1531-.2473*d319-.2144*d117-.2008*d723+.2497*d420+.1679*d824+.6075e-1*d1329-.2425*d521+.8422e-1*d1026-.3668e-1*d1127-.1069*d1430+.1852*d016+.2354*d218-.1285*d925-.1227e-1*d1228;

data[16]=0.25*((.7071)*(s016+s824+s420+s1228+s218+s1026+s622+s1430-s117-s925-s521-s1329-s319-s1127-s723-s1531));

data[17]=-.1069*d622+.2008*d1531-.3668e-1*d319-.1285*d117+.1489*d723-.1227e-1*d420-.1852*d824+.2425*d1329+.6074e-1*d521+.2473*d1127-.2354*d1026-.2260*d1430+.1679*d016+.2144*d925+.8422e-1*d218-.2497*d1228;

data[18]=0.25*((-.9801e-1)*(s218-s1026)+(-.9952)*(s622-s1430)+(-.7730)*(s420-s1228)+(-.9569)*(-s521+s1329)+(-.4714)*(-s319+s1127)+(.6344)*(s016-s824)+(.2903)*(-s117+s925)+(-.8819)*(-s723+s1531));

data[19]=-.6075e-1*d622+.2354*d1531-.1227e-1*d117+.2260*d319-.8422e-1*d723-.2473*d420+.2008*d824+.1679*d1329+.1852*d521+.2144*d1026-.1069*d1127-.2425*d1430+.1489*d016-.1285*d218-.2497*d925-.3668e-1*d1228;

data[20]=0.25*((-.9808)*(-s319-s1127+s723+s1531)+(-.8315)*(s218+s1026-s622-s1430)+(-.1951)*(-s117-s925+s521+s1329)+(.5556)*(s016+s824-s420-s1228));

data[21]=.2008*d622+.2497*d1531+.1069*d117+.1679*d319+.1227e-1*d723+.6074e-1*d420-.2144*d824-.2354*d521-.8422e-1*d1329-.3668e-1*d1026-.1852*d1127-.1489*d1430+.1285*d016-.2473*d218+.2260*d925+.2425*d1228;

data[22]=0.25*((-.9569)*(s218-s1026)+(-.7730)*(s521-s1329)+(.9952)*(s723-s1531)+(-.8819)*(-s420+s1228)+(-.9802e-1)*(-s319+s1127)+(.4714)*(s016-s824)+(.2903)*(-s622+s1430)+(-.6344)*(-s117+s925));

data[23]=-.2497*d622+.2425*d1531+.2008*d117-.1285*d319+.6075e-1*d723+.2354*d420+.2260*d824-.2473*d1329+.3668e-1*d521+.2144*d1127-.1679*d1026+.1227e-1*d1430+.1069*d016-.1489*d925-.1852*d218+.8422e-1*d1228;

data[24]=0.25*((-.9239)*(-s117-s925-s521-s1329+s319+s1127+s723+s1531)+(.3827)*(s016+s824+s420+s1228-s218-s1026-s622-s1430));

data[25]=.1852*d622+.2144*d1531+.2473*d117-.2425*d319-.1285*d723-.1069*d420-.2354*d824-.1489*d1329+.2008*d521+.6074e-1*d1127+.2497*d1026+.1679*d1430+.8422e-1*d016+.3668e-1*d925+.1227e-1*d218-.2260*d1228;

data[26]=0.25*((-.6344)*(s319-s1127)+(-.9569)*(s420-s1228)+(.9801e-1)*(-s521+s1329)+(.7730)*(s723-s1531)+(-.4714)*(-s218+s1026)+(.2903)*(s016-s824)+(-.8819)*(-s622+s1430)+(-.9952)*(-s117+s925));

data[27]=-.3668e-1*d622+.1679*d1531+.2354*d117-.1227e-1*d319+.1852*d723-.2144*d420+.2425*d824+.1069*d1329-.2260*d521-.1489*d1026-.2497*d1127+.2473*d1430+.6075e-1*d016+.8422e-1*d925+.2008*d218-.1285*d1228;

data[28]=0.25*((.5556)*(s319+s1127-s723-s1531)+(-.8315)*(-s117-s925+s521+s1329)+(-.9808)*(-s218-s1026+s622+s1430)+(.1951)*(s016+s824-s420-s1228));

data[29]=-.1285*d622+.1069*d1531+.2354*d319+.1679*d117-.2260*d723+.1489*d420-.2473*d824+.2497*d1329+.1227e-1*d521+.8422e-1*d1127-.6074e-1*d1026+.2144*d1430+.3668e-1*d016-.1852*d925+.2425*d218+.2008*d1228;

data[30]=0.25*((.9569)*(s319-s1127)+(.6344)*(s622-s1430)+(.8819)*(s521-s1329)+(-.9952)*(-s420+s1228)+(-.7730)*(-s218+s1026)+(.9802e-1)*(s016-s824)+(-.4714)*(-s117+s925)+(-.2903)*(-s723+s1531));

data[31]=.2354*d622+.3668e-1*d1531+.1489*d319+.6075e-1*d117+.2473*d723+.1852*d420+.2497*d824+.1285*d1329+.2144*d521+.2260*d1026+.2008*d1127+.8422e-1*d1430+.1227e-1*d016+.1069*d218+.2425*d925+.1679*d1228;

data=data+32;

}

//列DCT变换

    data = sourcedata;

    for(i=0; i<32; i++)

{

s016=data[0*32]+data[31<<5];

s824=data[16<<5]+data[15<<5];

s420=data[8<<5]+data[23<<5];

s1228=data[24<<5]+data[7<<5];

s218=data[4<<5]+data[27<<5];

    s1026=data[20<<5]+data[11<<5];

s622=data[12<<5]+data[19<<5];

    s1430=data[28<<5]+data[3<<5];

s117=data[2<<5]+data[29<<5];

s925=data[18<<5]+data[13<<5];

s521=data[10<<5]+data[21<<5];

s1329=data[26<<5]+data[5<<5];

s319=data[6<<5]+data[25<<5];

s1127=data[22<<5]+data[9<<5];

s723=data[14<<5]+data[17<<5];

s1531=data[30<<5]+data[1<<5];

    

d016=data[0*32]-data[31<<5];

d824=data[16<<5]-data[15<<5];

d420=data[8<<5]-data[23<<5];

d1228=data[24<<5]-data[7<<5];

d218=data[4<<5]-data[27<<5];

    d1026=data[20<<5]-data[11<<5];

d622=data[12<<5]-data[19<<5];

    d1430=data[28<<5]-data[3<<5];

d117=data[2<<5]-data[29<<5];

d925=data[18<<5]-data[13<<5];

d521=data[10<<5]-data[21<<5];

d1329=data[26<<5]-data[5<<5];

d319=data[6<<5]-data[25<<5];

d1127=data[22<<5]-data[9<<5];

d723=data[14<<5]-data[17<<5];

d1531=data[30<<5]-data[1<<5];

data[0]=.1768*(s1430+s1026+s824+s117+s925+s521+s1228+s016+s420+s1329+s319+s1127+s218+s723+s1531+s622);

data[1<<5]=.8422e-1*d622-.2473*d1531+.2425*d117+.2008*d319+.3668e-1*d723+.1679*d420-.1227e-1*d824-.2144*d1329+.1285*d521-.1489*d1127-.1069*d1026-.2354*d1430+.2497*d016+.2260*d218-.6074e-1*d925-.1852*d1228;

data[2<<5]=0.25*((.8819)*(s117-s925)+(-.7730)*(s622-s1430)+(-.4714)*(s521-s1329)+(.6344)*(s218-s1026)+(-.9801e-1)*(s420-s1228)+(-.9569)*(s723-s1531)+(.2903)*(s319-s1127)+(.9952)*(s016-s824));

data[3<<5]=-.2144*d622-.2260*d1531+.1852*d117-.8422e-1*d319-.1069*d723-.2008*d420+.3668e-1*d824-.2497*d521+.1227e-1*d1329+.2354*d1127+.2425*d1026-.1285*d1430+.2473*d016+.1679*d925+.6075e-1*d218+.1489*d1228;

data[4<<5]=0.25*((-.8315)*(s319+s1127-s723-s1531)+(-.1951)*(s218+s1026-s622-s1430)+(.5556)*(s117+s925-s521-s1329)+(.9808)*(s016+s824-s420-s1228));

data[5<<5]=.2473*d622-.1852*d1531-.2497*d319+.8422e-1*d117+.1679*d723-.1285*d420-.6074e-1*d824+.2260*d1329+.1069*d521-.2008*d1026+.1227e-1*d1127+.3668e-1*d1430+.2425*d016-.1489*d218-.2354*d925+.2144*d1228;

data[6<<5]=0.25*((.9802e-1)*(s117-s925)+(-.4714)*(-s622+s1430)+(-.8819)*(s218-s1026)+(-.9952)*(-s521+s1329)+(-.7730)*(s319-s1127)+(.6344)*(-s723+s1531)+(.9569)*(s016-s824)+(-.2903)*(-s420+s1228));

data[7<<5]=-.1679*d622-.1285*d1531-.6075e-1*d319-.3668e-1*d117-.2144*d723+.2260*d420+.8422e-1*d824+.2008*d1329+.1489*d521-.2425*d1127+.1227e-1*d1026+.1852*d1430+.2354*d016+.2473*d925-.2497*d218-.1069*d1228;

data[8<<5]=0.25*((-.3827)*(s117+s925+s521+s1329-s319-s1127-s723-s1531)+(.9239)*(s016+s824+s420+s1228-s218-s1026-s622-s1430));

data[9<<5]=.1227e-1*d622-.6074e-1*d1531-.1489*d117+.2144*d319+.2425*d723+.8422e-1*d420-.1069*d824-.2473*d521-.3668e-1*d1329+.1285*d1127+.1852*d1026+.2497*d1430+.2260*d016-.2008*d925-.1679*d218-.2354*d1228;

data[10<<5]=0.25*((-.7730)*(s117-s925)+(-.6344)*(s521-s1329)+(-.9569)*(-s622+s1430)+(-.4714)*(s420-s1228)+(-.9952)*(-s319+s1127)+(.2903)*(-s218+s1026)+(.9801e-1)*(-s723+s1531)+(.8819)*(s016-s824));

data[11<<5]=.1489*d622+.1227e-1*d1531+.1852*d319-.2260*d117-.2497*d723-.2425*d420+.1285*d824-.2354*d1329+.8422e-1*d521+.1679*d1127-.2473*d1026+.2008*d1430+.2144*d016+.1069*d925+.3668e-1*d218+.6075e-1*d1228;

data[12<<5]=0.25*((-.1951)*(-s319-s1127+s723+s1531)+(-.9808)*(s117+s925-s521-s1329)+(-.5556)*(-s218-s1026+s622+s1430)+(.8315)*(s016+s824-s420-s1228));

data[13<<5]=-.2425*d622+.8422e-1*d1531-.2497*d117-.1069*d319+.2354*d723-.3668e-1*d420-.1489*d824+.1679*d521-.1852*d1329+.1285*d1026-.2260*d1127+.6074e-1*d1430+.2008*d016+.1227e-1*d925+.2144*d218+.2473*d1228;

data[14<<5]=0.25*((-.9802e-1)*(s622-s1430)+(.4714)*(s723-s1531)+(.2903)*(-s521+s1329)+(-.9569)*(s117-s925)+(-.6344)*(-s420+s1228)+(.8819)*(-s319+s1127)+(.7730)*(s016-s824)+(-.9952)*(-s218+s1026));

data[15<<5]=.2260*d622+.1489*d1531-.2473*d319-.2144*d117-.2008*d723+.2497*d420+.1679*d824+.6075e-1*d1329-.2425*d521+.8422e-1*d1026-.3668e-1*d1127-.1069*d1430+.1852*d016+.2354*d218-.1285*d925-.1227e-1*d1228;

data[16<<5]=0.25*((.7071)*(s016+s824+s420+s1228+s218+s1026+s622+s1430-s117-s925-s521-s1329-s319-s1127-s723-s1531));

data[17<<5]=-.1069*d622+.2008*d1531-.3668e-1*d319-.1285*d117+.1489*d723-.1227e-1*d420-.1852*d824+.2425*d1329+.6074e-1*d521+.2473*d1127-.2354*d1026-.2260*d1430+.1679*d016+.2144*d925+.8422e-1*d218-.2497*d1228;

data[18]=0.25*((-.9801e-1)*(s218-s1026)+(-.9952)*(s622-s1430)+(-.7730)*(s420-s1228)+(-.9569)*(-s521+s1329)+(-.4714)*(-s319+s1127)+(.6344)*(s016-s824)+(.2903)*(-s117+s925)+(-.8819)*(-s723+s1531));

data[19<<5]=-.6075e-1*d622+.2354*d1531-.1227e-1*d117+.2260*d319-.8422e-1*d723-.2473*d420+.2008*d824+.1679*d1329+.1852*d521+.2144*d1026-.1069*d1127-.2425*d1430+.1489*d016-.1285*d218-.2497*d925-.3668e-1*d1228;

data[20<<5]=0.25*((-.9808)*(-s319-s1127+s723+s1531)+(-.8315)*(s218+s1026-s622-s1430)+(-.1951)*(-s117-s925+s521+s1329)+(.5556)*(s016+s824-s420-s1228));

data[21<<5]=.2008*d622+.2497*d1531+.1069*d117+.1679*d319+.1227e-1*d723+.6074e-1*d420-.2144*d824-.2354*d521-.8422e-1*d1329-.3668e-1*d1026-.1852*d1127-.1489*d1430+.1285*d016-.2473*d218+.2260*d925+.2425*d1228;

data[22<<5]=0.25*((-.9569)*(s218-s1026)+(-.7730)*(s521-s1329)+(.9952)*(s723-s1531)+(-.8819)*(-s420+s1228)+(-.9802e-1)*(-s319+s1127)+(.4714)*(s016-s824)+(.2903)*(-s622+s1430)+(-.6344)*(-s117+s925));

data[23<<5]=-.2497*d622+.2425*d1531+.2008*d117-.1285*d319+.6075e-1*d723+.2354*d420+.2260*d824-.2473*d1329+.3668e-1*d521+.2144*d1127-.1679*d1026+.1227e-1*d1430+.1069*d016-.1489*d925-.1852*d218+.8422e-1*d1228;

data[24<<5]=0.25*((-.9239)*(-s117-s925-s521-s1329+s319+s1127+s723+s1531)+(.3827)*(s016+s824+s420+s1228-s218-s1026-s622-s1430));

data[25<<5]=.1852*d622+.2144*d1531+.2473*d117-.2425*d319-.1285*d723-.1069*d420-.2354*d824-.1489*d1329+.2008*d521+.6074e-1*d1127+.2497*d1026+.1679*d1430+.8422e-1*d016+.3668e-1*d925+.1227e-1*d218-.2260*d1228;

data[26<<5]=0.25*((-.6344)*(s319-s1127)+(-.9569)*(s420-s1228)+(.9801e-1)*(-s521+s1329)+(.7730)*(s723-s1531)+(-.4714)*(-s218+s1026)+(.2903)*(s016-s824)+(-.8819)*(-s622+s1430)+(-.9952)*(-s117+s925));

data[27<<5]=-.3668e-1*d622+.1679*d1531+.2354*d117-.1227e-1*d319+.1852*d723-.2144*d420+.2425*d824+.1069*d1329-.2260*d521-.1489*d1026-.2497*d1127+.2473*d1430+.6075e-1*d016+.8422e-1*d925+.2008*d218-.1285*d1228;

data[28<<5]=0.25*((.5556)*(s319+s1127-s723-s1531)+(-.8315)*(-s117-s925+s521+s1329)+(-.9808)*(-s218-s1026+s622+s1430)+(.1951)*(s016+s824-s420-s1228));

data[29<<5]=-.1285*d622+.1069*d1531+.2354*d319+.1679*d117-.2260*d723+.1489*d420-.2473*d824+.2497*d1329+.1227e-1*d521+.8422e-1*d1127-.6074e-1*d1026+.2144*d1430+.3668e-1*d016-.1852*d925+.2425*d218+.2008*d1228;

data[30<<5]=0.25*((.9569)*(s319-s1127)+(.6344)*(s622-s1430)+(.8819)*(s521-s1329)+(-.9952)*(-s420+s1228)+(-.7730)*(-s218+s1026)+(.9802e-1)*(s016-s824)+(-.4714)*(-s117+s925)+(-.2903)*(-s723+s1531));

data[31<<5]=.2354*d622+.3668e-1*d1531+.1489*d319+.6075e-1*d117+.2473*d723+.1852*d420+.2497*d824+.1285*d1329+.2144*d521+.2260*d1026+.2008*d1127+.8422e-1*d1430+.1227e-1*d016+.1069*d218+.2425*d925+.1679*d1228;

data ++;

}

return sourcedata;

}