wavemain.cpp

this is source code for image compression following jpeg2000 staandard

/*---------------------------------------------------------------------------*/
// Image Compression Toolbox v1.2
// written by
// Satish Kumar S
// [e-mail]
//
// Copyright 1999 Satish Kumar S
//
// Permission is granted to use this software for research purposes as
// long as this notice stays attached to this software.
/*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include "Image.h"
#include "Wavelet.h"
#include "aricoder.h"
#include "quantizer.h"

#define sqrt2 1.414213
#define NSTEPS 5
#define NBANDS (NSTEPS*3+1)
#define MAXREAL 1e+100

double HarrCoeffs[] = {1.0/sqrt2, 1.0/sqrt2};
double AntoniniAnalysis[] =  {  3.782845550699535e-02,
			     -2.384946501937986e-02,
			     -1.106244044184226e-01,
			      3.774028556126536e-01,
			      8.526986790094022e-01,
			      3.774028556126537e-01,
			     -1.106244044184226e-01,
			     -2.384946501937986e-02,
			      3.782845550699535e-02 };
double AntoniniSynthesis [] = { -6.453888262893856e-02,
			      -4.068941760955867e-02,
			       4.180922732222124e-01,
			       7.884856164056651e-01,
			       4.180922732222124e-01,
			      -4.068941760955867e-02,
			      -6.453888262893856e-02 };
//----------------
struct stBlock
{
    int x, y;
    int width, height;
    int ncoeffs, nonzerocoeffs;
    unsigned int *symbols;
    unsigned int zeroTreeSymbol;
    double *coeffs;
    double maxcoeff, mincoeff, meancoeff;
    UniformQuant *uq;
};
//----------------
Filter fHarrAnalysis(2, 0, HarrCoeffs);
Wavelet wHarr(&fHarrAnalysis);

//----------------
Filter fAntoniniAnalysis(9, -4, AntoniniAnalysis);
Filter fAntoniniSynthesis(7, -3, AntoniniSynthesis);
Wavelet wAntonini(&fAntoniniAnalysis, &fAntoniniSynthesis);

double nBitsForZero =0;

//------------------------------------------------------------//
//------------------------------------------------------------//
void AllocSubBands(Image &image, struct stBlock *subband, int nsteps)
{
    int prevllWidth  = image.width;
    int prevllHeight = image.height;
    int i;

    for(i=nsteps-1; i>=0; i--)
    {
        int llWidth  = (prevllWidth+1)/2;
        int llHeight = (prevllHeight+1)/2;
        int hhWidth  = prevllWidth - llWidth;
        int hhHeight = prevllHeight - llHeight;

        // top right block
        subband[i*3+1].x = llWidth;
        subband[i*3+1].y = 0;
        subband[i*3+1].width = hhWidth;
        subband[i*3+1].height = llHeight;

        // bottom left block
        subband[i*3+2].x = 0;
        subband[i*3+2].y = llHeight;
        subband[i*3+2].width = llWidth;
        subband[i*3+2].height = hhHeight;

        // bottom right block
        subband[i*3+3].x = llWidth;
        subband[i*3+3].y = llHeight;
        subband[i*3+3].width = hhWidth;
        subband[i*3+3].height = hhHeight;

        prevllWidth = llWidth;
        prevllHeight = llHeight;
    }
    // ll coefficient band.
    subband[0].x = 0;
    subband[0].y = 0;
    subband[0].width = prevllWidth;
    subband[0].height = prevllHeight;

    // alloc mem & copy the coefficients to individual arrays.
    for(i=0; i<NBANDS; i++)
    {
        subband[i].ncoeffs = subband[i].width * subband[i].height;
        subband[i].coeffs = new double[ subband[i].ncoeffs ];
        subband[i].symbols = new unsigned int[ subband[i].ncoeffs ];
        assert(subband[i].coeffs && subband[i].symbols);
        subband[i].uq = NULL;
    }
}

void FreeSubbands(struct stBlock *subband, int nbands)
{
    for(int i=0; i<nbands; i++)
    {
        delete[] subband[i].coeffs;
        delete[] subband[i].symbols;
    }
}
//------------------------------------------------------------//
//------------------------------------------------------------//
void InitSubBands(Image &image, struct stBlock *subband)
{
    int i,j;
    
    // copy the data to the indiv. subbands find the max & min values in each subband.
    for(i=0; i<NBANDS; i++)
    {
        image.GetData(subband[i].x, subband[i].y, subband[i].width, subband[i].height, 
            subband[i].coeffs);

        // find the max & min in the coefficient array.
        double maxcoeff = -MAXREAL;
        double mincoeff = MAXREAL;
        double *coeffs  = subband[i].coeffs;
        subband[i].nonzerocoeffs = 0;
        for(j=subband[i].ncoeffs; j>0; j--)
        {
            if (maxcoeff < *coeffs) maxcoeff = *coeffs;
            else if (mincoeff > *coeffs) mincoeff = *coeffs;
            if (*coeffs != 0.0) subband[i].nonzerocoeffs++;
            coeffs++;
        }
        subband[i].maxcoeff = maxcoeff;
        subband[i].mincoeff = mincoeff;
        // except for lowpass subband, assume the pdfs are centered around zero.
        if (i > 0)
        {
            if (fabs(maxcoeff) > fabs(mincoeff))
                subband[i].mincoeff = -maxcoeff;
            else
                subband[i].maxcoeff = -mincoeff;
        }
        subband[i].meancoeff = (subband[i].maxcoeff+subband[i].mincoeff)/2;
    }
}

//------------------------------------------------------------//
//------------------------------------------------------------//
void CalcRateDistortion(struct stBlock &subband, int nBits, double &Rate, double &Distortion, ArithEncoder &arienc)
{
    int nLevels = (1<<nBits);

    // initialize the quantizer.
    UniformQuant uq(subband.maxcoeff, subband.mincoeff, nLevels);

    // calculate the rate & distortion.
    double *coeffs = subband.coeffs;
    Distortion = 0.0;
    Rate = 0.0;

    for(int ncoeffs = subband.ncoeffs; ncoeffs>0; ncoeffs--)
    {
        unsigned int symbol = uq.Symbol(*coeffs);
        double recons = uq.Value(symbol);

        Distortion += pow((*coeffs - recons),2);
        coeffs++;
        if (nBits)
        {
            double tRate = arienc.EncodingCost(symbol);
            //double tRate=nBits;
            if (*coeffs==0.0) nBitsForZero += tRate;
            Rate += tRate;
        }
    }
    Distortion /= subband.ncoeffs;
    Rate /= subband.ncoeffs;
}

//------------------------------------------------------------//
// BitAllocation is done using a modified version of the algo
// given in "Optimal Bit Allocation via the General BFOS algo"
//           by Eve A. Riskin
//------------------------------------------------------------//
void CalcBitAllocation(struct stBlock *subband, int *B, int M, int q, double TargetRate)
{
    // calculate the bitallocation for each subband/class.
    double *p;          // the probability of each class
    double *r;          // the rate for each class for a given no of bits.
    double *d;          // the distortion for each class for a given no. of bits.
    int *allow;
    int i,j;

    p = new double[M];
    r = new double[M * (q+1)];      // 2-d array
    d = new double[M * (q+1)];      // 2-d array
    allow = new int[M];

    // initialize all p(i)'s
    double totcoeffs = 0;
    for(i=0; i<M; i++)
        totcoeffs += subband[i].ncoeffs;
    for(i=0; i<M; i++)
        p[i] = ((double)subband[i].ncoeffs)/totcoeffs;

    // first calc all di(j)'s
    for(i=0; i<M; i++)
    {
        printf("For class %d\r",i);
        for(j=0; j<=q; j++)
        {
            // Initialize the arithmetic coder with a dummy bitstream.
            ArithEncoder arienc(NULL);
            arienc.SetNumSyms((1<<j));
            arienc.InitEncode();

            double tRate, tDist;
    
            CalcRateDistortion(subband[i], j, tRate, tDist, arienc);
            r[i*(q+1)+j] = tRate;
            d[i*(q+1)+j] = tDist;
        }
    }

    // initialise all B(i)'s
    for(i=0; i<M; i++)
        B[i] = q;

    double D;       // overall average distortion.
    double R;       // overall average rate.
    double *lowestSi = new double[M];         // lowest Si in each class.
    int *lowestSiClass_n = new int[M];  // the corresponding n.
    double TargetRateInBytes = TargetRate * totcoeffs;  // assuming 1pixel = 1byte in src.
    for(i=0; i<M; i++) allow[i]=1;

    do
    {
        // find all Si's and the class(es) for which Si is the lowest.
        double lowestSiInAll = MAXREAL;
        for(i=0; i<M; i++)
            lowestSi[i] = MAXREAL;
        for(i=0; i<M; i++)
        {
            for(int n=1; n<=B[i]; n++)
            {
                double Si = ( d[i*(q+1) + B[i]-n] - d[i*(q+1) + B[i]] )/(double)n;
                if(Si <= lowestSi[i] && allow[i])
                {
                    lowestSiClass_n[i] = n;
                    lowestSi[i] = Si;
                    if (Si <= lowestSiInAll)
                        lowestSiInAll = Si;
                }
            }
        }
        if (lowestSiInAll == MAXREAL)       // cant find one??
            break;
        for(i=0; i<M; i++)
            if (lowestSi[i] == lowestSiInAll)
                B[i] = B[i] - lowestSiClass_n[i];   // remove bits from that class.

        // calculate the rate and distortion of the new code.
        D=0; R=0;
        for(i=0; i<M; i++)
        {
            D += p[i] * d[i*(q+1) + B[i]];
            R += r[i*(q+1) + B[i]] * subband[i].ncoeffs;
        }
        if (R<TargetRateInBytes-50)
        {
            for(i=0; i<M; i++)
                if (lowestSi[i] == lowestSiInAll)
                {
                    B[i] = B[i] + lowestSiClass_n[i];
                    allow[i] = 0;
                }
            R = TargetRateInBytes+1;
        }
        else
            for(i=0; i<M; i++) allow[i]=1;
    }while(R>TargetRateInBytes);

    // print the statistics.
    double totdist = 0,totrate=0;
    for(i=0; i<M; i++)
    {
        printf("Class % 2d\t%2d bits\n", i, B[i]);
        totdist += d[i*(q+1)+B[i]] * subband[i].ncoeffs;
        totrate += r[i*(q+1)+B[i]] * subband[i].ncoeffs;
    }
    double rms = sqrt(totdist/totcoeffs);
    double psnr = 20.0 * log(255.0/rms)/log(10.0);

    printf("Requested Rate = %d\n",(int)TargetRateInBytes);
    printf("Actual rate = %d\n",(int)totrate);
    printf("RMS error = %11.2f\n",(float)rms);
    printf("PSNR = %11.2f\n",(float)psnr);

    delete[] allow;
    delete[] p;
    delete[] r;
    delete[] d;
}

//------------------------------------------------------------//
//------------------------------------------------------------//
void normal_decode(char *fname)
{
    FILE *fp = fopen(fname,"rb");
    BitStream bs(fp);

    int imwidth = bs.ReadBits(16);
    int imheight = bs.ReadBits(16);
    int nsteps = bs.ReadBits(16);
    int nbands = nsteps*3+1;
    
    Image image;
    image.InitEmptyImage(imwidth, imheight);

    struct stBlock *subband = new struct stBlock[nbands];
    AllocSubBands(image, subband, nsteps);

    for(int i=0; i<nbands; i++)
    {
        // read in the quantizer data.
        UniformQuant *uq = new UniformQuant();
        uq->ReadHeader(bs);
        subband[i].uq = uq;
        int nlevels = uq->nLevels;

        unsigned int *symbols = subband[i].symbols;
        // only a single level? then no bits will be stored, and the symbol is always 0.
        if (nlevels<=1)
        {
            for(int j=0; j<subband[i].ncoeffs; j++)
                *symbols++ = 0;
        }
        else    // normal case.
        {
            // Setup the arithmetic decoder.
            ArithDecoder adc(&bs);
            adc.SetNumSyms(nlevels);
            adc.InitDecode();

            // now readin the symbols.
            for(int j=0; j<subband[i].ncoeffs; j++)
            {
                *symbols = adc.Decode();
                symbols++;
            }
        }
        // and dequantize them.
        double *coeffs = subband[i].coeffs;
        symbols = subband[i].symbols;

        for(int j=0; j<subband[i].ncoeffs; j++)
        {
            *coeffs++ = uq->Value(*symbols++);
        }

        image.SetData(subband[i].x, subband[i].y, subband[i].width, subband[i].height, subband[i].coeffs);
    }

    // close the file.
    fclose(fp);

    // take inverse transform.
    Wavelet *cw = &wAntonini;
    cw->Inverse(image, nsteps);
    image.SaveRawFile("wletn.raw");

    FreeSubbands(subband, nbands);
    for(i=0; i<nbands; i++)
        delete subband[i].uq;
    delete[] subband;
}

//------------------------------------------------------------//
//------------------------------------------------------------//
void normal_encode(char *fname)
{
    int i,j;
    Image image;