\documentclass[11pt,leqno]{article}
\usepackage{amssymb, amsmath, amsthm} 
\usepackage{latexsym} 
\usepackage{graphics} 
\usepackage[hyphens]{url}
%\def\UrlBreaks{\do\/\do-}
\expandafter\def\expandafter\UrlBreaks\expandafter{\UrlBreaks% save the current one
  \do\a\do\b\do\c\do\d\do\e\do\f\do\g\do\h\do\i\do\j%
  \do\k\do\l\do\m\do\n\do\o\do\p\do\q\do\r\do\s\do\t%
  \do\u\do\v\do\w\do\x\do\y\do\z\do\A\do\B\do\C\do\D%
  \do\E\do\F\do\G\do\H\do\I\do\J\do\K\do\L\do\M\do\N%
  \do\O\do\P\do\Q\do\R\do\S\do\T\do\U\do\V\do\W\do\X%
  \do\Y\do\Z\do\*\do\-\do\~\do\'\do\"\do\-}%

\usepackage{xcolor}
\usepackage{hyperref} 

\usepackage{tikz}
\usetikzlibrary{babel,angles,quotes}
\usetikzlibrary{shadings}
% \pgfdeclareradialshading{mgradi}{\pgfpoint{0.5cm}{0.5cm}}%
%   {rgb(0cm)=(0.9,0,0);
%    rgb(0.7cm)=(0.7,0,0);
%    rgb(1cm)=(0.5,0,0);
%    rgb(1.05cm)=(1,1,1)}
\pgfdeclareradialshading{mgradi}{\pgfpoint{0bp}{0bp}}%
  {rgb(0bp)=(0.0,0,0.0);
   rgb(0.2bp)=(0.04,0.04,0.04);
   rgb(0.4bp)=(0.16,0.16,0.16);
   rgb(0.6bp)=(0.36,0.36,0.26);
   rgb(0.8bp)=(0.64,0.64,0.64);
   rgb(1bp)=(1,1,1)}
%\pgfuseshading{mgradi}

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


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

\decimalpoint

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

\theoremstyle{plain}

\newtheorem{theorem}{Teorema}
\newtheorem{corollary}[theorem]{Corolario}
\newtheorem{lemma}[theorem]{Lema}

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 	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{\ejer}{\ \par\vspace{-7pt}{\bf\number\value{ejeno})\ }\addtocounter{ejeno}{1}}


\usepackage{listings}


\lstdefinelanguage{Sage}[]{Python}
{morekeywords={False,sage,True},sensitive=true}

\lstset{
  frame=none,
  showtabs=False,
%  showtabs=True,
  showspaces=False,
  showstringspaces=False,
%  showstringspaces=True,
  commentstyle={\ttfamily\color{green!30!black}},
  keywordstyle={\ttfamily\color{blue!30!black}\bfseries},
  stringstyle={\ttfamily\color{gray!30!black}\bfseries},
%  language=Sage,
%  basicstyle={\fontsize{10pt}{10pt}\ttfamily},
  basicstyle=\footnotesize,
%  basicstyle={\fontsize{10pt}{10pt}\rmfamily},
  aboveskip=0.6em,
  belowskip=0.5em,
%  numbers=left,
%  numberstyle=\footnotesize
  xleftmargin=25pt
}

%\lstset{framexleftmargin=5mm, frame=shadowbox, rulesepcolor=\color{blue}}



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


En esta hoja última hoja vamos a ver dos aplicaciones relacionadas con cosas que has estudiado hasta ahora. Por los posibles apuros de tiempo que puedas tener, he intentando que haya tareas un poco abiertas, dependiendo de cuánto quieras experimentar con el ordenador. 

\paragraph{1. El formato de imágenes JPEG.}
La mayoría de las fotos se almacenan en nuestros dispositivos electrónicos en un formato que corresponde a las extensiones \texttt{.jpg} o \texttt{.jpeg}. En él participa de forma crucial la DCT que estudiaste en la primera hoja. También aparece la codificación óptima a la que me referí en la última. Comencemos por este último punto que tiene interés independiente. 


\ejer
Lee cómo se lleva a cabo la codificación de Huffman en la sección 4.2.2 de \cite{ch_sp} y escríbelo con tus propias palabras. No hace falta que incluyas la prueba de que es óptima, solo el método. Hay también una explicación elemental en el texto clásico \cite{renyi}, por si lo localizas y quieres darle un vistazo. Para ver que lo has entendido, establece la codificación de Huffman para un conjunto $S=\{s_1,\dots,s_8\}$ con probabilidades dadas por esta tabla:
\[
 \begin{array}{r|c|c|c|c|c|c|c|c|}
  &s_1& s_2& s_3& s_4& s_5& s_6& s_7& s_8
  \\
  \hline
  \text{Prob.}&
  2/5& 7/40& 1/10& 1/10& 1/10& 1/20& 1/20& 1/40
  \\
  \hline
 \end{array}
\]

\ejer
Lee acerca del formato JPEG en la sección 2.2.2 de \cite{ch_sp}, limitándote solo al caso de imágenes en blanco y negro y explícalo con tus propias palabras. Por supuesto, puedes manejar otras fuentes de información. En el apartado {\sl ``JPEG codec example''} de \cite{wiki_jpeg} hay ejemplos numéricos y en el resto del artículo algunas cosas técnicas. 
Si te sobra espacio o tienes curiosidad, en 
\cite{austin} se comenta brevemente el uso del color y el formato JPEG2000, que es mejor, pero ha caído en el olvido. 

\medskip

Ahora te propongo hacer dos experimentos. 
En el primero, si no sabes usar {\sl software} (Photoshop, GIMP,\dots) para crear las imágenes, limítate a dar una explicación del fenómeno con lo que has aprendido y, si no hay más remedio, te paso yo una tabla. Nota que en informática los {\sl pixels} se empiezan a contar desde~$0$, no desde~$1$, se suele denominar $(0,0)$ al de la esquina superior izquierda. 

\ejer
Crea una imagen de tamaño $256\times 256$ consistente en un cuadrado  $128\times 128$ negro sobre fondo blanco tal que uno de los vértices esté justamente en la posición $(x_0,y_0)= (32,32)$. Cambia $(x_0,y_0)$ por otros valores y haz una pequeña tabla mostrando el tamaño en bytes del fichero correspondiente cuando se guarda en formato JPEG. Explica por qué si $x_0$ o $y_0$ es múltiplo de~$8$ el tamaño se reduce, y más todavía si lo son ambos.  
\medskip

El siguiente experimento es muy abierto y consiste en que juegues como quieras con un programa que he hecho. El código lo puedes copiar de la fuente \LaTeX{}.


\ejer
Dada una imagen y una proporción $0<p<1$,
el siguiente programa {\sf matlab} (con una función)
estudia el efecto de reducir a cero un número de coeficientes DCT correspondiente a esa proporción o hacerlo directamente con los {\sl pixels} (cero corresponde a negro). Practica su efecto con diferentes elecciones de $p$ y tipos de imágenes (te sugiero que no sean de gran tamaño), por ejemplo con gradaciones (fotos) o con bordes (dibujos geométricos). Comenta los resultados que te parezcan significativos incluyendo las imágenes en tu trabajo. Nota que descartar casi todos los coeficientes DCT, por ejemplo con $p=0.9$, todavía da una imagen reconocible, pero casi cualquier $p$ tiene un efecto sobre los bordes abruptos. Trata de dar alguna explicación matemática de lo que observes.


\lstset{language=matlab}
\begin{lstlisting}
imag = imread('nombre_imagen.extension');
imag = rgb2gray(imag);
imag = im2double(imag);

figure(1) % Imagen original, pero en blanco y negro
imshow( imag )

p = 0.8;

figure(2) % 100p por ciento de los coef. DCT acero (los menores)
img = quant_image(imag, p);
imshow( img )

figure(3) % 100p por ciento de los pixels al azar -> negro
img = imag.* (rand(size(imag))>p);
imshow( img )

figure(4) %  100p por ciento de los pixels (los más oscuros) -> negro
img = imag.* (imag> quantile( reshape(imag,1,[]), p));
imshow( img )
\end{lstlisting}
\noindent
Este programa llama a
\verb!quant_image! que debes poner en
un fichero aparte \verb!quant_image.m! (en el mismo directorio).
\lstset{language=matlab}
\begin{lstlisting}
function [img] = quant_image(imag, p)
  T = dctmtx(8);
  % Aplica la DCT
  dct_fun = @(block_struct) T * block_struct.data * T';
  img = blockproc(imag, [8 8], dct_fun, 'PadPartialBlocks',true);
  
  % Calcula el umbral
  thrs = quantile( reshape(img,1,[]), p);
  % Pone a cero los que no superen el umbral
  mask_fun = @(block_struct) (block_struct.data>thrs) .* block_struct.data;
  img = blockproc(img, [8 8], mask_fun, 'PadPartialBlocks',true);
  
  % Aplica la DCT inversa
  idct_fun = @(block_struct) T' * block_struct.data * T;
  img = blockproc(img, [8 8], idct_fun, 'PadPartialBlocks',true);
end 
\end{lstlisting}


\paragraph{2. Tomografías radiales.}
Casi cualquier objeto es semitransparente para los rayos $X$. Como es natural, estos rayos se atenúan en función de la densidad de los puntos que atraviesan. Al estar interesados en secciones de objetos sólidos, consideramos que la densidad es una función $\rho:\R^2\longrightarrow\R$.  
Denotando $\ell_x$ a la recta vertical que pasa por $(x,0)$, al tomar una radiografía unidimensional desde arriba podremos deducir la función 
\[
 S(x)=\int_{\ell_x}\rho=\int_{-\infty}^\infty \rho(x,y)\, dy.
\]
Diremos que $S$ es la \emph{sombra} de la densidad (me he inventado el nombre, si quieres cámbialo). 
En las TAC (\emph{tomografías axiales computarizadas}) se plantea el problema matemático de obtener la densidad a partir de las sombras (radiografías unidimensionales) en diferentes direcciones. Matemáticamente se necesitarían radiografías en todas la direcciones y en la práctica se usan muchas (\cite{smith} es un viejo artículo sobre ello). El caso en que $\rho$ es radial, esto es, solo depende de $x^2+y^2$, es especial porque las sombras en diferentes direcciones coinciden (¿lo ves claro?) y basta $S(x)$ para obtener $\rho$. Es decir, se deduce la estructura interna con una sola radiografía. Este caso simplificado es el que vamos a abordar. Se podría resolver por medio de la llamada \emph{transformada de Abel}~\cite{wiki_tabel} (igual que el caso general se puede resolver con la de Radon~\cite{wiki_tradon}). Aquí seguiremos una estrategia distinta para que practiques con algunas de las cosas que has estudiado y hagas algunos cálculos numéricos. 
\medskip

Por ejemplo, las densidades $\rho_1(x,y)$, $\rho_2(x,y)$ y $\rho_3(x,y)$, dadas por 
\[
 \begin{cases}
  1&\text{si }x^2+y^2<1,
  \\
  0&\text{si }x^2+y^2>1,
 \end{cases}
 \qquad
 \begin{cases}
  0&\text{si }x^2+y^2<\frac 14,
  \\
  1&\text{si }\frac 14<x^2+y^2<1,
  \\
  0&\text{si }x^2+y^2>1,
 \end{cases}
 \qquad
 \begin{cases}
  1-x^2-y^2&\text{si }x^2+y^2<1,
  \\
  0&\text{si }x^2+y^2>1,
 \end{cases}
\]
corresponden a un disco unidad homogéneo, a otro en el que se le ha vaciado su parte central y a otro que es más denso cuando más cerca del centro estemos. 

\hspace{-1cm}\begin{tikzpicture}[x=1pt,y=1pt,scale=30]
          \clip (-2,2) rectangle (13.3,6);
          \draw[black,fill=black] (1,4.7) circle (1);
          \draw[black,fill=black] (6,4.7) circle (1);
          \draw[white,fill=white] (6,4.7) circle (0.5);
          \shade[shading=mgradi] (11,4.7) circle (25);
          \node at (1,3) {\scriptsize$\displaystyle S_1(x)=\begin{cases}
                                                   2\sqrt{1-x^2}&\text{ si $|x|<1$},
                                                   \\
                                                   0&\text{ si $|x|\ge 1$},
                                                  \end{cases}$};
          \node at (6,3) {\scriptsize$\displaystyle S_2(x)=S_1(x)-\frac 12 S_1(2x)$};
          \node at (11,3) {\scriptsize$\displaystyle S_3(x)=\begin{cases}
                                                   \frac 43(1-x^2)^{3/2}&\text{ si $|x|<1$},
                                                   \\
                                                   0&\text{ si $|x|\ge 1$}.
                                                  \end{cases}$};
\end{tikzpicture}

Por cierto, si quieres copiar el código que genera estas imágenes de la fuente \LaTeX{}, quizá tengas que cargar el paquete TikZ y alguna de sus librerías. Mira la cabecera de este documento.

\ejer
Comprueba que las funciones sombra en estas tres situaciones  son las indicadas en la figura anterior. 
\medskip

Es conocido que si $\rho$ es radial, su transformada de Fourier (en dos variables) también lo es~\cite{folland}. Este hecho se utiliza en el siguiente problema.

\ejer
Explica por qué para $\rho$ radial, $\rho(x,y)=\rho(t,0)$ con $t=\sqrt{x^2+y^2}$ y $S$ es siempre una función par. Justifica en esta situación las siguientes igualdades sin preocuparte de la regularidad:
\begin{eqnarray*}
 \rho(t,0)
 &=&
 \int_0^{2\pi} \int_0^\infty 
 r\,\widehat{\rho}(r\cos\theta,r\sen\theta)\, e(tr\cos\theta)\, dr d\theta
 \\
 &=&
 \int_0^{\pi} \int_{-\infty}^\infty 
 |r|\,\widehat{\rho}(r,0)\, e(tr\cos\theta)\, dr d\theta
 =
 \int_0^{\pi} \int_{-\infty}^\infty 
 |r|\,\widehat{S}(r)\, e(tr\cos\theta)\, dr d\theta.
\end{eqnarray*}
{\sf Indicación:} Utiliza la fórmula de inversión en ambas variables y pasa a coordenadas polares. ¿Por qué $\widehat{\rho}(r,0)$ es real y par?
\medskip

En principio la fórmula anterior resuelve matemáticamente el problema: a partir de $S$ podemos recuperar $\rho$. Sin embargo, sería duro aproximar numéricamente la transformada de Fourier de $S$ y después hacer una integral doble con ella. Vamos a aprovechas algunas cosas que conoces para diseñar un método numérico más eficiente. 

\ejer
Recordando la serie de Fourier de $G(x)$ en la primera hoja, explica por qué 
se cumple 
\[
 |r|
 =
 N\sum_{k=0}^\infty a_k\cos\Big(\frac{2\pi k r}N\Big)
 \quad\text{para}\quad  r\in \Big[-\frac N2,\frac N2\Big]
 \quad\text{con}\quad
 a_k
 =
 \begin{cases}
  \frac14 &\text{si }k=0,
  \\ 
  -\frac{2}{\pi^2 k^2} &\text{si }k\text{ es impar},
  \\ 
  0 &\text{en el resto}.
 \end{cases}
\]

\ejer
Con ello y usando las propiedades de la transformada de Fourier, justifica por qué para $N$ grande cabe esperar la aproximación  
\[
 \int_{-\infty}^\infty 
 |r|\,\widehat{S}(r)\, e(tr\cos\theta)\, dr
 \approx
 \frac N2
 \sum_{k=0}^\infty a_k\, \big(S(t\cos\theta+k/N)+S(-t\cos\theta+k/N)\big).
\]

\ejer 
Truncando la serie hasta $2N$ y aplicando  a la integral en $\theta$ la simetría $\cos(\pi-\theta)=-\cos\theta$ junto con 
la regla del trapecio \cite{wiki_trapez}, obtén 
\[
 \rho(x,y)\approx
 \sum_{0\le \ell<\pi N}
 \sum_{k=0}^{2N} a_k\, S\Big(\frac kN + t\cos\frac{\ell}N\Big)
 \qquad\text{con}\quad t=\sqrt{x^2+y^2}.
\]

\ejer
Escribe un programa que efectúe el doble sumatorio anterior. Tomando $N\ge 20$ tan grande como permita tu ordenador y tu paciencia, representa gráficamente las aproximaciones de $\rho(x,0)$ con $x\in [0,3/2]$ correspondientes a $S=S_1,\, S_2,\, S_3$ superponiendo la gráfica exacta de $\rho(x,0)=\rho_j(x,0)$ para compararlas. 
{\sf Indicación:} Para no alargar mucho el tiempo de ejecución, al dibujar la gráfica de la aproximación $\rho(x,0)$ no evalúes en demasiados valores de $x$. Con $50$ o así, ya verás el aspecto que tiene. 
\medskip

Como he dicho antes, el caso de las tomografías radiales se puede resolver de otra forma con la transformada de Abel \cite{wiki_tabel}. No olvides mencionarlo  en tu trabajo no vaya a ser que alguien del tribunal lo sepa y piense que estás haciendo una tontería. El propósito aquí es ilustrar lo que viste en el primer capítulo y transformarlo en un método numérico. 

\


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


\paragraph{Tarea a entregar.}
Escribe un documento que combine las soluciones de los ejercicios anteriores. 
El resultado será el cuarto y último capítulo de tu TFG  bajo el título {\sl Dos aplicaciones} o cualquier variante que te parezca adecuada. La extensión la dejo a tu elección ya que puede variar mucho dependiendo de las imágenes que incluyas. En caso de que sea grande, parte de ellas podrían ir a un apéndice. 

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

\begin{thebibliography}{1}

\bibitem{austin}
D.~Austin.
\newblock What is{$\dots$} {JPEG}?
\newblock {\em Notices Amer. Math. Soc.}, 55(2):226--229, 2008.

\bibitem{ch_sp}
F.~Chamizo.
\newblock A course on signal processing.
\newblock
  {\small\url{http://matematicas.uam.es/~fernando.chamizo/libreria/fich/signal_processing.pdf}},
  Curso 2017/2018.

\bibitem{folland}
G.~B. Folland.
\newblock {\em Fourier analysis and its applications}.
\newblock Pacific Grove, CA: Wadsworth \& Brooks/Cole Advanced Books \&
  Software, 1992.

\bibitem{renyi}
A.~R{\'e}nyi.
\newblock {\em A diary on information theory}.
\newblock Akad\'emiai Kiad\'o (Publishing House of the Hungarian Academy of
  Sciences), Budapest, 1984.

\bibitem{smith}
K.~T. Smith, D.~C. Solmon, and S.~L. Wagner.
\newblock Practical and mathematical aspects of the problem of reconstructing
  objects from radiographs.
\newblock {\em Bull. Am. Math. Soc.}, 83:1227--1270, 1977.

\bibitem{wiki_tabel}
{Wikipedia contributors}.
\newblock Abel transform --- {Wikipedia}{,} the free encyclopedia.
\newblock
  {\small\url{https://en.wikipedia.org/w/index.php?title=Abel_transform&oldid=1239227961}},
  2024.
\newblock [Online; accessed 7-April-2025].

\bibitem{wiki_tradon}
{Wikipedia contributors}.
\newblock Radon transform --- {Wikipedia}{,} the free encyclopedia.
\newblock
  {\small\url{https://en.wikipedia.org/w/index.php?title=Radon_transform&oldid=1250073731}},
  2024.
\newblock [Online; accessed 7-April-2025].

\bibitem{wiki_jpeg}
{Wikipedia contributors}.
\newblock {JPEG} --- {Wikipedia}{,} the free encyclopedia.
\newblock
  \url{https://en.wikipedia.org/w/index.php?title=JPEG&oldid=1284340128}, 2025.
\newblock [Online; accessed 9-April-2025].

\bibitem{wiki_trapez}
{Wikipedia contributors}.
\newblock Trapezoidal rule --- {Wikipedia}{,} the free encyclopedia.
\newblock
  {\small\url{https://en.wikipedia.org/w/index.php?title=Trapezoidal_rule&oldid=1279480977}},
  2025.
\newblock [Online; accessed 8-April-2025].

\end{thebibliography}


%\bibliography{bibtfg24MC}
%\bibliographystyle{plain}

\end{document}

