\documentclass[11pt,leqno]{article}
\usepackage{amssymb, amsmath} 
\usepackage{latexsym} 
\usepackage{graphics} 
\usepackage{url} 
\usepackage{color} 
\usepackage{hyperref} 

%\usepackage[]{showkeys}
%\usepackage[notcite,notref,color]{showkeys}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	SPANISH
%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[spanish]{babel}
\usepackage[latin1]{inputenc}
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	FORMAT
%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setlength{\textheight}{23.0cm}
\setlength{\textwidth}{15.75cm}
\setlength{\oddsidemargin}{0.5cm}
\setlength{\evensidemargin}{0.5cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	LISTINGS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}

\definecolor{dgreencolor}{rgb}{0.2,0.4,0.0}

\newcommand{\code}[1]{\text{\small\tt #1}}

\lstset{
	aboveskip = 12pt,
	belowskip = 10pt,
	basicstyle = {\fontsize{7pt}{7pt}\rmfamily},
	commentstyle = {\ttfamily\color{dgreencolor}},
	stringstyle = {\ttfamily\color{red}\bfseries},
	keywordstyle = {\ttfamily\color{blue}\bfseries},
	inputencoding = latin1,	
	showstringspaces = true, %only spaces in strings
	numbers = left,
	numberstyle = {\fontsize{7pt}{7pt}\ttfamily},
	numbersep = 8pt, %separation numbers - listing
	breaklines = true,
	breakatwhitespace = true, %breaks only at white spaces
	postbreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\color{red}\hookrightarrow\space}},
	breakindent = 40pt, %indentation after line break
	fontadjust = yes,
	xleftmargin = 2.5cm
}



%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	COUNTERS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcounter{ejeno}
\setcounter{ejeno}{1}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	MACROS
%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{{\mathbb R}} 
\newcommand{\C}{{\mathbb C}} 
\newcommand{\N}{{\mathbb N}} 
\newcommand{\Z}{{\mathbb Z}} 
\newcommand{\Q}{{\mathbb Q}} 

\newcommand{\pvec}[1]{\vec{#1}\mkern2mu\vphantom{#1}'}

\newcommand{\ejer}{\ \par\vspace{-7pt}{\bf\number\value{ejeno})\ }\addtocounter{ejeno}{1}}



\begin{document}
\pagestyle{myheadings}
\markboth{}{\rm\textsf{\jobname}\qquad\hfill}

Antes de nada, como te dije, cambiamos el punto 6 del temario planeado a uno relacionado. El título podría ser ``Algunos métodos de restauración de imágenes'' u otra variante que te guste. El punto 7 lo omitimos, por lo que esta hoja es la última.

Necesitarás usar \code{matlab}. Supongo que tienes acceso {\sl online}  porque en principio todos los alumnos de la UAM lo tenéis. No te asustes si no recuerdas demasiado de \code{matlab} porque los programas principales te los daré yo y tus tareas serán fundamentalmente experimentar con ellos ajustando los parámetros.  También puedes usar \code{octave}, si lo prefieres, en ese caso deberás cargar los paquetes \code{image} y \code{signal}, incluyendo al principio de los programas de esta hoja \code{pkg load image} y \code{pkg load signal}.

\

Para comenzar debes aprender las bases de lo que se llama análisis de Fourier discreto


\ejer
Lee \S2.2.1 \cite{chamizo_sp} hasta el final de la página 63. Conviene que primero des un vistazo a la parte de \S1.2.1 en la  página 10 donde se explica la notación $e(x)$ y otras cosas. 
\medskip

En \code{matlab} la DCT se calcula con \code{dct} y su inversa con \code{idct}. Los análogos bidimensionales que aparecerán en los programas son  \code{dct2} e \code{idct2}.

Tomemos una función $1$-periódica que no sea una combinación lineal de senos y cosenos, por ejemplo
\[
 f(x)=\sen\Big({\frac 12 \sen(2\pi x)-\frac 14 \cos(6\pi x)}\Big).
\]
Como $f$ es suave, las frecuencias altas deberían tener poca importancia al analizar~$f$. 
El siguiente programa, discretiza $f(x)$ en la variable \code{y} usando \code{N} nodos, omite los últimos~\code{N-M} valores de su DCT (los de las frecuencias altas) e invierte el resultado para obtener \code{yrec}. Si dibujas las gráficas de \code{y} e \code{yrec}, por ejemplo con \code{plot(x,y,x,yrec)}, verás que se parecen mucho a pesar de que nos hemos olvidado de casi el $90\%$ de las frecuencias (esta es la clave para algunos métodos de compresión, aunque se suelen usar \emph{wavelets} en vez de cosenos \cite{walker}).
\begin{lstlisting}[language=Matlab]
N = 128
M = 15
x = linspace(0,1,N);
y = sin(sin(2*pi*x)/2-cos(6*pi*x)/4);

mmask = zeros(1,N);
mmask(1:M) = 1;
yrec = idct( dct(y).*mmask );
\end{lstlisting}

\ejer
Haz una tabla o una gráfica que muestre el tamaño de $\|\code{y}-\code{yrec}\|_\infty$ cuando \code{M} vale $10$, $20$,{\dots}, $100$. 

\ejer
Utiliza a lo más una página y media para: 1) Describir las fórmulas para la DFT, la DCT y sus inversas; 2) dar una prueba de la fórmula de inversión para la DFT (con matrices unitarias o como tú prefieras); 3) incluir el ejercicio anterior para ilustrar que las funciones suaves tienen DCT (y DFT) que decae rápidamente.  

\

Al igual que la transformada de Fourier habitual, la DFT y la DCT se extienden a más variables. Simplemente se aplican las fórmulas a cada variable por separado. Una imagen en blanco y negro (tonos de gris) es una función que al pixel $(i,j)$  le asigna un número que indica su tono de gris. 
Este número normalmente está entre~$0$ y~$255$ aunque en los programas \code{matlab} posteriores se escala entre~$0$ y~$1$ con la instrucción \code{im2double}, donde~$0$ es negro y~$1$ es blanco\footnote{El caso de imágenes en color no lo consideraremos pero es similar con tres números que corresponden a los \emph{canales de color} R ({\sl red}), G ({\sl green}) y  B ({\sl blue}).}. 
Entonces tiene sentido considerar la DCT de una foto en blanco y negro. Habitualmente habrá discontinuidades correspondientes a aristas o sombras y por tanto no es de esperar que haya un decaimiento drástico como en el ejemplo anterior. Como un comentario al margen, en el formato JPEG \cite[\S2.2.2]{chamizo_sp} \cite{salomon} se hace una división en trozos de $8\times 8$ píxeles para que en la mayoría de ellos sea fácil que la porción de foto sea suave. 
\smallskip

En el resto de la hoja experimentarás con tres métodos para restaurar ciertos defectos de imágenes y dos de ellos tienen que ver con el análisis de Fourier discreto.

\paragraph{1. Filtrado homomorfo.}
Vamos a ver cómo usar la DCT para homogeneizar la iluminación de una imagen. 
A veces una iluminación descompensada responde a fines artísticos, como en los cuadros de Caravaggio, pero otras veces no es intencionada y queda antiestético que, por ejemplo, según vamos de derecha a izquierda de una foto esté más oscurecida por la posición del Sol u otra fuente puntual. 
Para eliminar gradientes de iluminación, se parte del \emph{modelo de iluminación-reflectividad} y se usa lo que tiene el pomposo nombre de \emph{filtrado homomorfo}. Puedes encontrar la descripción al final de \S 2.2.4 en \cite{chamizo_sp}, desde ``The second example\dots'' en la página 83. Mira también en mi web la página:
\begin{center}
 \url{http://matematicas.uam.es/~fernando.chamizo/dark/d_homom.html}
\end{center}
donde hay imágenes de tamaño mayor y un programa del que usarás una parte. 



\ejer
Escriba unas pocas líneas resumiendo el modelo y el método de filtrado. Si te es posible, da un vistazo a alguna de las referencias. 
\medskip

Ahora vamos con la tarea que consiste en experimentar con el filtro
usando un programa sencillo hecho por mí.

Copia la función {\small\verb!hom_ff!} que está al final de la página web anterior, en un fichero llamado {\small\verb!hom_ff.m!}. Como se explica en la línea de comentarios, los argumentos segundo tercero y cuarto son $\sigma$, $\gamma_L$ y $\gamma_H$ y el último es el nombre del fichero donde quieres que se guarde el resultado. Para utilizar esta función, considera este programa:
\begin{lstlisting}[language=Matlab]
ima = imread('original.png');
ima = rgb2gray(ima);  % solo si original.png es en color
imwrite( ima ,'original_gris.png'); % solo si original.png es en color
ima = im2double(ima);
epsi = 1/512;
hom_ff( ima, 60, 0, 3, epsi, 'resultado.png')
\end{lstlisting}
Esto aplica a la imagen de partida \code{original.png} el filtro con $\sigma=60$, $\gamma_L=0$ y $\gamma_H=3$ para obtener \code{resultado.png}. Solo funciona con imágenes cuadradas y te recomiendo que las uses de tamaño $256\times 256$ o $512\times 512$ (en el primer caso quizá sea bueno duplicar \code{epsi}). El filtro está pensado para imágenes en blanco y negro (tonos de gris) pero lo mismo no sabes cómo convertir una imagen tomada con tu cámara a este formato con Photoshop o GIMP. Por eso he incluido las líneas 2 y 3 que pasan el original a blanco y negro y guardan una copia. Si ya lo estuviera, como en el ejemplo de mi web, esas dos líneas sobran. He puesto la extensión \code{png} pero puedes usar \code{jpg}, que es lo común para fotos, y otros formatos que te permita \code{matlab}. Prueba con una imagen cuadrada cualquiera que el anterior programa no te da errores, aunque el resultado sea absurdo. La imagen debe estar en el mismo directorio que el programa o en un {\sl path} que reconozca \code{matlab}. 

\ejer
Haz una foto (preferible) o busca una imagen en la web que tenga un gradiente suave de iluminación, es decir, que esté más oscurecida a un lado que a otro sin cambios demasiado bruscos. Recórtala para que sea cuadrada (recuerda los tamaños recomendados) y empieza a practicar con los parámetros $\sigma$, $\gamma_L$ y $\gamma_H$ en el programa anterior hasta conseguir que \code{resultado.png} tenga una iluminación más nivelada que el original, o su versión en blanco y negro. Reproduce los dos o tres resultados que te parezcan más atractivos. 
En una imagen artística con un gradiente exagerado, los parámetros pueden ser muy diferentes de los del ejemplo con texto de mi web y los resultados no muy satisfactorios porque el modelo es solo aproximado.
%\footnote{Nada más que para jugar, he tomado {\sl``La crucifixión de San Pedro''} de Caravaggio que es 
%de los cuadros del autor con claroscuro algo menos exagerado y, pagando con pérdida de calidad, tomando
%$\sigma=0.8$, $\gamma_L=0$ y $\gamma_H=1$ he logrado ver detalles que no se distinguían en mi pantalla 
%en la imagen original.}.

Seguramente tengas que probar con varias fotos y parámetros antes de encontrar resultados que te gusten. En eso consiste la tarea. Cuanto mejor entiendas el filtro, menos al azar serán las pruebas. 


\paragraph{2. El filtro de mediana.}
El \emph{ruido sal y pimienta} \cite{wiki_s_and_p} en una imagen (siempre en esta hoja, en blanco y negro) es una contaminación que consiste en que aleatoriamente algunos píxeles pasan a ser blancos (sal) y otros negros (pimienta). Esto puede ocurrir cuando hay una saturación esporádica de la señal al transmitir la información. El \emph{filtro de mediana} lo reduce de forma bastante espectacular sustituyendo el valor de cada píxel por la mediana de los valores en el cuadrado $(2r+1)\times (2r+1)$ que lo rodea con $r$ píxeles en cada dirección. Habitualmente se toma~$r=1$. 
\[
 r=1\quad\to\quad
 \begin{matrix}
  \square\square\square
  \\[-5pt]
  \square\blacksquare\square
  \\[-5pt]
  \square\square\square
 \end{matrix}
 \qquad\ \qquad\ \qquad
 r=2\quad\to\quad
 \begin{matrix}
  \square\square\square\square\square
  \\[-5pt]
  \square\square\square\square\square
  \\[-5pt]
  \square\square\blacksquare\square\square
  \\[-5pt]
  \square\square\square\square\square
  \\[-5pt]
  \square\square\square\square\square
 \end{matrix}
\]

El siguiente programa contamina una imagen con un ruido de sal y pimienta de densidad~\code{d} y 
aplica el filtro con $r$ igual a \code{r}. Tras el apartado anterior ya tendrás al menos una imagen en tonos de gris. Si prefieres partir de una en color, como antes,  quita el \code{\%} de comentario de las líneas~2 y~3.
\begin{lstlisting}[language=Matlab]
ima = imread('original.png');
%ima = rgb2gray(ima);
%imwrite( ima ,'original_gris.png');
ima = im2double(ima);


% Imagen con ruido
d = 0.15; %densidad de ruido
imagn = imnoise(ima,'salt & pepper',d);
imwrite(imagn, 'noisy.png')

% Filtro mediano
r = 1;
imagm = medfilt2(imagn,[2*r+1,2*r+1]);
imwrite(imagm, 'resultado.png')
\end{lstlisting}

\ejer
Practica con algunos ejemplos manteniendo \code{r=1}. 
Tienes que comparar la imagen contaminada \code{noisy.png} con la limpieza conseguida en \code{resultado.png}.
Examina el resultado cuando \code{d} varía del $15\%$ indicado de píxeles dañados a otras proporciones. Escribe unas líneas explicando por qué el filtro de mediana funciona bien con el ruido de sal y pimienta cuando su densidad no es muy grande. 
\medskip

En principio, si la densidad del ruido es mayor, sería natural tomar \code{r} mayor que $1$ pero eso da problemas. 

\ejer
Comprueba con ejemplos qué tipo de problemas aparecen y justifica a qué se deben. 
\medskip

Añadiendo las siguientes líneas al programa:
\begin{lstlisting}[language=Matlab]
% Filtro de mediana solo para píxeles negros o blancos
mask = (imagn.*(1-imagn)==0);
imagm2 = (1-mask).*imagn + mask.*medfilt2(imagn, [2*r+1,2*r+1] );
imwrite(imagm2, 'resultado2.png');
\end{lstlisting}
se aplica el filtro de mediana  solo a los píxeles que valen cero o uno que son los únicos que pueden estar afectados por el ruido. 

\ejer
Compruebe que el nuevo filtro con \code{r=2} funciona mejor que el anterior para \code{r=1} y \code{r=2} con densidades algo mayores. Busca una razón para esta mejora. 



\paragraph{3. Reducción de ruido periódico.}
Ahora nuestro objetivo es limpiar una imagen de ruido en forma de {\sl tonos puros} periódicos que se superponen a los valores dando lugar a ciertas oscilaciones. En la electrónica analógica esto no es tan extraño por las posibles interferencias de componentes funcionando a baja frecuencia. Quizá te suene haber visto algo de este tipo en imágenes de televisión antiguas. 
En la página de mi sitio web:
\begin{center}
 \url{http://matematicas.uam.es/~fernando.chamizo/dark/d_periodic.html}
\end{center}
traté el problema y escribí un programa. La idea es que el ruido se tiene que manifestar en picos sueltos en la DCT y un pico se puede eliminar con un filtro de mediana. Lo malo es que la DCT tiene de por sí picos, por ejemplo para la frecuencia nula, y hay que buscar algo que a ser posible no toque en absoluto estos picos buenos. Los correspondientes al ruido deben ser muy finos, porque la frecuencia está muy localizada, por ello impuse como regla para detectar un pico de ruido que el valor en un punto sea mucho mayor que la mediana de los de alrededor.

\ejer
Por favor, relee el párrafo anterior y también lo mencionado en la página web hasta que entiendas el procedimiento. Si necesitas más detalles, mira el código de la función \code{denoip}, incluido al final de la página. 
\medskip

Considera el siguiente programa que aplica la idea para cierto \code{r} que da el tamaño del filtro de mediana y para cierto \code{alpha} que dice cuántas veces mayor debe ser un valor que la mediana para que el pico se considere malo.  La imagen debe ser de tamaño $512\times 512$. 
\begin{lstlisting}[language=Matlab]
ima = imread('original.png');
%ima = rgb2gray(ima);
%imwrite( ima ,'original_gris.png');
ima = im2double(ima);

% Imagen con ruido
t = linspace(0,1 ,512);
[X, Y] = meshgrid(t, t);
iman = ima + 0.2*sin(16*pi*(X-Y));
imwrite( im2double(mat2gray(iman)) ,'noisy.png' );

% Limpieza
r = ...
alpha = ...
res = denoip( iman, r, alpha );
imwrite( res ,'resultado.png' );
\end{lstlisting}



\ejer
Escoge una imagen y experimenta con valores para los puntos suspensivos en las líneas~13 y~14 hasta que te parezca que no puedes reducir más el ruido. No estoy seguro pero sospecho que los imágenes que no tienen grandes zonas totalmente blancas o negras deberían ser más adecuadas. Si no consigues una reducción apreciable, prueba con diferentes imágenes. 
También, si quieres, puedes cambiar la contaminación de la línea 9 por otra. 






\

\begin{center}
\noindent\rule{0.95\textwidth}{1.1pt}
\end{center}


\paragraph{Tarea a entregar.}
El documento a entregar debe comenzar la parte teórica del análisis de Fourier discreto con un ejemplo, descrita en el tercer ejercicio. Después, debes incluir los experimentos del filtrado homomorfo, el filtro de mediana y la eliminación de interferencias periódicas, todos ellos precedidos con algunas explicaciones con el grado de profundidad que prefieras. 

No te pongo un límite para el tamaño del documento porque, más que en otras ocasiones, hay muchas posibilidades para modificarlo según las necesidades de tu TFG, por ejemplo escalando las imágenes o reduciendo el número de ejemplos. Para las partes de código que están en mi web, puedes dar la referencia o copiarlo, lo que prefieras. 

\begin{center}
\noindent\rule{0.95\textwidth}{1.1pt}
\end{center}

\begin{thebibliography}{{Wik}20}

\bibitem[Cha20]{chamizo_sp}
F.~Chamizo.
\newblock A course on signal processing.
\newblock
  {\small\url{http://matematicas.uam.es/~fernando.chamizo/libreria/libreria.html}},
  2020.

\bibitem[Sal02]{salomon}
D.~Salomon.
\newblock {\em A guide to data compression methods}.
\newblock Springer-Verlag, 2002.

\bibitem[Wal08]{walker}
J.~S. Walker.
\newblock {\em A primer on wavelets and their scientific applications}.
\newblock Studies in Advanced Mathematics. Chapman \& Hall/CRC, Boca Raton, FL,
  second edition, 2008.

\bibitem[{Wik}20]{wiki_s_and_p}
{Wikipedia contributors}.
\newblock Salt-and-pepper noise --- {Wikipedia}{,} the free encyclopedia.
\newblock
  \url{https://en.wikipedia.org/w/index.php?title=Salt-and-pepper_noise&oldid=943959344},
  2020.
\newblock [Online; accessed 21-March-2021].

\end{thebibliography}


%\bibliography{bibtfg20}
%\bibliographystyle{alpha}

\end{document}


