//小波基
const double hCoef[10][20] =
{
    { .707106781187,  .707106781187},
    { .482962913145,  .836516303738,  .224143868042, -.129409522551 },
    { .332670552950,  .806891509311,  .459877502118, -.135011020010, -.085441273882,  .035226291882 },
    { .230377813309,  .714846570553,  .630880767930, -.027983769417,-.187034811719,  .030841381836,  .032883011667, -.010597401785 },
    { .160102397974,  .603829269797,  .724308528438,  .138428145901, -.242294887066,-.032244869585,  .077571493840, -.006241490213, -.012580751999,  .003335725285 },
    { .111540743350,  .494623890398,  .751133908021,  .315250351709, -.226264693965,
    -.129766867567,  .097501605587,  .027522865530, -.031582039318,  .000553842201,
    .004777257511, -.001077301085 },
    { .077852054085,  .396539319482,  .729132090846,  .469782287405, -.143906003929,
    -.224036184994,  .071309219267,  .080612609151, -.038029936935, -.016574541631,
    .012550998556,  .000429577973, -.001801640704,  .000353713800 },
    { .054415842243,  .312871590914,  .675630736297,  .585354683654, -.015829105256,
    -.284015542962,  .000472484574,  .128747426620, -.017369301002, -.044088253931,
    .013981027917,  .008746094047, -.004870352993, -.000391740373,  .000675449406,
    -.000117476784 },
    { .038077947364,  .243834674613,  .604823123690,  .657288078051,  .133197385825,
    -.293273783279, -.096840783223,  .148540749338,  .030725681479, -.067632829061,
    .000250947115,  .022361662124, -.004723204758, -.004281503682,  .001847646883,
    .000230385764, -.000251963189,  .000039347320 },
    { .026670057901,  .188176800078,  .527201188932,  .688459039454,  .281172343661,
    -.249846424327, -.195946274377,  .127369340336,  .093057364604, -.071394147166,
    -.029457536822,  .033212674059,  .003606553567, -.010733175483,  .001395351747,
    .001992405295, -.000685856695, -.000116466855,  .000093588670, -.000013264203 }
};
//一維小波變換
BOOL DWTStep_1D(double* pDbSrc, int nCurLevel,int nInv, int nStep,int nSupp)
{
    double s = sqrt((double)2);
    // 獲得小波基的指標
    double* h = (double*)hCoef[nSupp-1];
    // 確認當前層數有效
    ASSERT(nCurLevel>=0);
    // 計算當前層數的長度
    int CurN = 1<<nCurLevel;
    if (nInv) CurN <<= 1;
    // 確認所選擇的小波基和當前層數的長度有效
    if (nSupp<1 || nSupp>10 || CurN<2*nSupp)
        return FALSE;
    // 分配臨時記憶體用於存放結果
    double *ptemp = new double[CurN];
    if (!ptemp) return FALSE;
    double    s1, s2;
    int    Index1, Index2;
    // 判斷是進行DWT還是IDWT
    if (!nInv)
    {    // DWT
        Index1=0;
        Index2=2*nSupp-1;
        // 進行卷積，其中s1為低頻部分，s2為高頻部分的結果
        for (int i=0; i<CurN/2; i++)
        {    
            s1 = s2 = 0;
            double t = -1;
            for (int j=0; j<2*nSupp; j++, t=-t)
            {
                s1 += h[j]*pDbSrc[(Index1 & CurN-1) * nStep];
                s2 += t*h[j]*pDbSrc[(Index2 & CurN-1) * nStep];
                Index1++;
                Index2--;
            }
            // 將結果存放在臨時記憶體中
            ptemp[i] = s1/s;
            ptemp[i+CurN/2] = s2/s;
            Index1 -= 2*nSupp;
            Index2 += 2*nSupp;
            Index1 += 2;
            Index2 += 2;
        }
    }
    // 否則進行IDWT
    else
    {    // IDWT
        Index1 = CurN/2;
        Index2 = CurN/2-nSupp+1;
        // 進行卷積，其中其中s1為低頻部分，s2為高頻部分的結果
        for (int i=0; i<CurN/2; i++)
        {
            s1 = s2 = 0;
            int Index3 = 0;
            for (int j=0; j<nSupp; j++)
            {
                s1 += h[Index3]*pDbSrc[(Index1 & CurN/2-1) * nStep]
                +h[Index3+1]*pDbSrc[((Index2 & CurN/2-1) + CurN/2) * nStep];
                s2 += h[Index3+1]*pDbSrc[(Index1 & CurN/2-1) * nStep]
                -h[Index3]*pDbSrc[((Index2 & CurN/2-1) + CurN/2) * nStep];

                Index3+=2;
                Index1--,        Index2++;
            }
            // 將結果存入臨時記憶體
            ptemp[2*i] = s1*s;
            ptemp[2*i+1] = s2*s;
            Index1 += nSupp;
            Index2 -= nSupp;
            Index1++;
            Index2++;
        }
    }
    // 將結果存入源影象中
    for (int i=0; i<CurN; i++)
        pDbSrc[i*nStep] = ptemp[i];
    // 釋放臨時記憶體，並返回
    delete[] ptemp;
    return TRUE;
}
//二維小波變換
BOOL DWTStep_2D(double* pDbSrc, int nCurWLevel, int nCurHLevel,
                int nMaxWLevel, int nMaxHLevel, int nInv, int nStep, int nSupp)
{
    // 計算影象的長度和寬度（2次冪對齊）
    int W = 1<<nMaxWLevel, H = 1<<nMaxHLevel;
    // 計算當前分解的影象的長度和寬度
    int CurW = 1<<nCurWLevel, CurH = 1<<nCurHLevel;
    // 判斷是進行DWT還是IDWT
    if (!nInv)
    {    
        int i = 0;
        // 對行進行一維DWT
        for (i=0; i<CurH; i++)
            if (!DWTStep_1D(pDbSrc+(int)i*W*nStep, nCurWLevel, nInv, nStep, nSupp)) return FALSE;
        // 對列進行一維DWT
        for (i=0; i<CurW; i++)
            if (!DWTStep_1D(pDbSrc+i*nStep, nCurHLevel, nInv, W*nStep, nSupp)) return FALSE;
    }
    // 否則進行IDWT
    else
    {
        // 計算當前變換的影象的長度和寬度
        CurW <<= 1;
        CurH <<= 1;
        int i = 0;
        // 對列進行IDWT
        for (i=0; i<CurW; i++)
            if (!DWTStep_1D(pDbSrc+i*nStep, nCurHLevel, nInv, W*nStep, nSupp)) return FALSE;
        // 對行進行IDWT
        for (i=0; i<CurH; i++)
            if (!DWTStep_1D(pDbSrc+(int)i*W*nStep, nCurWLevel, nInv, nStep, nSupp)) return FALSE;
    }
    // 返回
    return TRUE;
}
//資料轉換
BYTE FloatToByte(double f)
{
    if (f<=0) return (BYTE)0;
    else if (f>=255) return (BYTE)255;
    else return (BYTE)(f+0.5);
}
//資料轉換
char FloatToChar(double f)
{
    if (f>=0)
        if (f>=127.0)
            return (char)127;
        else return (char)(f+0.5);
    else
        if (f<=-128)
            return (char)-128;
        else return -(char)(-f+0.5);
}
//資料轉換
int Log2(int n)
{
    int rsl = 0;
    while (n >>= 1) rsl++;
    return rsl;
}
//影象小波變換
BOOL DIBDWTStep(LPSTR lpDIBBits,double*m_pDbImage, int nWidth,int nHeight, int nInv,int m_nDWTCurDepth,int m_nSupp)
{
    // 迴圈變數
    int i, j;
    // 獲取變換的最大層數
    int nMaxWLevel = Log2(nWidth);
    int nMaxHLevel = Log2(nHeight);
    int nMaxLevel;
    if (nWidth == 1<<nMaxWLevel && nHeight == 1<<nMaxHLevel)
        nMaxLevel = min(nMaxWLevel, nMaxHLevel);
    // 臨時變數
    double    *pDbTemp;
    BYTE    *pBits;
    // 如果小波變換的儲存記憶體還沒有分配，則分配此記憶體
    if(!m_pDbImage){            
        m_pDbImage = new double[nWidth*nHeight];
        if (!m_pDbImage)    return FALSE;
        // 將影象資料放入m_pDbImage中
        for (j=0; j<nHeight; j++)
        {
            pDbTemp = m_pDbImage + j*nWidth;
            pBits = (unsigned char *)lpDIBBits + (nHeight-1-j)*nWidth;        
            for (i=0; i<nWidth; i++)
                pDbTemp[i] = pBits[i];
        }
    }
    // 進行小波變換（或反變換）
    if (!DWTStep_2D(m_pDbImage, nMaxWLevel-m_nDWTCurDepth, nMaxHLevel-m_nDWTCurDepth,
        nMaxWLevel, nMaxHLevel, nInv, 1, m_nSupp))
        return FALSE;
    // 如果是反變換，則當前層數減1
    if (nInv)
        m_nDWTCurDepth --;
    // 否則加1
    else
        m_nDWTCurDepth ++;
    // 然後，將資料拷貝回原CDib中，並進行相應的資料轉換
    int lfw = nWidth>>m_nDWTCurDepth, lfh = nHeight>>m_nDWTCurDepth;
    for (j=0; j<nHeight; j++)
    {
        pDbTemp = m_pDbImage + j*nWidth;
        pBits = (unsigned char *)lpDIBBits + (nHeight-1-j)*nWidth;
        for (i=0; i<nWidth; i++)
        {
            if (j<lfh && i<lfw)
                pBits[i] = FloatToByte(pDbTemp[i]);
            else
                pBits[i] = BYTE(FloatToChar(pDbTemp[i]) ^ 0x80);            
        }
    }
    // 返回
    return TRUE;
}