Transformación wavelet de imágenes naturales

A continuación paso a describir la rutina <strong>WIDec</strong> que he programado en MATLAB y se encarga de realizar una descomposición de una de imagen en una wavelet dada.

En realidad el usuario realiza un input de una colección de imágenes, la rutina busca el archivo wtype (que es un listado de nombres de wavelets) y realiza la descomposición de todas y cada una de las imágenes en todas y cada una de las wavelet. Los resultados los almacena en una carpeta con el nombre de la wavelet (carpeta que tiene que estar creada previamente).

El archivo de resultados contiene la imagen original, la imagen reconstruida y todos los coeficientes wavelet.

xA=double(I);% convertimos los valores de la imagen en double

for t=1:nlevel;% nlevel es el número de niveles en que se puede descomponer la imagen, en función de su tamaño (es decir si el tamaño de I es 2^n x 2^n entonces nlevel = n-1)

[A{t},H{t},V{t},D{t}]=dwt2(xA,wtype{nw},'mode','per'); %hacemos la transformación wavelet para cada nivel (wtype{nw}: es una de las wavelet del fichero)

xA=A{t};

end

%%%%% 4 - the image and the details are saved to a file

[iI]=iImage(H,V,D,wtype{nw});%reconstructed image. Realiza una transformada wavelet inversa.

Imágenes sintéticas

Voy a comprobar el funcionamiento de la rutina WIDec, que se encarga de hacer la descomposición en una wavelet definida por el usuario, de una imagen 2ⁿ x 2ⁿ dada.

Para ello, además de revisar el código, aplicaré la rutina a una imagen sintética* y analizaré las diferencias entre sus coeficientes wavelet y los que se utilizaron para la construcción de la imagen sintética.

*evidentemente la wavelet ha de ser la misma tanto para la construcción como para la descomposición.

El resultado es que tanto la imagen original como la reconstruida son iguales. Lo mismo se puede decir de los coeficientes wavelet hasta orden 10^-9 aproximadamente.

Imágenes naturales

La rutina también funciona con imágenes naturales.

He realizado una descomposición de un recorte de 512 x 512 sobre la imagen natural imk00327 de la base de imágenes de Van Hateren con la wavelet db2 y calculado la suma sobre el valor absoluto de la diferencia entre la imagen original y la reconstruida, que da cero.

Para poder comparar las imágenes ha sido necesario normalizar las matrices que las describen, lo que he hecho llevando los valores de cada una de ellas al correspondiente valor de la paleta de grises 0 – 255.

Conclusión

La rutina de descomposición de imágenes WIDec funciona correctamente, lo cual era de esperar teniendo en cuenta que me he apoyado directamente en la biblioteca de funciones de MATLAB.

Generación de imágenes sintéticas

Para nuestro trabajo es conveniente, en determinados momentos, utilizar como banco de pruebas para nuestras rutinas a imágenes construidas artificalmente (sintéticas), cuyas propiedades por diseño conocemos y se ajustan a nuestros fines.

Vamos a construir imágenes de acuerdo a dos procesos, log-normal y log-poisson. Lo esencial del código que he utilizado es el siguiente (MATLAB):

%%%%% Parameters

d = 2;% space dimension

hinf = -1/2 + 1;% most singular exponent

%log-poisson process

k = 1/2;

Dinf = 1;

beta = 1 + hinf/(d - Dinf);

l = (d - Dinf)*(-log(k));

%log-normal process

hmean = 0.0;

sigmah = 1.0;

%hmean-hinf;% we use the log-normal process multifractal spectrum in terms of hinf, D(hinf)=1

mean = hmean*log(k);

sigma = sigmah*(-log(k)/2)^(1/2);

%%%%% 1 - Initialize detail coefficients

%Los coeficientes en el nivel menos resuelto se inicializan (aleatoriamente) en los valores +1 ó -1, excepto el coeficiente A de las aproximaciones cuyo valor es 0.

D{nlevel}(1,1) = (-1)^binornd(1,0.5);

V{nlevel}(1,1) = (-1)^binornd(1,0.5);

H{nlevel}(1,1) = (-1)^binornd(1,0.5);

A{nlevel}(1,1) = 0;

%2 - Generate child coefficients acording to log process

%Ahora propagamos los valores iniciales para cada coeficiente a través de la cascada, del nivel menos resuelto al nivel mas resuelto, utilizando el proceso log-poisson

for t = (nlevel-1):-1:1;

for ti = 1:uint16(2^(nlevel-t));

for tj = 1:uint16(2^(nlevel-t));

D{t}(ti,tj) = 1/2*D{t+1}(ti/2,tj/2)*(k^hinf)*(beta^poissrnd(l))*(-1)^binornd(1,0.5);

H{t}(ti,tj) = 1/2*H{t+1}(ti/2,tj/2)*(k^hinf)*(beta^poissrnd(l))*(-1)^binornd(1,0.5);

V{t}(ti,tj) = 1/2*V{t+1}(ti/2,tj/2)*(k^hinf)*(beta^poissrnd(l))*(-1)^binornd(1,0.5);

end end end

%log normal process

%En el caso log-normal la rutina es idéntica, simplemente empleamos la fórmula correspondiente al proceso

D{t}(ti,tj) = 1/2*D{t+1}(ti/2,tj/2)*lognrnd(mean,sigma)*(-1)^binornd(1,0.5);

H{t}(ti,tj) = 1/2*H{t+1}(ti/2,tj/2)*lognrnd(mean,sigma)*(-1)^binornd(1,0.5);

V{t}(ti,tj) = 1/2*V{t+1}(ti/2,tj/2)*lognrnd(mean,sigma)*(-1)^binornd(1,0.5);

%%%%% 3 - reconstruct the image with the given wavelet

%Para poder reconstruir la imagen simplemente nos faltaría obtener los coeficientes de aproximación A al nivel menos resuelto. Vamos a ir obteniendo cada coeficiente A a cada nivel de resolución realizando una transformada wavelet inversa de los (4) coeficientes del nivel inmediatamente menos resuelto, hasta llegar al nivel más resuelto.

%Es en éste momento cuando decidimos utilizar una wavelet determinada, información que se debe conservar junto con la imagen, puesto que un análisis posterior de la imagen sintética realizada con otra wavelet distinta a aquella con la que se ha construido no será óptima.

for t=(nlevel-1):-1:1;

nl = uint16(2^(nlevel-t-1));

cA = A{t+1}(1:nl,1:nl);

cH = H{t+1}(1:nl,1:nl);

cV = V{t+1}(1:nl,1:nl);

cD = D{t+1}(1:nl,1:nl);

A{t}(1:2*nl,1:2*nl)=idwt2(cA,cH,cV,cD,wtype,'mode','per');

end

%Por último tenemos que:

I = idwt2(A{1},H{1},V{1},D{1},wtype,'mode','per');

%Y si queremos normalizar la imagen (0-255)

imI=uint8(floor((I-min(min(I)))*255/(max(max(I))-min(min(I)))));

Con ésta rutina he obtenido imágenes como éstas:

Imagen sintética de 512 x 512 puntos generada con un proceso log-normal y una wavelet db2:

Y, por otro lado, una imagen análoga utilizando un proceso log-poisson:



Podemos además analizar los histogramas de cada una de las dos imágenes (es decir, cuantos píxeles hay para cada valor de la paleta de grises 0-255), para analizar que distribución tienen y comprobar si se corresponden con las distribuciones que pudieran esperarse de ambos procesos:

[caption id="attachment_79" align="alignright" width="595" caption="log-poisson"][/caption]

[caption id="attachment_78" align="alignleft" width="595" caption="log-normal"][/caption]