Contents

1 Introducción

• La DWT (Discrete Wavelet Transform) es una herramienta matemática que se utiliza, entre otras cosas, para representar señales por niveles de resolución. Así, por ejemplo, si calculamos la DWT de una imagen obtendremos un conjunto de coeficientes DWT que descomprimidos por subbandas sucesivas permitirán aumentar progresivamente la resolución espacial de la imagen descomprimida.
• También es posible descomprimir una imagen usando, por ejemplo, siempre el máximo nivel de resolución, y en este caso, lo que ocurre durante la recepción progresiva de la imagen es que la calidad de la misma aumenta conforme se reciben más y más datos.
• A la primera forma de descompresión se le conoce como escalabilidad espacial y a la segunda, como escalabilidad por calidad.

2 Procedimiento

• En esta práctica vamos a comprimir una imagen quantificando los coeficientes wavelet de la imagen y luego usando un codificador entrópico estándar (gzip) para mapear los coeficientes cuantificados en un código de longitud variable. En otras palabras, el compresor sería:
  1. Leer una imagen desde disco.
  2. Calcular la DWT de la imagen.
  3. Quantificar los coeficientes DWT. Recuérdese que la ganancia de los coeficientes wavelet no es la misma en todas las subbandas. Esto significa que a la hora de cuantificarlos, no deberíamos usar el mismo nivel de cuantificación sobre coeficientes con diferentes ganancias (que pertenecen a diferentes subbandas).
  4. Escribir en el disco los coeficientes cuantificados (no es necesario mapear los coeficientes cuantificados sobre un código más eficiente desde el punto de vista del gasto de memoria porque la siguiente etapa se va a encargar jústamente de ésto).
  5. Comprimir los coeficientes con gzip.

  Y su correspondiente descompresor:

  1. Descomprimir los coeficientes (con gzip).
  2. Leer desde el disco la matriz de coeficientes cuantificados.
  3. Recuperar los valores aproximados de los coeficientes. Recuérdese que han sido cuantificados, teniendo en cuenta sus ganancias.
  4. Calcular la DWT inversa.
  5. Escribir la imagen a disco.

• Por comodidad, trabajaremos con imágenes PPM (Portable PixMap) que son fáciles de tratar. Sin embargo, los coeficientes se escriben en un formato raw (crudo), con una cabecera sencilla (creada por nosotros) que permita descomprimir la imagen sin necesidad de especificar los parámetros de descompresión. Dicha cabecera deberá especificar: (1) las dimensiones (alto y ancho, en este orden) de la imagen y (2) el número de niveles de la DWT.

3 Estructura de una imagen PPM binaria

• El compresor leerá imágenes PPM y escribirá una matriz de coeficientes wavelet cuantificados usando un formato a medida (el descrito en la sección anterior). El descompresor leerá dicho formato y escribirá imágenes PPM.
• Una imagen binaria cruda RGB (Red Green Blue) encapsulada según el formato PPM tiene la siguiente estructura:
  /* Note: uintN = N bits unsigned integer */  
   
  struct PPM {  
    uint8[3] magic_number = "P6\n";  
    uint8* width;                   /* In ASCII!!! */  
    uint space = " ";  
    uint8* height;                  /* In ASCII!!! */  
    uint8 new_line_1 = "\n";  
    uint8* maximum_component_value; /* 255 or 65535. Again, in ASCII!!!  
                                    (Components values are always positive). */  
    uint8 new_line_2 = "\n";  
    if (maximum_component_Value == 255) {  
      RGB8[width][height] matrix;  
    } else {  
      RGB16[Width][Height] matrix;  
    }  
  };  
   
  struct RGB8 {  
    uint8 R, G, B;  
  };  
   
  struct RGB16 {  
    uint16 R, G, B;  
  };

4 La transformada wavelet 1D Lineal (5/3)

• A continuación se presenta el pseudo-código de la DWT unidimensional, directa e inversa. Existen 4 métodos porque las secuencias de muestras (pixels en el caso de las imágenes) podrían tener un número par o impar de elementos, y en cada caso, el cálculo de la DWT depende de esta circunstancia.
• Nótese que el tipo de dato usado para representar la señal de entrada signal y las subbandas de frecuencia low y high tienen 32 bits de precisión (con signo) y las imágenes tienen sin embargo sólo 8 bits de profundidad. Esto deberá tenerse en cuenta cuando se lean y escriban las imágenes desde y a disco.
  53_analyze_even(int32* signal, int32* low, int32* high, uint32 samples) {  
   
    for (uint32 i=0; i < (samples/2-1); i++) {  
      uint32 i2 = i*2;  
      high[i] = signal[i2+1] - (signal[i2] + signal[i2+2])/2;  
    }  
    high[i] = signal[samples-1] - signal[samples-2];  
   
    low[0] = signal[0] + high[0]/2;  
    for (uint32 i=1; i < (samples/2); i++) {  
      uint32 i2 = i*2;  
      low[i] = signal[i2] + (high[i] + high[i-1])/4;  
    }  
   
  }  
   
  53_synthesize_even(int32* signal, int32* low, int32* high, uint32 samples) {  
   
    signal[0] = low[0] - h[0]/2;  
    for (unit32 i=1; i < (samples/2) ;i++) {  
      int i2 = i*2;  
      signal[i2] = low[i] - (high[i] + high[i-1])/4;  
    }  
   
    for(uint32 i=0; i < (samples/2-1); i++) {  
      uint32 i2 = i*2;  
      signal[i2+1] = high[i] + (signal[i2] + signal[i2+2])/2;  
    }  
    signal[samples-1] = high[i] + signal[samples-2];  
   
  }  
   
  53_analyze_odd(int32* signal, int32* low, int32* high, uint32 samples) {  
   
    for(uint32 i=0; i < (samples/2); i++) {  
      uint32 i2 = i*2;  
      high[i] = signal[i2+1] - (signal[i2] + signal[i2+2])/2;  
    }  
   
    low[0] = signal[0] + high[0]/2;  
    for(uint32 i=1; i < (samples/2) ;i++) {  
      uint32 i2 = i*2;  
      low[i] = signal[i2] + (high[i]+high[i-1])/4;  
    }  
    low[i] = signal[samples-1] + high[i-1]/2;  
   
  }  
   
  synthesize_odd(int32* signal, int32* low, int32* high, uint32 samples) {  
   
    signal[0] = low[0] - high[0]/2;  
    for(uint32 i=1; i < (samples/2); i++) {  
      uint32 i2 = i*2;  
      signal[i2] = low[i] - (high[i] + high[i-1])/4;  
    }  
    signal[samples-1] = low[i] - high[i-1]/2;  
   
    for(uint32 i=0; i < (samples/2); i++) {  
      uint32 i2 = i*2;  
      signal[i2+1] = high[i] + (signal[i2] + signal[i2+2])/2;  
    }  
   
  }

5 La transformada wavelet 2D

• Es posible calcular la transformada 2D a partir de la 1D usando la propiedad de la separabilidad de la DWT. En concreto, la DWT 2D se calcula aplicanto la DWT 1D primero a todas las filas de la imagen y a continuación a todas las columnas. El orden (primero filas o primero columnas) es irrelevante, siempre y cuando al descomprimir se realize en el orden inverso.
   
  /* Two "global" vectors of unsigned integers. */  
  uint32* in_line;  
  uint32* out_line;  
   
  /* Asign memory to those vectors taking into account the maximun  
     row or column size. */  
  line_allocation(uint32 width, uint32 height) {  
   
    uint32 max_length = 0;  
    if (max_length < width) {  
       max_length = width;  
    } else if (max_length < height) {  
       max_length = heidht;  
    }  
    in_line = memory_allocation(uint32, max_length);  
    out_line = memory_allocation(uint32, max_length);  
   
  }  
   
  /*  
    Structure in the memory of the matrix of pixels (that is, a vector of  
    pointers that points to vectors of pixels)):  
   
                        +----> cols (columns)  
    image               |  
    |                   v  
    v                  rows  
    +---------------+    +-------------+-------------+     +------------------+  
    |      image[0] |--->| image[0][0] | image[0][1] | ... | image[0][cols-1] |  
    +---------------+    +-------------+-------------+     +------------------+  
    |      image[1] |---> ...  
    +---------------+  
    :               :  
    +---------------+    +------------------+     +-----------------------+  
    | image[rows-1] |--->| image[rows-1][0] | ... | image[rows-1][cols-1] |  
    +---------------+    +------------------+     +-----------------------+  
  */  
   
  analyze_2d(int32** image, uint32 rows, uint32 cols, uint8 levels) {  
   
    for(uint8 l=0; l<levels; l++) {  
   
      uint32 x = cols;  
      cols /= 2;  
      if (cols == 0) cols = 0;  
   
      uint32 y = rows;  
      rows /= 2;  
      if (rows == 0) rows = 0;  
   
      /* Iterate over the columns */  
      if (x & 1) { /* There is a odd number of columns */  
        for(uint32 j = 0; j < y; j++) {  
          for(uint32 i = 0; i < x; i++) in_line[i] = image[j][i];  
          analyze_1d_odd(in_line, image[j], image[j] + width + 1, x);  
        }  
      } else { /* There is a even number of columns */  
        for(uint32 j = 0; j < y; j++) {  
          for(uint32 i = 0; i < x; i++) in_line[i] = image[j][i];  
          analyze_1d_even(in_line, image[j], image[j] + width, x);  
        }  
      }  
   
      /* Iterate over the rows */  
      if (y & 1) { /* There is a odd number of rows */  
        for(uint32 i = 0; i < x; i++) {  
          for(uint32 j = 0; j < y; j++) in_line[j] = image[j][i];  
          analyze_1d_odd(in_line, out_line, out_line+height+1, y);  
          for(uint32 j = 0; j < y; j++) image[j][i] = out_line[j];  
        }  
      } else { /* There is a even number of rows */  
        for(uint32 i=0; i<x; i++) {  
          for(uint32 j=0;j <y; j++) in_line[j]=image[j][i];  
          analyze_1d_even(in_line, out_line, out_line + height, y);  
          for(uint32 j=0;j <y; j++) image[j][i]=out_line[j];  
        }  
      }  
    }  
  }  
   
  synthesize_2d(int32 **image, uint32 width, uint32 height, uint8 levels) {  
    int nx, ny, lv;  
    int i,j;  
   
    uint32 x = width >> levels;  
    uint32 y = height >> levels;  
   
    for (uint8 lv = levels-1; lv >= 0; lv--) {  
   
      uint32 mx = x;  
      x = width >> lv;  
      if(x==0) x = 1;  
   
      uint32 my = y;  
      y = height >> lv;  
      if(y==0) y = 1;  
   
      /* Columns transformation */  
      if(y & 1) { /* Odd number of rows */  
        for(uint32 i = 0; i < x; i++) {  
          for(uint32 j = 0; j < y; j++) in_line[j] = image[j][i];  
          synthesize_1d_odd(out_line, in_line, in_line + my + 1, y);  
          for(uint32 j = 0; j < y; j++) image[j][i] = out_line[j];  
        }  
      } else { /* Even number of rows */  
        for(uint32 i = 0; i < x; i++) {  
          for(uint32 j = 0; j < y; j++) in_line[j] = image[j][i];  
          synthesize_1d_even(out_line, in_line, in_line + my, y);  
          for(uint32 j = 0; j < y; j++) image[j][i] = out_line[j];  
        }  
      }  
   
      /* Rows transformation */  
      if(x & 1) { /* Odd number of cols */  
        for(uint32 j = 0; j < y; j++) {  
          memcpy(in_line, image[j], x*sizeof(int32));  
          synthesize_1d_odd(image[j], in_line, in_line + mx + 1, x);  
        }  
      } else { /* N’umero par de columas */  
        for(j=0;j<ny;j++) {  
          memcpy(in_line, image[j], x*sizeof(int32));  
          synthesize_1d_even(image[j], in_line, in_line + mx, x);  
        }  
      }  
    }  
  }  
  

6 Entregables

• Deberá implementarse un compresor como el descrito, junto a su correspondiente descompresor. En la memoria aparecerá el código desarrollado y comentado, así como una “demo” de su uso.
• La nota de la práctica dependerá de la exposición (en la memoria) del trabajo realizado así como del ratio bit-rate/distorsión introducido por el compresor. Nótese que dicho ratio dependerá del paso de cuantificación elegido (a mayor paso, mayor distorsión aunque mayor tasa de compresión).