wavelet transform in C

You can find something here:
http://www.c3.lanl.gov/feemads-sw/doxygen/impl/html/classFEx_1_1DWT.html

Also:
http://www.ddj.com/ftp/1992/1992.04

Hi
 Thanks for your information but that is too hard to be understand.  I just need to transform an image with irregular size.

for a quick fix just scale the image to square first. It will do well enough... better if you scale up smaller if you scale down

Hi,
  The following is the sample code for this purpose but without the main function.

  How can I add the main function in order to make it runable??

  The image size is 720x576 / 512x512 / 3027 x 2048 / 1465 x 1999

  This is my first times to touch this kind of program and how can we handle the large image size (e.g. 3027x2048??)

so I hope finally
input image -> Transform -> transformed output

Hope someone can help.
*********************************************************

/* describes an image by its memory location and its organization */
struct image {
  long *pic;
  long modulo;
  long w,h;  /* dimensions */
};

/* load a PGM encoded image into a long-array, return width, height
 * and row-offset suitable for 5/3 filtering of "cascade" decomposition
 * levels
 */
long *loadpgm(FILE *in,long *w,long *h,long *modulo,int blocksize,int cascade)
{
  long d;
  int ascii;
  long x,y;
  long *p,*pic;
  long *mem;
  char buf[256];
  long rh;
 
  fscanf(in,"%s\n",buf);
  if (!strcmp(buf,"P2")) {
    ascii = 1;
  } else if (!strcmp(buf,"P5")) {
    ascii = 0;
  } else {
    fprintf(stderr,"stdin is not a valid PGM image.\n");
    exit(10);
  }
 
  do {
    fgets(buf,256,in);
  } while(buf[0] == '#'); /* skip comments */

  // read width, height.
  *w = 0, *h = 0;
  sscanf(buf,"%ld %ld",w,h);
  if (w <= 0 || h <= 0) {
    fprintf(stderr,"input image dimensions are invalid.\n");
    exit(10);
  }

  do {
    fgets(buf,256,in);
  } while(buf[0] == '#'); /* skip comments */

  d = 0;
  sscanf(buf,"%ld",&d);
  if (d != 255) {
    fprintf(stderr,"this program handles only bit precisions of 255\n");
    exit(10);
  }

  /* Get enough memory to handle this */
  *modulo  = *w + 2*(blocksize + cascade * 3);
  rh       = *h + 2*(blocksize + cascade * 3);
  mem      = malloc(*modulo * rh * sizeof(long));
  memset(mem,0,*modulo * rh * sizeof(long));
  pic      = mem +  (blocksize + cascade * 3) * (*modulo + 1);
  p        = pic;

  fprintf(stderr,"loading the image...\n");

  for(y=0;y<*h;y++) {
    for(x=0;x<*w;x++) {
      long v;
      /* Now read in one line after another */
      if (ascii) {
fscanf(in,"%ld ",&v);
      } else {
v = fgetc(in);
      }
      if (v < 0) {
fprintf(stderr,"unexpected end of file\n");
exit(10);
      }
      *p++   = v - 128;
    }
    p += *modulo - *w;
  }

  return pic;
}

/* Perform 5/3 filtering, building an array of images from the
 * image data passed in.
 * NOTE: This throws intermediate bands away. The modifications
 * to keep these should be rather simple.
 */
struct image *filterimage(long *pic,long w,long h,long modulo,
  long cascade,long blocksize)
{
  long x,y;
  long *p,*q;
  struct image *picmap;
 
  /* Generate an array for the four images we get from filtering:
   * LL,HL,LH,HH from straight wavelet filtering
   */
  picmap = calloc(4,sizeof(struct image));

  fprintf(stderr,"Running the wavelet transform....\n");
  while(cascade > 0) {
    cascade--;  /* Down in the next lower cascade now */
    waveletfilter(pic,1,modulo,w,h);
    waveletfilter(pic,modulo,1,h,w);
    /* Sort the output such that we have direct addressing of the
     * high-high pass or whichever band we want. This is only
     * necessary for all but the last cascade as we now build
     * the statistics for all bands at once.
     */    
    if (cascade == 0) {
      long outmod,outh,band;
      /* Get dimensions of filtered image now */
      outmod  = (w + (blocksize << 2))>>1;
      outh    = (h + (blocksize << 2))>>1;
 
      for(band = 0;band<4;band++) {
picmap[band].pic    = calloc(outmod * outh,sizeof(long));
picmap[band].pic   += (blocksize * outmod) + blocksize;
picmap[band].modulo = outmod;
picmap[band].w      = w>>1;
picmap[band].h      = h>>1;
hvextract(pic,modulo,picmap[band].pic,outmod,w,h,band);
      }
      return picmap;
    }
    /* Otherwise get only the low-low-pass here. */
    /* otherwise, we would need to extract other bands in here
     * as well.
     */
    hvextract(pic,modulo,pic,modulo,w,h,0);
    h >>= 1;
    w >>= 1;
  }
}

/* Write an image (in longs) out to a pgm file
 */
void savepgm(FILE *out,long *img,long imgmod,long w,long h,long min,long max)
{
  long x,y,v;

  fprintf(out,"P5\n%ld %ld\n255\n",w,h);
  for(y = 0;y<h;y++) {
    for (x = 0;x<w;x++) {
      v = (img[x+y*imgmod] - min) * 255 / (max - min);
      if (v < 0)   v = 0;
      if (v > 255) v = 255;
      putc(v,out);
    }
  }
}

/* Make a simplified wavelet transform with the 5/3 filter
 */
void waveletfilter(long *o,int incx,int incy,int lx,int ly)
{
  int i,j;
  long *row,*p;

  row = o;
  for(i=0;i<ly;i++) {
    p = row - incx;
    for(j =-1;j<(lx + 1)>>1;j++) {
      p[0]    -= (p[-incx] + p[incx])>>1;
      p       += incx<<1;
    }
    p = row;
    for(j = 0;j<(lx + 1)>>1;j++) {
      p[0] += (p[-incx] + p[incx] + 2)>>2;
      p    += incx<<1;
    }
    row += incy;
  }
}

Hi,
 I added the main program as belows but I am not familiar with the PGM header.  Who can help to modify the following to be workable??

***************************************************
#include <stdio.h>
#include <stdlib.h>
/* describes an image by its memory location and its organization */
struct image {
  long *pic;
  long modulo;
  long w,h;  /* dimensions */
};

/* load a PGM encoded image into a long-array, return width, height
 * and row-offset suitable for 5/3 filtering of "cascade" decomposition
 * levels
 */
long *loadpgm(FILE *in,long *w,long *h,long *modulo,int blocksize,int cascade)
{
  long d;
  int ascii;
  long x,y;
  long *p,*pic;
  long *mem;
  char buf[256];
  long rh;
 
  fscanf(in,"%s\n",buf);
  if (!strcmp(buf,"P2")) {
    ascii = 1;
  } else if (!strcmp(buf,"P5")) {
    ascii = 0;
  } else {
    fprintf(stderr,"stdin is not a valid PGM image.\n");
    exit(10);
  }
 
  do {
    fgets(buf,256,in);
  } while(buf[0] == '#'); /* skip comments */

  // read width, height.
  *w = 0, *h = 0;
  sscanf(buf,"%ld %ld",w,h);
  if (w <= 0 || h <= 0) {
    fprintf(stderr,"input image dimensions are invalid.\n");
    exit(10);
  }

  do {
    fgets(buf,256,in);
  } while(buf[0] == '#'); /* skip comments */

  d = 0;
  sscanf(buf,"%ld",&d);
  if (d != 255) {
    fprintf(stderr,"this program handles only bit precisions of 255\n");
    exit(10);
  }

  /* Get enough memory to handle this */
  *modulo  = *w + 2*(blocksize + cascade * 3);
  rh       = *h + 2*(blocksize + cascade * 3);
  mem      = malloc(*modulo * rh * sizeof(long));
  memset(mem,0,*modulo * rh * sizeof(long));
  pic      = mem +  (blocksize + cascade * 3) * (*modulo + 1);
  p        = pic;

  fprintf(stderr,"loading the image...\n");

  for(y=0;y<*h;y++) {
    for(x=0;x<*w;x++) {
      long v;
      /* Now read in one line after another */
      if (ascii) {
fscanf(in,"%ld ",&v);
      } else {
v = fgetc(in);
      }
      if (v < 0) {
fprintf(stderr,"unexpected end of file\n");
exit(10);
      }
      *p++   = v - 128;
    }
    p += *modulo - *w;
  }

  return pic;
}

/* Perform 5/3 filtering, building an array of images from the
 * image data passed in.
 * NOTE: This throws intermediate bands away. The modifications
 * to keep these should be rather simple.
 */
struct image *filterimage(long *pic,long w,long h,long modulo,
  long cascade,long blocksize)
{
  long x,y;
  long *p,*q;
  struct image *picmap;
 
  /* Generate an array for the four images we get from filtering:
   * LL,HL,LH,HH from straight wavelet filtering
   */
  picmap = calloc(4,sizeof(struct image));

  fprintf(stderr,"Running the wavelet transform....\n");
  while(cascade > 0) {
    cascade--;  /* Down in the next lower cascade now */
    waveletfilter(pic,1,modulo,w,h);
    waveletfilter(pic,modulo,1,h,w);
    /* Sort the output such that we have direct addressing of the
     * high-high pass or whichever band we want. This is only
     * necessary for all but the last cascade as we now build
     * the statistics for all bands at once.
     */    
    if (cascade == 0) {
      long outmod,outh,band;
      /* Get dimensions of filtered image now */
      outmod  = (w + (blocksize << 2))>>1;
      outh    = (h + (blocksize << 2))>>1;
 
      for(band = 0;band<4;band++) {
picmap[band].pic    = calloc(outmod * outh,sizeof(long));
picmap[band].pic   += (blocksize * outmod) + blocksize;
picmap[band].modulo = outmod;
picmap[band].w      = w>>1;
picmap[band].h      = h>>1;
hvextract(pic,modulo,picmap[band].pic,outmod,w,h,band);
      }
      return picmap;
    }
    /* Otherwise get only the low-low-pass here. */
    /* otherwise, we would need to extract other bands in here
     * as well.
     */
    hvextract(pic,modulo,pic,modulo,w,h,0);
    h >>= 1;
    w >>= 1;
  }
}

/* Write an image (in longs) out to a pgm file
 */
void savepgm(FILE *out,long *img,long imgmod,long w,long h,long min,long max)
{
  long x,y,v;

  fprintf(out,"P5\n%ld %ld\n255\n",w,h);
  for(y = 0;y<h;y++) {
    for (x = 0;x<w;x++) {
      v = (img[x+y*imgmod] - min) * 255 / (max - min);
      if (v < 0)   v = 0;
      if (v > 255) v = 255;
      putc(v,out);
    }
  }
}

/* Make a simplified wavelet transform with the 5/3 filter
 */
void waveletfilter(long *o,int incx,int incy,int lx,int ly)
{
  int i,j;
  long *row,*p;

  row = o;
  for(i=0;i<ly;i++) {
    p = row - incx;
    for(j =-1;j<(lx + 1)>>1;j++) {
      p[0]    -= (p[-incx] + p[incx])>>1;
      p       += incx<<1;
    }
    p = row;
    for(j = 0;j<(lx + 1)>>1;j++) {
      p[0] += (p[-incx] + p[incx] + 2)>>2;
      p    += incx<<1;
    }
    row += incy;
  }
}


int main(int argc, char **argv)
{
    int i, j, k, cl, bits;
    double rate;
    IMAGE *img;
    ENCODER *enc;
    clock_t st_clk, ed_clk;
    char *infile, *outfile;
    FILE *fp;

    st_clk = clock();
    setbuf(stdout, 0);
    infile = outfile = NULL;
     if (infile == NULL) {
            infile = argv[i];
          } else {
            outfile = argv[i];
          }
      }
    }
    if (infile == NULL || outfile == NULL) {
      printf("infile:     Input file (must be in a raw PGM format)\n");
      printf("outfile:    Output file\n");
      exit(0);
    }
    img = loadpgm(infile);
    fp = fileopen(outfile, "wb");
    printf("%s -> %s (%dx%d)\n", infile, outfile, img->width, img->height);
}

OK, so what is DWT and PGM and some more background please. I think I will be able to help in C language part.

OK let me gather something and explain a bit here.

http://vsp.ee.nthu.edu.tw/EE5632/handout/ch10.ppt
have a look of this first (look at the 5/3 reversible transform as it is lossless)
I will show you a simple example later

not quite understand. what I understood is if sum of 2 integer doesn't overflow, you can save 1 bit.

no problem
I will give u an example, it is easy
don't worry

Original image:

50  51  54  56  < row1
53  52  50  57
45  52  52  49
43  48  47  41

Let takes the first row (row 1) as a an example

2-D filtering is just the tensor product of horizontal  and vertical filtering.
The order wouldn't matter if we'd had no rounding step, but as we do, the
order must be given. JPEG2000 expects for forward transformation first
vertical, then horizontal filtering.

1) First step: Extension of the data. This is what JPEG2000 Part I calls
the 1D_EXTD procedure. The 5/3 filter requires an overlap of at most two, so
to be on the safe side, we need to generate two additional pixels on each
edge of the data. Less would be sufficient, see table F-8 in the JPEG2000 specs,
or see my comment below:

54 51|50 51 54 56|54 51

Note that I just flipped the data over by using 50 and 56 as mirror points.
This is called WSS extension of the data. ("Whole Sample Symmetric").

2) The filtering process applies. For that, we need to know which pixels
are even, and which are odd.  For simplicity, let's assume
that the top left "50" is at the canvas origin (0,0) and is hence an
even pixel. (canvas offset = Xorg,Yorg = 0,0):

2.a) First lifting step: We need to process all odd pixels.  

a(2n+1) + floor((a(2n) + a(2n+2))/2) -> a(2n+1)

Hence: n =  0: 51 - (50 + 54)/2 = -1
       n =  1: 56 - (54 + 54)/2 = -2

We hence get:

54 -1|50 -1 54 -2|54 51

2.b) Perform the second lifting step
 
a(2n) + floor((a(2n-1) + a(2n+1) + 2)/4) -> a(2n)
 
Hence, in our example, n=0 and n=1.

n=0: 50 + ((-1 + -1 + 2)/4) = 50
n=1:    54 + ((-1 + -2 + 2)/4) = 53

The result is:

54 -1|50 -1 53 -1|54 51

Note that we didn't make use of the rightmost 51, so it was not necessary
to construct it. The result of the filtering process is now between the borders | and |, with the
even coefficients being the low-passes and the odd coefficients being
the high-passes, hence:

50 53 is the low pass and
-1 -1 is the high pass.


This process has now to be applied first vertically on columns, then horizontally.

ok Do it vertically:

**************************************
53  52  53  59
6   1   -3  5
46  52  51  49
-2  -4  -4  -8

then Horizontal.  (use the result from vertical transform)

53  -1  54  6
7   0   -1  8
48  4   53  -2
-2  -1  -4  0


reconstruct:

53  54  -1  6
48  53  -1  8
7   4   0  -2
-2  -1  -4  0


You can see how does the pic look like here:

http://www.gvsu.edu/math/wavelets/student_work/EF/how-works.html

 



 

>> 56 - (54 + 54)/2 = -2
isn't it should be 2 instead of -2. but anyway, I understand the concept, just to confirm no special meaning.

>> reconstruct:
what is the rule of reconstruction?

your main program have problem. big problem :-)

can you point me a test file, test image file. I'm working on the C code, will continue tomorrow.

I think I'll be able to give you at least a test C code tomorrow.

do u have email?? I can send u some images

send to my hotmail, the user id is the same as my id of EE. I don't type it out because somebody use tools to scan page to get email address.

I sent you a 512x512 image for testing
and if you want to try a big one, you can download it from

http://www.gorex.com.hk/EE/bigimage.zip

wait..
I am better send you a small one first

ok sent a small image size with 256 x 256

I have send the code to you via email, but the function "hvextract" which is used by "filterimage" function is missing from you original code. So the code cannot link to executable because of that.

and what is blocksize, cascade, min and max means?

so can I test it?

In fact, u can ignore all those variables.  

>> so can I test it?
no

>>In fact, u can ignore all those variables.
those variables are required to called the loadpgm and savepgm function.

and:

Comment from kennethxu  12/18/2002 09:41PM PST  
>> 56 - (54 + 54)/2 = -2
isn't it should be 2 instead of -2. but anyway, I understand the concept, just to confirm no special meaning.

>> reconstruct:
what is the rule of reconstruction?
 

ok can you just ignore the PGM stuff?  As it is not necessary (we can use RAW image) and it is much easier to implement.

Reconstruction rule is up to u.  But in my mind, you can put all the "big" number in the upper left region.

in this case, majority code will have to be re-write. I'm trying

is it complicate?

here is the code that can read the row image with different size and maybe u can just add the 5/3 wavelet inside the code.  All the input data will be put in image[]:
****************************************************
#include <stdio.h>
#include <stdlib.h>


main()
{
  FILE *infile,*outfile;
  char in[15], out[15];
 unsigned char **image, **res_image;
 int  col_size,row_size;

 /* Allocate a char* array to hold the rows */
   
   printf("Enter an input file\n");
   scanf("%s", in);
   printf("Enter an output file\n");
   scanf("%s", out);
   printf("Enter Col size\n");
   scanf("%d", &col_size);
    printf("Enter row size\n");
   scanf("%d", &row_size);
   
   

   image  = (unsigned char **) calloc(row_size,sizeof(char *));
   for(row=0; row<row_size; row++)
   {
    image[row]  = (unsigned char *) calloc(col_size,sizeof(char));
   }

   res_image  = (unsigned char **) calloc(row_size,sizeof(char *));
   for(row=0; row<row_size; row++)
   {
    res_image[row]  = (unsigned char *) calloc(col_size,sizeof(char));
   }


 /* Allocate a buffer for each individual row */
 //for(r = 0; r < row_size; r++) {
 //     image[r] = malloc(col_size);
        //res_image[r] = malloc(row_size);
 //}


  infile=fopen(in,"rb");
  outfile=fopen(out,"wb");

  if(infile==NULL)
  {
  printf("Unable to open file %s",in);
  exit(1);
  }
  for(r=0; r<row_size; r++)
  fread(image[r],col_size,sizeof(char),infile);


---->>  wavelet transform .......


  for(row=0; row<row_size; row++)  // write the filtered data out
   for(col=0; col<col_size; col++)
     fputc(res_image[row][col], outfile);
 fclose(infile);

 /* Free allocated memory */
 for(r = 0; r < row_size; r++) {
     free(image[r]);
     free(res_image[r]);
 }
 free(image);
 free(res_image);
}

tested with your imaged:

wavelet lena.raw lena2.raw 256 256

Hi,
  I tried it!  But the output is wrong.  The output image shuld look like four sub-image in 4 regions.  Refer the following and you will know how it looks like.

http://www.gvsu.edu/math/wavelets/student_work/EF/how-works.html

(Figure 7)


Also, is it possible to let user assign the levels??  as we can use different levels to transform.

In figure 7 above link, the left pic is the original, middle is one level transform, and the last one is 2-level transform.

Not hurry, just do it when you have time
Merry Xmas :o)

I coded as what you have described, it does exactly the step you said it sould.

I don't know what do you mean by level and etc. My code is very easy to extend, so you can modify it if you know more detail about transform algorithm.

Ok thanks.  I do it myself.

>> The output image shuld look like four sub-image in 4 regions.

I guess this might because the reconstruct rule is not correct, you asked me to move big numbers to up-left. in fact those small number also matters.

let me know if you need help to understand my code.

hm...yes maybe the reconstruction rule got problem.  Your code is fine and I can understand.