wavelet.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 "Image.h"
#include "Wavelet.h"

//------------------------------------------------------------------
// Implementation of the Wavelet class.
//  some functions are taken from the 'Wavelet Image compression Kit'
//  written and distributed by Goeff Davis 
//  "http://research.microsoft.com/~geoffd"
//------------------------------------------------------------------

#define min(a,b) (((a)<(b))?(a):(b))

/*---------------------------------------------------------------------------*/
// Do symmetric extension of data using prescribed symmetries
//   Original values are in output[npad] through output[npad+size-1]
//   New values will be placed in output[0] through output[npad] and in
//      output[npad+size] through output[2*npad+size-1] (note: end values may
//      not be filled in)
//   left_ext = 1 -> extension at left bdry is   ... 3 2 1 | 0 1 2 3 ...
//   left_ext = 2 -> extension at left bdry is ... 3 2 1 0 | 0 1 2 3 ...
//   right_ext = 1 or 2 has similar effects at the right boundary
//
//   symmetry = 1  -> extend symmetrically
//   symmetry = -1 -> extend antisymmetrically
void Wavelet::SymmetricExtension(double *output, int size, int left_ext, int
				   right_ext, int symmetry, int npad)
{
    int i;
    int first = npad, last = npad + size-1;

    if (symmetry == -1)
    {
        if (left_ext == 1)
          output[--first] = 0;
        if (right_ext == 1)
          output[++last] = 0;
    }

    int originalFirst = first;
    int originalLast = last;
    int originalSize = originalLast-originalFirst+1;

    int period = 2 * (last - first + 1) - (left_ext == 1) - (right_ext == 1);

    if (left_ext == 2)
        output[--first] = symmetry * output[originalFirst];
    if (right_ext == 2)
        output[++last] = symmetry * output[originalLast];

    // extend left end
    int nextend = min (originalSize-2, first);
    for (i = 0; i < nextend; i++)
        output[--first] = symmetry * output[originalFirst+1+i];

    // should have full period now -- extend periodically
    while (first-- > 0)
        output[first] = output[first+period];

    // extend right end
    nextend = min (originalSize-2, 2*npad+size-1 - last);
    for (i = 0; i < nextend; i++)
        output[++last] = symmetry * output[originalLast-1-i];

    // should have full period now -- extend periodically
    while (last++ < 2*npad+size-1)
        output[last] = output[last-period];
}

/*---------------------------------------------------------------------------*/
// Do periodic extension of data using prescribed symmetries
//   Original values are in output[npad] through output[npad+size-1]
//   New values will be placed in output[0] through output[npad] and in
//      output[npad+size] through output[2*npad+size-1] (note: end values may
//      not be filled in)
void Wavelet::PeriodicExtension(double *output, int size, int npad)
{
    int first = npad, last = npad + size-1;

    // extend left periodically
    while (first-- > 0)
        output[first] = output[first+size];

    // extend right periodically
    while (last++ < 2*npad+size-1)
        output[last] = output[last-size];
}

Filter::Filter(int tsize, int tstartIndex, double *tcoeffs)
{
    coeffs = new double[tsize];
    if (!coeffs) 
        return;

    size = tsize;
    startIndex = tstartIndex;

    if (tcoeffs)
        for(int i=0;i<tsize;i++)
            coeffs[i] = tcoeffs[i];
}

Wavelet::Wavelet(Filter *aLow, Filter *sLow)
{
    int i, sign;
    anaLow = new Filter(aLow->size, aLow->startIndex, aLow->coeffs);

    // If no synthesis coeffs are given, assume wavelet is orthogonal.
    if (sLow == NULL)
    {
        // For orthogonal wavelets, compute the high pass filter using
        // the relation g_n = (-1)^n h_{1-n}^*
        // (or equivalently g_{1-n} = (-1)^{1-n} h_n^*).
        anaHigh = new Filter(anaLow->size, 2 - anaLow->size - anaLow->startIndex);

        // compute (-1)^(1-n) for n=startIndex
        sign = (anaLow->startIndex % 2)?1:-1;

        for(i=0; i<anaLow->size; i++)
        {
            anaHigh->coeffs[1 - i - anaLow->startIndex - anaHigh->startIndex] = 
                sign * anaLow->coeffs[i];
            sign *= -1;
        }

        // copy the analysis filters to the synthesis filters.
        synLow = new Filter(anaLow->size, anaLow->startIndex, anaLow->coeffs);
        synHigh = new Filter(anaHigh->size, anaHigh->startIndex, anaHigh->coeffs);
    }
    else    // assume wavelet is biorthogonal.
    {
        synLow = new Filter(sLow->size, sLow->startIndex, sLow->coeffs);

        // For orthogonal wavelets, compute the high frequency filter using
        // the relation g_n = (-1)^n complement (h~_{1-n}) and
        //              g~_n = (-1)^n complement (h_{1-n})
        // (or equivalently g_{1-n} = (-1)^{1-n} complement (h~_n))

        anaHigh = new Filter(synLow->size, 2 - synLow->size - synLow->startIndex);
        synHigh = new Filter(anaLow->size, 2 - anaLow->size - anaLow->startIndex);

        // fill in coeffs for analysis high.
        sign = (synLow->startIndex % 2)?1:-1;
        for(i=0; i<synLow->size; i++)
        {
            anaHigh->coeffs[1 - i - synLow->startIndex - anaHigh->startIndex] =
                sign * synLow->coeffs[i];
            sign *= -1;
        }

        // fill in coeffs for synthesis high.
        sign = (anaLow->startIndex % 2)?1:-1;
        for(i=0; i<anaLow->size; i++)
        {
            synHigh->coeffs[1 - i - anaLow->startIndex - synHigh->startIndex] =
                sign * anaLow->coeffs[i];
            sign *= -1;
        }
    }
}

int Wavelet::Transform(Image &image, int nsteps)
{
    int i, j, k;
    double *indata, *outdata, *imagedata;
    int *lpwidtharr, *lpheightarr;

    // Find out how much we need to pad (the greater of the two filter lengths-1).
    int npad = anaLow->size;
    if (npad < anaHigh->size)
        npad = anaHigh->size;

    // allocate data for the largest side.
    if (image.width > image.height)
    {
        indata = new double[image.width + npad * 2 + image.width];
        outdata = indata + image.width + npad * 2;
    }
    else
    {
        indata = new double[image.height + npad * 2 + image.height];
        outdata = indata + image.height + npad * 2;
    }
    if (!indata)
        return 1;

    // allocate space to store the size of the lowpass regions in each step.
    lpwidtharr = new int[nsteps *2];
    if (!lpwidtharr)
    {
        delete[] indata;
        return 1;
    }
    lpheightarr = lpwidtharr + nsteps;

    // and calculate the sizes.
    lpwidtharr[nsteps-1]  = image.width;
    lpheightarr[nsteps-1] = image.height;
    for(i=nsteps-2; i>=0; i--)
    {
        lpwidtharr[i] = (lpwidtharr[i+1]+1)/2;
        lpheightarr[i]= (lpheightarr[i+1]+1)/2;
    }
    
    // now do the transform.
    for(k=nsteps-1; k>=0; k--)
    {
        int lpwidth  = lpwidtharr[k];
        int lpheight = lpheightarr[k];

        // first apply row-wise.
        imagedata = image.data;
        for(i=0; i<lpheight; i++)
        {
            // copy 1 row worth of indata to the buffer.
            for(j=0; j<lpwidth; j++)
                indata[npad+j] = imagedata[j];

            TransformStep(indata, lpwidth, npad, anaLow, anaHigh, outdata);
            
            // copy back
            for(j=0; j<lpwidth; j++)
                imagedata[j] = outdata[j];

            imagedata += image.width;
        }

        // then apply column-wise.
        imagedata = image.data;
        for(i=0; i<lpwidth; i++)
        {
            // copy 1 row worth of indata to the buffer.
            for(j=0; j<lpheight; j++)
                indata[npad+j] = imagedata[j*image.width];

            TransformStep(indata, lpheight, npad, anaLow, anaHigh, outdata);

            // copy back
            for(j=0; j<lpheight; j++)
                imagedata[j * image.width] = outdata[j];

            imagedata += 1;
        }
    }
    return 0;
}

void Wavelet::TransformStep(double *data, int datasize, int npad, Filter *lpf, Filter *hpf, double *dest)
{
    int length, start;
    double *coeffs;
    int i,j;

    // fill in the left and right padding samples.
    int leftext, rightext;
    leftext = (anaLow->size % 2) ? 1 : 2;
    rightext = leftext;
    SymmetricExtension(data, datasize, leftext, rightext, 1, npad);

    // now apply the lowpass filter.
    length = lpf->size;
    start  = lpf->startIndex;
    coeffs = lpf->coeffs;

    int lpsize = (datasize+1)/2;
    for(i=0; i<lpsize; i++)
    {
        dest[i] = 0;
        for(j=0; j<length; j++)
            dest[i] += data[npad + i*2 + start + j] * coeffs[j];
    }
    // now apply the highpass filter.
    length = hpf->size;
    start  = hpf->startIndex;
    coeffs = hpf->coeffs;
    
    for(i=lpsize; i<datasize; i++)
    {
        dest[i] = 0;
        for(j=0; j<length; j++)
            dest[i] += data[npad + (i-lpsize)*2 + start + j] * coeffs[j];
    }
}

int Wavelet::Inverse(Image &image, int nsteps)
{
    int i, j, k;
    double *indata, *outdata, *imagedata;
    int *lpwidtharr, *lpheightarr;

    // Find out how much we need to pad (the greater of the two filter lengths-1).
    int npad = synLow->size;
    if (npad < synHigh->size)
        npad = synHigh->size;

    // allocate data for the largest side.
    if (image.width > image.height)
    {
        indata = new double[(image.width + npad * 2) * 2];
        outdata = indata + image.width + npad * 2;
    }
    else
    {
        indata = new double[(image.height + npad * 2) * 2];
        outdata = indata + image.height + npad * 2;
    }
    if (!indata)
        return 1;

    // allocate space to store the size of the lowpass regions in each step.
    lpwidtharr = new int[nsteps *2];
    if (!lpwidtharr)
    {
        delete[] indata;
        return 1;
    }
    lpheightarr = lpwidtharr + nsteps;

    // and calculate the sizes.
    lpwidtharr[nsteps-1]  = image.width;
    lpheightarr[nsteps-1] = image.height;
    for(i=nsteps-2; i>=0; i--)
    {
        lpwidtharr[i] = (lpwidtharr[i+1]+1)/2;
        lpheightarr[i]= (lpheightarr[i+1]+1)/2;
    }
    
    // now do the transform.
    for(k=0; k<=nsteps-1; k++)
    {
        int lpwidth  = lpwidtharr[k];
        int lpheight = lpheightarr[k];

        // apply column-wise.
        imagedata = image.data;
        for(i=0; i<lpwidth; i++)
        {
            // copy 1 row worth of indata to the buffer.
            for(j=0; j<lpheight; j++)
                indata[npad+j] = imagedata[j*image.width];

            InverseStep(indata, lpheight, npad, synLow, synHigh, outdata);

            // copy back
            for(j=0; j<lpheight; j++)
                imagedata[j * image.width] = outdata[npad + j];

            imagedata += 1;
        }

        // apply row-wise.
        imagedata = image.data;
        for(i=0; i<lpheight; i++)
        {
            // copy 1 row worth of indata to the buffer.
            for(j=0; j<lpwidth; j++)
                indata[npad+j] = imagedata[j];

            InverseStep(indata, lpwidth, npad, synLow, synHigh, outdata);
            
            // copy back
            for(j=0; j<lpwidth; j++)
                imagedata[j] = outdata[npad + j];

            imagedata += image.width;
        }
    }
    return 0;
}

void Wavelet::InverseStep(double *data, int datasize, int npad, Filter *lpf, Filter *hpf, double *dest)
{
    int length;
    double *coeffs, *tempin;
    int i,j;

    int lpsize = (datasize+1)/2;
    int hpsize = datasize - lpsize;

    tempin = new double[npad*2+datasize];

    // copy the lowpass samples.
    for(i=0;i<lpsize;i++)
        tempin[npad+i] = data[npad+i];

    // fill in the padding samples.
    int leftext, rightext;
    leftext = (anaLow->size % 2) ? 1 : 2;               // odd filter length? leftext = 1
    rightext = (datasize % 2) ? 1 : 2;                  // odd data length? rightext = 1
    SymmetricExtension(tempin, lpsize, leftext, rightext, 1, npad);

    // initialise the result to zero.
    for(i=0; i<datasize+npad*2;i++)
        dest[i] = 0.0;

    // apply the lowpass filter.
    length = lpf->size;
    coeffs = lpf->coeffs;

    int firstIndex = lpf->startIndex;
    int lastIndex = length - 1 + firstIndex;

    for (i = -lastIndex/2; i <= (datasize-1-firstIndex)/2; i++)
        for (j = 0; j < length; j++)
            dest[npad + 2*i + firstIndex + j] += tempin[npad+i] * coeffs[j];

    // now to the highpass samples.
    for(i=0; i<hpsize; i++)
        tempin[npad+i] = data[npad+lpsize+i];

    // fill in the padding samples to the left and right.
    leftext = 2;
    int symmetry;
    if (anaLow->size % 2)
    {
        rightext = (datasize % 2) ? 2 : 1;
        symmetry = 1;
    }
    else
    {
        rightext = (datasize % 2) ? 1 : 2;
        symmetry = -1;
    }
    SymmetricExtension(tempin, hpsize, leftext, rightext, symmetry, npad);

    // apply the highpass filter.
    length = hpf->size;
    coeffs = hpf->coeffs;

    firstIndex = hpf->startIndex;
    lastIndex = length - 1 + firstIndex;

    for (i = -lastIndex/2; i <= (datasize-1-firstIndex)/2; i++)
        for (j = 0; j < length; j++)
            dest[npad + 2*i + firstIndex + j] += tempin[npad+i] * coeffs[j];
}