5  Transformadas y Compresión

En los capítulos anteriores, todas las operaciones se realizaron en el dominio espacial, en el que los algoritmos actúan directamente sobre los valores de intensidad de los píxeles.

En este capítulo se presentará un enfoque complementario: el dominio de la frecuencia, en el cual la imagen se representa mediante las variaciones espaciales de intensidad, y no solo por los valores individuales de los píxeles.

El concepto de frecuencia espacial describe la rapidez con la que la intensidad varía a lo largo de la imagen. Las variaciones lentas corresponden a bajas frecuencias, mientras que los bordes, los detalles finos y los ruidos corresponden a altas frecuencias.

Esta representación se basa en el hecho de que cualquier imagen digital discreta puede descomponerse en una combinación de funciones ortogonales. La Transformada de Fourier utiliza una base de exponenciales complejas bidimensionales (equivalentes a senoides con orientación y frecuencia específicas). Otras transformadas, como la Transformada de Cosenos (DCT) y la Transformada Wavelet (DWT), utilizan diferentes familias de funciones de base — cosenos bidimensionales en el caso de la DCT, y funciones con soporte compacto en el caso de las wavelets.

Entre las principales aplicaciones de esta representación se destacan:

1. Filtrado en el dominio de la frecuencia, para atenuar o realzar determinadas bandas de frecuencia;
2. Análisis multirresolución mediante transformadas wavelet, que representa estructuras en diferentes escalas;
3. Compresión de imágenes, mediante la reducción del número de coeficientes necesarios para representar la imagen.

5.1 Objetivos

Al concluir este capítulo, usted será capaz de:

• Interpretar el espectro de Fourier de una imagen, distinguiendo magnitud, fase y componentes de frecuencia;
• Aplicar el Teorema de la Convolución para realizar filtrado en el dominio de la frecuencia utilizando la Transformada Rápida de Fourier (FFT);
• Diseñar y analizar filtros en el dominio de la frecuencia, comprendiendo el funcionamiento de filtros pasa-baja, pasa-alta y notch;
• Comprender el análisis multirresolución mediante transformadas wavelet y su aplicación en la representación jerárquica de imágenes;
• Describir el proceso de compresión de imágenes, incluyendo la Transformada Discreta del Coseno (DCT) y la cuantización de los coeficientes;
• Seleccionar formatos de almacenamiento de imágenes, como JPEG, PNG y WebP, de acuerdo con los requisitos de la aplicación.

5.2 Configuración del Entorno

import os, urllib.request

os.makedirs("tmp/state", exist_ok=True)  # artefactos de build de la ruta C++

url = "https://raw.githubusercontent.com/fzampirolli/pdi-vc/master/morph/config.py"
if not os.path.exists("config.py"):
    urllib.request.urlretrieve(url, "config.py")

# Kernel Python incluso en la ruta C++. cpp=True descarga morph.hpp + stb; las
# celdas %%writefile *.cpp de este capítulo compilan CON OpenCV
# (-DMM_USE_OPENCV + pkg-config opencv4).
import config
config.setup(cpp=True)
from morph import mm
import numpy as np

✅ Entorno listo. Morph: 1.1.9 | OpenCV: 5.0.0

5.3 Transformada de Fourier Discreta 2D

El análisis de Fourier se basa en el principio de que cualquier señal periódica puede representarse como una suma de funciones senoidales con diferentes frecuencias, amplitudes y fases. Este concepto también se aplica a las imágenes digitales, permitiendo representarlas en el dominio de la frecuencia en lugar del dominio espacial.

La Figura 5.1 ilustra esta descomposición para una señal unidimensional. En el caso de una imagen, la Transformada Discreta de Fourier (DFT) convierte la matriz de intensidades \(f(x,y)\) en un conjunto de coeficientes que describe la contribución de las diferentes frecuencias espaciales presentes en la imagen.

%%writefile tmp/fig_decomposicao_1d.cpp
#define MM_OUT "tmp/fig_decomposicao_1d.png"
//| label: fig-decomposicao-1d
//| fig-cap: "Decomposição de Fourier 1D: uma onda quadrada (linha tracejada) é aproximada pela soma das primeiras senoides (linhas coloridas). Quanto mais termos, melhor a aproximação."
//| echo: false
//| output: true

#include <opencv2/opencv.hpp>
#include "morph.hpp"
#include <vector>
#include <string>
#include <cmath>
#include <filesystem>

int main() {
    // Crear vector x: 400 puntos de 0 a 2*pi
    std::vector<double> x(400);
    for (int i = 0; i < 400; ++i) {
        x[i] = 2.0 * M_PI * i / 399.0;
    }

    // Onda cuadrada: 1 para x < pi, -1 en caso contrario
    std::vector<double> square(400);
    for (int i = 0; i < 400; ++i) {
        square[i] = (x[i] < M_PI) ? 1.0 : -1.0;
    }

    // Suma acumulada y armónicos
    std::vector<double> soma(400, 0.0);
    std::vector<std::vector<double>> harms;
    for (int n = 1; n <= 3; ++n) {
        std::vector<double> h(400);
        for (int i = 0; i < 400; ++i) {
            h[i] = (4.0 / M_PI) * (1.0 / (2 * n - 1)) * std::sin((2 * n - 1) * x[i]);
            soma[i] += h[i];
        }
        harms.push_back(h);
    }

    // Curvas y etiquetas
    std::vector<std::vector<double>> ys = {square, harms[0], harms[1], harms[2], soma};
    std::vector<std::string> labels = {
        "Onda quadrada ideal", "1a harmonica", "3a harmonica", "5a harmonica", "Soma (3 primeiras)"
    };
    std::vector<cv::Scalar> colors = {
        cv::Scalar(40, 40, 40), cv::Scalar(60, 160, 80), cv::Scalar(180, 120, 60),
        cv::Scalar(150, 80, 160), cv::Scalar(60, 60, 220)
    };

    // Generar gráfico
    mm::Image chart = mm::lineChart(
        x, ys, labels, colors,
        "Sintese de Fourier: de senos a uma onda quadrada",
        "Posicao", "Intensidade"
    );

    // Mostrar resultado
    std::vector<mm::Image> images = {chart};
    std::vector<std::string> titles = {"Decomposicao de Fourier 1D"};
    mm::show(images, MM_OUT, titles, 1);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(chart, "tmp/fig_decomposicao_1d_0.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_decomposicao_1d.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_decomposicao_1d.cpp -o tmp/fig_decomposicao_1d -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_decomposicao_1d \
  && test -f "tmp/fig_decomposicao_1d.png" \
  || echo "⚠ mm::show não gravou tmp/fig_decomposicao_1d.png"

[1] Decomposicao de Fourier 1D

try:
    mm.show(
        [
            mm.read("tmp/fig_decomposicao_1d_0.png"),
        ],
        titles=[
            'Decomposicao de Fourier 1D',
        ],
        cols=1,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_decomposicao_1d_0.png (ver a versao Python)")

Figura 5.1: Decomposição de Fourier 1D: uma onda quadrada (linha tracejada) é aproximada pela soma das primeiras senoides (linhas coloridas). Quanto mais termos, melhor a aproximação.

5.3.1 Simulador: Reconstruyendo Señales con Senoides

Antes de estudiar imágenes bidimensionales, el simulador de la Figura 5.2 ilustra el principio del análisis de Fourier para señales unidimensionales: una forma de onda puede aproximarse mediante la suma de senoides con diferentes frecuencias y amplitudes.

A medida que se añaden nuevos términos, la suma de las senoides (curva negra) se aproxima a la forma de onda de referencia (trazada). El gráfico inferior presenta el espectro de amplitudes, indicando la contribución de cada frecuencia a la reconstrucción de la señal.

TipActividad

Explore el simulador y responda:

1. ¿Cuántos términos son necesarios para obtener una buena aproximación de la onda cuadrada?
2. ¿Cuál de las tres formas de onda converge más rápidamente? Justifique su respuesta.
3. ¿Cómo se altera el espectro de amplitudes al cambiar la onda cuadrada por la triangular?

NotaRespuestas

1. ¿Cuántos términos son necesarios para una buena aproximación de la onda cuadrada?

Con aproximadamente 15 a 20 términos, la forma de la onda ya se aproxima bien a la referencia. Sin embargo, cerca de las discontinuidades permanece una pequeña oscilación, conocida como fenómeno de Gibbs, que no desaparece incluso con la adición de más términos.

2. ¿Cuál forma converge más rápidamente? ¿Por qué?

La onda triangular converge más rápidamente, pues las amplitudes de sus armónicos decaen más rápido que las de la onda cuadrada y la onda de diente de sierra. Como consecuencia, pocos términos ya producen una buena aproximación.

3. ¿Cómo cambia el espectro entre la onda cuadrada y la triangular?

Ambas poseen únicamente armónicos impares, pero, en la onda triangular, las amplitudes disminuyen mucho más rápidamente. Así, pocos armónicos son suficientes para reconstruir la señal con buena precisión.

∿ Simulador: Descomposición de Fourier 1D suma de senoides

Términos

1

Error RMS

–

Forma Objetivo

cuadrada

Forma Objetivo

N.º de Términos

1

Visualización

Figura 5.2: Simulador interactivo de la descomposición de Fourier 1D: visualización de la suma de senoides con diferentes frecuencias, amplitudes y fases. Añada términos y observe la convergencia hacia formas de onda arbitrarias.

5.3.2 Interpretación del espectro de frecuencia

Al aplicar la Transformada Discreta de Fourier (DFT) a una imagen y visualizar el módulo de sus coeficientes (ver Figura 5.5), se obtiene el espectro de magnitud, que muestra la distribución de las frecuencias espaciales presentes en la imagen.

El coeficiente ubicado en el origen de la DFT, denominado componente DC (Direct Current), corresponde a la frecuencia nula y representa la intensidad media de la imagen. Por convención, dicho coeficiente se almacena en la esquina superior izquierda del espectro. Para facilitar su interpretación, se aplica la operación FFT Shift, que desplaza la componente DC al centro de la imagen. Tras este desplazamiento, las bajas frecuencias se concentran en la región central, mientras que las altas frecuencias se sitúan cerca de los bordes, tal como resume la Tabla 5.1.

Tabla 5.1: Correspondencia entre las regiones del espectro de magnitud tras la aplicación del FFT Shift.

Región del espectro	Componentes predominantes	Ejemplos en la imagen
Centro (bajas frecuencias)	Variaciones espaciales lentas	Iluminación, regiones homogéneas y formas globales
Región intermedia (frecuencias medias)	Variaciones de escala intermedia	Texturas y patrones repetitivos
Bordes (altas frecuencias)	Variaciones espaciales rápidas	Contornos, detalles finos y ruido

Esta organización facilita la interpretación del espectro y el diseño de filtros. La atenuación de las bajas frecuencias reduce las variaciones globales de intensidad, mientras que la atenuación de las altas frecuencias suaviza la imagen al disminuir los detalles finos y parte del ruido.

5.3.3 El Experimento de la Rejilla: Construyendo una Imagen a partir de un Único Coeficiente

Antes de presentar la formulación matemática de la Transformada Discreta de Fourier (DFT), es útil analizar su inversa, denominada Transformada Discreta Inversa de Fourier (IDFT). Considere un espectro en el que todos los coeficientes sean nulos, excepto uno. Un ejemplo de esta construcción se presenta en el código de la Figura 5.3 y puede explorarse interactivamente en el simulador de la Figura 5.4..

La imagen reconstruida es una senoide bidimensional. La posición del coeficiente en el espectro determina su orientación y su frecuencia espacial, mientras que su magnitud y su fase definen, respectivamente, su amplitud y su desplazamiento espacial. Así, cada coeficiente de la DFT representa una componente senoidal, y la imagen original puede reconstruirse mediante la suma de todas estas componentes.

%%writefile tmp/fig_05_grade_2d.cpp
#define MM_OUT "tmp/fig_05_grade_2d.png"
//| label: fig-05-grade-2d
//| fig-cap: "Toda frequência no espectro (ponto isolado) corresponde a uma onda senoidal 2D rotacionada no domínio espacial."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include "morph.hpp"
#include <vector>
#include <string>
#include <filesystem>

int main() {
    int N_grid = 100;
    cv::Mat espectro_vazio = cv::Mat::zeros(N_grid, N_grid, CV_64FC2);

    // Acendendo um único ponto (frequência) fora do centro
    int u0 = 10, v0 = 5;
    espectro_vazio.at<cv::Vec2d>(N_grid/2 - v0, N_grid/2 - u0) = cv::Vec2d(1000, 0);

    // Retornando para o domínio espacial (IDFT)
    cv::Mat espectro_shifted;
    cv::Mat planes[2];
    cv::split(espectro_vazio, planes);
    // Manual fftshift
    int cx = N_grid / 2;
    int cy = N_grid / 2;
    cv::Mat q0(espectro_vazio, cv::Rect(0, 0, cx, cy));
    cv::Mat q1(espectro_vazio, cv::Rect(cx, 0, cx, cy));
    cv::Mat q2(espectro_vazio, cv::Rect(0, cy, cx, cy));
    cv::Mat q3(espectro_vazio, cv::Rect(cx, cy, cx, cy));
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    cv::Mat espectro_ifft;
    cv::idft(espectro_vazio, espectro_ifft, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT);

    // Manual ifftshift (swap again)
    cv::Mat q0b(espectro_ifft, cv::Rect(0, 0, cx, cy));
    cv::Mat q1b(espectro_ifft, cv::Rect(cx, 0, cx, cy));
    cv::Mat q2b(espectro_ifft, cv::Rect(0, cy, cx, cy));
    cv::Mat q3b(espectro_ifft, cv::Rect(cx, cy, cx, cy));
    q0b.copyTo(tmp);
    q3b.copyTo(q0b);
    tmp.copyTo(q3b);
    q1b.copyTo(tmp);
    q2b.copyTo(q1b);
    tmp.copyTo(q2b);

    cv::Mat onda_2d = espectro_ifft;

    cv::Mat onda_vis, espectro_vis;
    cv::normalize(onda_2d, onda_vis, 0, 255, cv::NORM_MINMAX, CV_8U);
    cv::Mat espectro_mag;
    cv::split(espectro_vazio, planes);
    cv::magnitude(planes[0], planes[1], espectro_mag);
    cv::normalize(espectro_mag, espectro_vis, 0, 255, cv::NORM_MINMAX, CV_8U);

    // Destaque visual do ponto
    cv::Mat espectro_color;
    cv::cvtColor(espectro_vis, espectro_color, cv::COLOR_GRAY2BGR);
    cv::circle(espectro_color, cv::Point(N_grid/2 - u0, N_grid/2 - v0), 2, cv::Scalar(0, 0, 255), -1);

    mm::show(std::vector<mm::Image>{mm::Image(espectro_color), mm::Image(onda_vis)},
             MM_OUT,
             {"Espectro (1 ponto ativo)", "Onda 2D Resultante (IDFT)"},
             2);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(espectro_color, "tmp/fig_05_grade_2d_0.png");
mm::write(onda_vis, "tmp/fig_05_grade_2d_1.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_grade_2d.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_grade_2d.cpp -o tmp/fig_05_grade_2d -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_grade_2d \
  && test -f "tmp/fig_05_grade_2d.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_grade_2d.png"

[1] Espectro (1 ponto ativo)
[2] Onda 2D Resultante (IDFT)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_grade_2d_0.png"),
            mm.read("tmp/fig_05_grade_2d_1.png"),
        ],
        titles=[
            'Espectro (1 ponto ativo)',
            'Onda 2D Resultante (IDFT)',
        ],
        cols=2,
        figsize=(10, 4),
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_grade_2d_0.png (ver a versao Python)")

Figura 5.3: Toda frequência no espectro (ponto isolado) corresponde a uma onda senoidal 2D rotacionada no domínio espacial.

∿ Simulador: Síntesis de Frecuencia 2D (IDFT) Espacio de Fourier

Frecuencia u

10

Frecuencia v

5

Distancia R

11.18

Ángulo θ

26.6°

Espectro (Haz clic para mover el punto)

➔

Onda 2D Resultante (Dominio Espacial)

Desplazamiento u (X): 10

Desplazamiento v (Y): 5

Figura 5.4: Simulador interactivo de la síntesis de Fourier 2D. Cambie la posición horizontal (\(u\)) y vertical (\(v\)) del coeficiente en el espectro de frecuencias centrado y observe cómo la distancia respecto al centro determina la frecuencia espacial (grosor) y el ángulo determina la orientación de la onda sinusoidal generada.

5.3.4 Definición Matemática

Considere una imagen \(f(x,y)\) con dimensiones \(M \times N\). Su Transformada Discreta de Fourier 2D (DFT) se define por:

\[ F(u,v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)\, e^{-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)} \tag{5.1}\]

donde \(u = 0, 1, \ldots, M-1\) y \(v = 0, 1, \ldots, N-1\) representan las frecuencias discretas en las direcciones horizontal y vertical, respectivamente. El término exponencial corresponde a una senoide bidimensional, cuya frecuencia y orientación están determinadas por los índices \((u,v)\).

La Transformada Discreta Inversa de Fourier 2D (IDFT) reconstruye la imagen original a partir de sus coeficientes:

\[ f(x,y) = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u,v)\, e^{j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)} \tag{5.2}\]

Las Ecuaciones Ecuación 5.1 y Ecuación 5.2 muestran que la DFT y la IDFT forman un par de transformaciones: la primera convierte la imagen al dominio de la frecuencia, mientras que la segunda reconstruye exactamente la imagen original a partir de sus coeficientes.

NotaSobre el símbolo \(j\)

El término \(j\) denota la unidad imaginaria, definida por \(j^2 = -1\). En ingeniería y procesamiento de señales, se adopta \(j\) en lugar de \(i\) para evitar conflictos con la notación de corriente eléctrica. Su uso en la exponencial compleja, regida por la fórmula de Euler (\(e^{j\theta} = \cos\theta + j\sin\theta\)), permite representar de forma compacta la amplitud y la fase de cada frecuencia espacial presente en la imagen.

Nota¿Qué es el componente DC?

El coeficiente \(F(0,0)\), denominado componente DC (Direct Current), es igual a la suma de las intensidades de todos los píxeles de la imagen (ver Figura 5.5):

\[ F(0,0)=MN\,\bar{f}, \]

donde \(\bar{f}\) es la intensidad media de la imagen. Por ello, el componente DC representa el nivel medio de intensidad y, en la mayoría de las imágenes naturales, posee la mayor magnitud del espectro.

Los demás coeficientes representan variaciones en torno a esa media. Tras la aplicación del FFT Shift, el componente DC se desplaza al centro del espectro, concentrando las bajas frecuencias en la región central y las altas frecuencias en los bordes.

Anatomía del Espectro de Fourier 2D (tras fftshift)

Figura 5.5: Diagrama conceptual del espectro de Fourier 2D centrado.

5.3.5 Magnitud y Fase

Cada coeficiente de la Transformada Discreta de Fourier (DFT) es un número complejo y puede escribirse como

\[ F(u,v)=R(u,v)+j\,I(u,v), \]

donde \(R(u,v)\) e \(I(u,v)\) corresponden, respectivamente, a las partes real e imaginaria del coeficiente. De la Ecuación Ecuación 5.1, se obtienen

\[ R(u,v)= \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y) \cos\!\left( 2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right) \right), \]

y

\[ I(u,v)= - \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y) \sin\!\left( 2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right) \right). \]

A partir de esta representación, se definen dos magnitudes fundamentales:

• Magnitud, que indica la intensidad de la componente de frecuencia,

\[ |F(u,v)|=\sqrt{R(u,v)^2+I(u,v)^2}; \]

• Fase, que determina la alineación (o desplazamiento) espacial de la componente,

\[ \phi(u,v)=\operatorname{atan2}\!\left(I(u,v),\,R(u,v)\right). \]

Así, cada coeficiente también puede escribirse en su forma polar,

\[ F(u,v)=|F(u,v)|\,e^{j\phi(u,v)}. \]

El espectro de Fourier puede, por tanto, visualizarse mediante dos imágenes distintas: el espectro de magnitud, normalmente utilizado para analizar la distribución de las frecuencias, y el espectro de fase, que describe la organización espacial de las componentes senoidales.

Aunque el espectro de magnitud sea el más utilizado para la inspección visual, la fase contiene gran parte de la información estructural de la imagen. La combinación de magnitud y fase permite reconstruir exactamente la imagen original mediante la IDFT.

5.3.6 ¿Qué transportan la magnitud y la fase?

Una demostración clásica consiste en combinar la magnitud de una imagen con la fase de otra y reconstruir el resultado. Este experimento evidencia que:

• La fase preserva la estructura espacial de la imagen, incluida la posición de los objetos, sus contornos y su geometría. Pequeñas alteraciones en la fase pueden provocar grandes cambios visuales.
• La magnitud controla cómo se distribuye la energía entre las frecuencias espaciales, influyendo principalmente en el contraste y la textura.

Cuando una imagen se reconstruye con la magnitud de A y la fase de B, el resultado tiende a parecerse más a B que a A, evidenciando que la fase es el principal componente responsable de la organización espacial de la escena. Sin embargo, la magnitud sigue siendo importante, ya que modula el contraste de las estructuras reconstruidas. Así, una reconstrucción fiel depende de la combinación coherente entre magnitud y fase.

Un ejemplo de este comportamiento se presenta en la Figura 5.6.

NotaAnalogía con el Audio: Limitaciones y Precauciones

La fase de una señal desempeña papeles distintos en audio e imágenes:

• Audio estéreo o multicanal: la fase relativa entre los canales es fundamental para la percepción de la posición de las fuentes sonoras, mediante las diferencias interaurales de tiempo (ITD, Interaural Time Differences).
• Audio monoaural: la fase absoluta ejerce poca influencia perceptual directa.
• Imágenes (DFT): la fase es el principal factor responsable de la organización espacial de la escena, mientras que la magnitud modula el contraste y la distribución de la energía entre las frecuencias.

En ambos dominios, la magnitud está relacionada con la intensidad de los componentes de frecuencia: en audio, influye en el timbre y la intensidad percibida; en imágenes, influye en el contraste y la textura.

%%writefile tmp/fig_05_dft_intro.cpp
#define MM_OUT "tmp/fig_05_dft_intro.png"
// Compile with: g++ -std=c++17 -O2 program.cpp -o program $(pkg-config --cflags --libs opencv4)
#include <opencv2/opencv.hpp>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <iostream>
#include <string>
#include <cmath>
#include <vector>
#include <filesystem>
#include "morph.hpp"

// Helper: swap quadrants for fftshift
void fftshift(cv::Mat& mat) {
    int cx = mat.cols / 2;
    int cy = mat.rows / 2;
    cv::Mat q0(mat, cv::Rect(0, 0, cx, cy));
    cv::Mat q1(mat, cv::Rect(cx, 0, cx, cy));
    cv::Mat q2(mat, cv::Rect(0, cy, cx, cy));
    cv::Mat q3(mat, cv::Rect(cx, cy, cx, cy));
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}

// Helper: compute magnitude spectrum in dB
cv::Mat spectrumMag(cv::Mat& img) {
    cv::Mat img32f;
    img.convertTo(img32f, CV_32F);
    cv::Mat planes[] = {img32f, cv::Mat::zeros(img32f.size(), CV_32F)};
    cv::Mat complexImg;
    cv::merge(planes, 2, complexImg);
    cv::dft(complexImg, complexImg);
    cv::Mat mag, phase;
    cv::magnitude(complexImg, complexImg, mag);

    // Shift quadrants
    fftshift(mag);

    // Compute log(1+|F|) and normalize
    mag += 1.0;
    cv::log(mag, mag);
    cv::normalize(mag, mag, 0, 255, cv::NORM_MINMAX, CV_8U);
    return mag;
}

int main() {
    // ── Experimento: A Importância da Fase ───────────────────────────────────────
    // ── Carregamento da imagem ────────────────────────────────────────────────────
    std::string url = "https://upload.wikimedia.org/wikipedia/commons/2/25/GAZI.MD.AHAD_11.jpg";
    std::string caminho = "imagens/coins.jpg";

    if (!std::filesystem::exists(caminho)) {
        std::filesystem::create_directories("imagens");
        cv::Mat img_obj = cv::imread(url, cv::IMREAD_COLOR);
        mm::write(img_obj, caminho);
    }

    cv::Mat img_color = cv::imread(caminho, cv::IMREAD_COLOR);
    cv::Mat img_gray_mat;
    cv::cvtColor(img_color, img_gray_mat, cv::COLOR_BGR2GRAY);
    mm::Image img_gray(img_gray_mat);

    cv::Mat img_a;
    cv::resize(img_gray_mat, img_a, cv::Size(400, 400));

    // Criar uma imagem B sintética (padrão geométrico)
    cv::Mat img_b = cv::Mat::zeros(400, 400, CV_8U);
    cv::rectangle(img_b, cv::Point(100, 100), cv::Point(300, 300), cv::Scalar(255), -1);
    cv::circle(img_b, cv::Point(200, 200), 150, cv::Scalar(128), 10);

    // FFT de A e B
    cv::Mat img_a_f;
    img_a.convertTo(img_a_f, CV_32F);
    cv::Mat planesA[] = {img_a_f, cv::Mat::zeros(img_a.size(), CV_32F)};
    cv::Mat complexA;
    cv::merge(planesA, 2, complexA);
    cv::dft(complexA, complexA);

    cv::Mat img_b_f;
    img_b.convertTo(img_b_f, CV_32F);
    cv::Mat planesB[] = {img_b_f, cv::Mat::zeros(img_b.size(), CV_32F)};
    cv::Mat complexB;
    cv::merge(planesB, 2, complexB);
    cv::dft(complexB, complexB);

    // Separar magnitude e fase
    cv::Mat magA[2], magB[2], phaseA[2], phaseB[2];
    cv::Mat planesA_out[2], planesB_out[2];
    cv::split(complexA, planesA_out);
    cv::split(complexB, planesB_out);
    cv::magnitude(planesA_out[0], planesA_out[1], magA[0]);
    cv::magnitude(planesB_out[0], planesB_out[1], magB[0]);
    cv::phase(planesA_out[0], planesA_out[1], phaseA[0]);
    cv::phase(planesB_out[0], planesB_out[1], phaseB[0]);

    // Troca de fase: reconstruir com magnitude de A e fase de B
    cv::Mat complex_recAB[2];
    cv::Mat realAB, imagAB;
    cv::polarToCart(magA[0], phaseB[0], realAB, imagAB);
    complex_recAB[0] = realAB;
    complex_recAB[1] = imagAB;
    cv::Mat rec_complexAB;
    cv::merge(complex_recAB, 2, rec_complexAB);
    cv::Mat rec_A_mag_B_fase_mat;
    cv::idft(rec_complexAB, rec_A_mag_B_fase_mat, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT);

    // Troca de fase: reconstruir com magnitude de B e fase de A
    cv::Mat complex_recBA[2];
    cv::Mat realBA, imagBA;
    cv::polarToCart(magB[0], phaseA[0], realBA, imagBA);
    complex_recBA[0] = realBA;
    complex_recBA[1] = imagBA;
    cv::Mat rec_complexBA;
    cv::merge(complex_recBA, 2, rec_complexBA);
    cv::Mat rec_B_mag_A_fase_mat;
    cv::idft(rec_complexBA, rec_B_mag_A_fase_mat, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT);

    // Converter para mm::Image
    cv::Mat rec_A_mag_B_fase_8u, rec_B_mag_A_fase_8u;
    rec_A_mag_B_fase_mat.convertTo(rec_A_mag_B_fase_8u, CV_8U);
    rec_B_mag_A_fase_mat.convertTo(rec_B_mag_A_fase_8u, CV_8U);

    mm::Image img_a_im(img_a);
    mm::Image img_b_im(img_b);
    mm::Image rec_A_mag_B_fase(rec_A_mag_B_fase_8u);
    mm::Image rec_B_mag_A_fase(rec_B_mag_A_fase_8u);

    mm::show({img_a_im, img_b_im, rec_A_mag_B_fase, rec_B_mag_A_fase},
             MM_OUT, {"Imagem A", "Imagem B", "Mag(A) + Fase(B)", "Mag(B) + Fase(A)"}, 4);

    std::cout << "💡 La fase preserva bordes y contornos; la magnitud controla contraste y" << std::endl;
    std::cout << "textura. En audio estéreo, la fase afecta la localización espacial; en" << std::endl;
    std::cout << "imágenes, determina la organización de la escena." << std::endl;

    // Variables to persist
    img_gray = img_gray;  // already set

    
// [pdi:state-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp/state");
mm::write(img_gray, "tmp/state/img_gray_20.png");
// [pdi:state-io:end]

// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(img_a, "tmp/fig_05_dft_intro_0.png");
mm::write(img_b, "tmp/fig_05_dft_intro_1.png");
mm::write(rec_A_mag_B_fase, "tmp/fig_05_dft_intro_2.png");
mm::write(rec_B_mag_A_fase, "tmp/fig_05_dft_intro_3.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_dft_intro.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dft_intro.cpp -o tmp/fig_05_dft_intro -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dft_intro \
  && test -f "tmp/fig_05_dft_intro.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dft_intro.png"

[1] Imagem A
[2] Imagem B
[3] Mag(A) + Fase(B)
[4] Mag(B) + Fase(A)
💡 La fase preserva bordes y contornos; la magnitud controla contraste y
textura. En audio estéreo, la fase afecta la localización espacial; en
imágenes, determina la organización de la escena.

try:
    mm.show(
        [
            mm.read("tmp/fig_05_dft_intro_0.png"),
            mm.read("tmp/fig_05_dft_intro_1.png"),
            mm.read("tmp/fig_05_dft_intro_2.png"),
            mm.read("tmp/fig_05_dft_intro_3.png"),
        ],
        titles=[
            'Imagem A',
            'Imagem B',
            'Mag(A) + Fase(B)',
            'Mag(B) + Fase(A)',
        ],
        cols=4,
        figsize=(16, 4),
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dft_intro_0.png (ver a versao Python)")

Figura 5.6: Experimento de troca de fase: Imagem A (moedas) e Imagem B (padrão geométrico) reconstruídas com magnitudes e fases trocadas. O resultado mostra que a estrutura visual é muito mais sensível à fase do que à magnitude: quando a fase de B é mantida, a imagem resultante preserva a organização espacial de B, mesmo com a magnitude de A. A magnitude, por sua vez, influencia principalmente o contraste e a textura. Observe que a qualidade da reconstrução não é perfeita — há artefatos visíveis —, evidenciando a interdependência entre fase e magnitude para uma representação fiel da imagem.

5.4 Teorema de la Convolución y Estrategias de Filtrado

El Teorema de la Convolución establece una relación fundamental entre los dominios espacial y de la frecuencia:

\[ f(x,y) \circledast h(x,y) \;\overset{\mathcal{F}}{\longleftrightarrow}\; F(u,v)\,H(u,v) \tag{5.3}\]

donde \(\circledast\) representa la convolución circular discreta. Así, la convolución entre una imagen \(f(x,y)\) y un filtro \(h(x,y)\) puede sustituirse por la multiplicación de sus espectros.

En la práctica, para obtener el mismo resultado de la convolución lineal realizada en el dominio espacial, se aplica relleno con ceros (zero-padding) antes de la Transformada Rápida de Fourier (FFT), evitando artefactos en los bordes de la imagen.

Sin embargo, no siempre el filtrado en el dominio de la frecuencia es la alternativa más eficiente. Para filtros como el gaussiano y el filtro de la media (Box Filter), la propiedad de separabilidad permite reducir significativamente el costo computacional de la convolución en el dominio espacial.

5.4.1 Kernel Separable vs. No separable

Un kernel separable puede escribirse como el producto externo de dos vectores unidimensionales,

\[ H = v\,h^T, \]

lo que permite que la convolución bidimensional se reemplace por dos convoluciones unidimensionales consecutivas: una en la dirección horizontal y otra en la vertical.

En cambio, un kernel no separable no admite esta descomposición y, por lo tanto, su convolución debe realizarse directamente sobre la vecindad bidimensional.

En la práctica, para un kernel de dimensión \(K \times K\), la convolución directa requiere \(K^2\) multiplicaciones por píxel, mientras que un kernel separable requiere solo \(2K\) multiplicaciones, reduciendo significativamente el costo computacional.

5.4.2 Análisis de Eficiencia Computacional

Considere una imagen de dimensiones \(M \times N\) y un filtro cuadrado de tamaño \(K \times K\). La Tabla 5.2 compara la complejidad de las principales estrategias de filtrado.

Tabla 5.2: Comparación de la complejidad de la convolución directa, separable y vía Transformada Rápida de Fourier (FFT).

Método de Filtrado	Complejidad Asintótica	Dependencia de \(K\)	Aplicación típica
Espacial no separable	\(\mathcal{O}(MNK^2)\)	Cuadrática	Kernels pequeños y no separables
Espacial separable	\(\mathcal{O}(MNK)\)	Lineal	Filtros Gaussiano y de la media
Vía FFT	\(\mathcal{O}(MN\log(MN))\)	Independiente de \(K\)	Kernels grandes

Para kernels pequeños, la convolución espacial, especialmente cuando el filtro es separable, suele ser más eficiente debido al bajo costo de las operaciones. A medida que el tamaño del kernel aumenta, el filtrado vía FFT se vuelve más ventajoso, ya que su costo prácticamente no depende de la dimensión del filtro.

5.4.3 Discusión de los resultados experimentales

El gráfico obtenido en el ensayo con la imagen de las monedas (\(2560 \times 1920\)), presentado en la Figura 5.7, confirma el comportamiento previsto por el análisis de complejidad computacional.

1. Convolución no separable (\(\mathcal{O}(MNK^2)\))
   La convolución directa presenta un crecimiento cuadrático con el tamaño del kernel. Para valores pequeños de \(K\), el costo es bajo, pero aumenta rápidamente a medida que el kernel crece, volviéndose inviable para aplicaciones en tiempo real.

2. Filtrado mediante FFT (\(\mathcal{O}(MN \log(MN))\))
   El costo de la FFT depende únicamente del tamaño de la imagen, siendo independiente de \(K\). Por ello, su rendimiento permanece aproximadamente constante al variar el kernel, lo que la hace ventajosa para filtros grandes o no separables.

3. Convolución separable (\(\mathcal{O}(MNK)\))
   La descomposición del kernel en dos filtros unidimensionales reduce significativamente el costo computacional. En la práctica, este enfoque tiende a ser el más eficiente para filtros separables, especialmente en implementaciones optimizadas.

En general, la elección del método depende del tamaño y la estructura del kernel. Los filtros separables son más eficientes en el dominio espacial, mientras que la FFT resulta más ventajosa para kernels grandes o múltiples convoluciones en el dominio de la frecuencia.

\[ g = \mathcal{F}^{-1}\bigl[\mathcal{F}(f)\cdot \mathcal{F}(h)\bigr] \quad \text{(FFT)} \qquad g = f \circledast h \quad \text{(convolución directa)} \qquad g = (f \circledast v) \circledast h^T \quad \text{(separable)} \tag{5.4}\]

donde:

• \(f(x,y)\) representa la imagen de entrada;
• \(h(x,y)\) es el kernel bidimensional del filtro;
• \(v\) y \(h^T\) son, respectivamente, los vectores vertical y horizontal que componen el kernel separable.

%%writefile tmp/fig_05_conv_eficiencia.cpp
#define MM_OUT "tmp/fig_05_conv_eficiencia.png"
//| label: fig-05-conv-eficiencia
//| fig-cap: "Comparação de eficiência: Convolução Não Separável (Espacial 2D), Separável (Espacial 1D) e via FFT."
//| echo: false
//| output: true

#include <vector>
#include <string>
#include <cmath>
#include "morph.hpp"
#include <filesystem>

int main() {
    // Trilha C++: curvas de CUSTO relativo (contagem de operacoes) — a forma das
    // curvas (K^2 vs 2K vs log N) e o que importa; a trilha py mede tempos reais.
    std::vector<double> K = {3, 7, 11, 15, 21, 31, 41, 51};
    double N2 = 256.0 * 256.0;
    std::vector<double> t_nao_sep, t_sep, t_fft;
    for (double k : K) {
        t_nao_sep.push_back(N2 * k * k / 1e6);
        t_sep.push_back(N2 * 2.0 * k / 1e6);
        t_fft.push_back(N2 * std::log2(N2) / 1e6);
    }

    std::vector<std::vector<double>> xs = {K, K, K};
    std::vector<std::vector<double>> ys = {t_nao_sep, t_sep, t_fft};
    std::vector<std::string> labels = {
        "Nao Separavel  O(N^2 K^2)",
        "Separavel  O(N^2 . 2K)",
        "Via FFT  O(N^2 log N)"
    };

    mm::Image chart = mm::lineChart(
        xs, ys, labels, {},
        "Custo relativo: Espacial vs Frequencia",
        "Tamanho do kernel  K", "Operacoes  (x10^6)"
    );
    mm::show({chart}, MM_OUT, {"Comparacao de complexidade"}, 1);
    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(chart, "tmp/fig_05_conv_eficiencia_0.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_conv_eficiencia.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_conv_eficiencia.cpp -o tmp/fig_05_conv_eficiencia -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_conv_eficiencia \
  && test -f "tmp/fig_05_conv_eficiencia.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_conv_eficiencia.png"

[1] Comparacao de complexidade

try:
    mm.show(
        [
            mm.read("tmp/fig_05_conv_eficiencia_0.png"),
        ],
        titles=[
            'Comparacao de complexidade',
        ],
        cols=1,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_conv_eficiencia_0.png (ver a versao Python)")

Figura 5.7: Comparação de eficiência: Convolução Não Separável (Espacial 2D), Separável (Espacial 1D) e via FFT.

ImportanteEl problema de la convolución circular (wrap-around)

La Transformada Discreta de Fourier (DFT) asume que la imagen se extiende periódicamente en el espacio, es decir, que sus bordes se repiten indefinidamente.

En esta condición, la multiplicación en el dominio de la frecuencia corresponde a una convolución circular en el dominio espacial. Como consecuencia, regiones opuestas de la imagen (parte superior e inferior, izquierda y derecha) interactúan artificialmente, tal como se ilustra en Figura 5.8.

La aplicación de zero-padding antes de la FFT reduce este efecto al extender la imagen con valores nulos en los bordes, aproximando el resultado a la convolución lineal. Este comportamiento puede interpretarse a la luz del Teorema de la Convolución, presentado en Figura 5.9.

%%writefile tmp/fig_05_padding_error.cpp
#define MM_OUT "tmp/fig_05_padding_error.png"
//| label: fig-05-padding-error
//| fig-cap: "Sem *padding*, um deslocamento severo faz a imagem vazar para o lado oposto (convolução circular)."
//| echo: true
//| output: true
#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include <complex>
#include "morph.hpp"
#include <filesystem>


// Helpers para fftshift manual
static void fftshift_manual(const cv::Mat& src, cv::Mat& dst) {
    int cx = src.cols / 2;
    int cy = src.rows / 2;
    cv::Mat q0(src, cv::Rect(0, 0, cx, cy));
    cv::Mat q1(src, cv::Rect(cx, 0, src.cols - cx, cy));
    cv::Mat q2(src, cv::Rect(0, cy, cx, src.rows - cy));
    cv::Mat q3(src, cv::Rect(cx, cy, src.cols - cx, src.rows - cy));
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
    dst = src; // na prática, com q0..q3 cópias, pode precisar de clone
}

int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // img_gray já é válida (injetada automaticamente)

    // ── Definir M e N ─────────────────────────────────────────────────────────────
    int M = img_gray.h;
    int N = img_gray.w;

    // Simulação de um filtro de deslocamento brutal
    // H_shift[u,v] = exp(-1j * 2*pi * (u*delta_x/M + v*delta_y/N))
    cv::Mat img_gray_float;
    mm::Image img_gray_mm = img_gray;
    cv::Mat img_gray_mat = img_gray_mm;
    img_gray_mat.convertTo(img_gray_float, CV_64F);

    // FFT da imagem
    cv::Mat F_img_complex;
    cv::dft(img_gray_float, F_img_complex, cv::DFT_COMPLEX_OUTPUT);

    // Criar H_shift complexo (centrado no centro da imagem, como em numpy)
    // Em numpy, o loop i,j percorre a imagem sem fftshift — então faremos igual
    int delta = 120;
    cv::Mat H_shift_real(M, N, CV_64F, cv::Scalar(0));
    cv::Mat H_shift_imag(M, N, CV_64F, cv::Scalar(0));

    for (int u = 0; u < M; ++u) {
        for (int v = 0; v < N; ++v) {
            double theta = -2.0 * M_PI * (static_cast<double>(u) * delta / M + static_cast<double>(v) * delta / N);
            H_shift_real.at<double>(u, v) = std::cos(theta);
            H_shift_imag.at<double>(u, v) = std::sin(theta);
        }
    }

    // Multiplicação no domínio da frequência: F_img * H_shift
    cv::Mat F_img_planes[2];
    cv::split(F_img_complex, F_img_planes);
    cv::Mat mult_real = F_img_planes[0].mul(H_shift_real) - F_img_planes[1].mul(H_shift_imag);
    cv::Mat mult_imag = F_img_planes[0].mul(H_shift_imag) + F_img_planes[1].mul(H_shift_real);

    std::vector<cv::Mat> planes = {mult_real, mult_imag};
    cv::Mat F_filtered;
    cv::merge(planes, F_filtered);

    // IFFT
    cv::Mat img_vazada_float;
    cv::idft(F_filtered, img_vazada_float, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT);

    // Normalizar para visualização
    cv::Mat img_vazada_vis;
    cv::normalize(img_vazada_float, img_vazada_vis, 0, 255, cv::NORM_MINMAX, CV_8U);

    // Converter para mm::Image para mostrar
    mm::Image img_vazada_vis_img(img_vazada_vis);

    // Mostrar
    mm::show({img_gray, img_vazada_vis_img},
             MM_OUT,
             {"Original", "Filtragem s/ Padding (Vazamento)"}, 2);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(img_gray, "tmp/fig_05_padding_error_0.png");
mm::write(img_vazada_vis, "tmp/fig_05_padding_error_1.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_padding_error.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_padding_error.cpp -o tmp/fig_05_padding_error -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_padding_error \
  && test -f "tmp/fig_05_padding_error.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_padding_error.png"

[1] Original
[2] Filtragem s/ Padding (Vazamento)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_padding_error_0.png"),
            mm.read("tmp/fig_05_padding_error_1.png"),
        ],
        titles=[
            'Original',
            'Filtragem s/ Padding (Vazamento)',
        ],
        cols=2,
        figsize=(10, 4),
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_padding_error_0.png (ver a versao Python)")

Figura 5.8: Sem padding, um deslocamento severo faz a imagem vazar para o lado oposto (convolução circular).

%%writefile tmp/fig_05_conv_teorema.cpp
#define MM_OUT "tmp/fig_05_conv_teorema.png"
//| label: fig-05-conv-teorema
//| fig-cap: "Teorema da Convolución: filtrar en el dominio del espacio (Gaussiana) equivale a multiplicar el espectro por $H(u,v)$ en la frecuencia. Las salidas coinciden visualmente."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include "morph.hpp"
#include <vector>
#include <string>
#include <filesystem>

int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // Misma suavización por dos caminos: espacio vs. frecuencia.
    int M = img_gray.h;
    int N = img_gray.w;
    cv::Mat H = mm::gaussFilter(M, N, 30);      // H(u,v): transferencia pasa-bajas gaussiana
    mm::Image f_freq = mm::freqFilter(img_gray, H);    // convolución vía multiplicación en frecuencia
    mm::Image f_esp = mm::gaussian(img_gray, 31, 5);  // misma suavización hecha en el espacio

    mm::show(
        std::vector<mm::Image>{img_gray, f_esp, f_freq},
        MM_OUT,
        std::vector<std::string>{"Original", "Espacio: Gaussiana", "Frecuencia: H(u,v).F(u,v)"},
        3
    );
    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(img_gray, "tmp/fig_05_conv_teorema_0.png");
mm::write(f_esp, "tmp/fig_05_conv_teorema_1.png");
mm::write(f_freq, "tmp/fig_05_conv_teorema_2.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_conv_teorema.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_conv_teorema.cpp -o tmp/fig_05_conv_teorema -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_conv_teorema \
  && test -f "tmp/fig_05_conv_teorema.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_conv_teorema.png"

[1] Original
[2] Espacio: Gaussiana
[3] Frecuencia: H(u,v).F(u,v)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_conv_teorema_0.png"),
            mm.read("tmp/fig_05_conv_teorema_1.png"),
            mm.read("tmp/fig_05_conv_teorema_2.png"),
        ],
        titles=[
            'Original',
            'Espaco: Gaussiana',
            'Frequencia: H(u,v).F(u,v)',
        ],
        cols=3,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_conv_teorema_0.png (ver a versao Python)")

Figura 5.9: Teorema da Convolução: filtrar no domínio do espaço (Gaussiana) equivale a multiplicar o espectro por \(H(u,v)\) na frequência. As saídas coincidem visualmente.

NotaSobre la diferencia numérica

La diferencia residual del orden de \(10^{-13}\) no viola el Teorema de la Convolución, sino que refleja limitaciones computacionales inherentes a la aritmética de coma flotante (doble precisión, ~\(10^{-16}\)) y al orden de las operaciones entre los dos métodos:

• Convolución espacial: suma ponderada de vecinos con redondeos sucesivos.
• Convolución en frecuencia: involucra tres transformadas FFT y una multiplicación compleja, sujeta a errores de truncamiento y cuantización.

Por lo tanto, la igualdad teórica es exacta, pero la implementación numérica produce una diferencia prácticamente nula (error relativo < \(10^{-12}\)), lo que confirma el teorema dentro de la precisión de la máquina.

5.5 Filtros en el dominio de la frecuencia

Un filtro en el dominio de la frecuencia puede interpretarse como una función de transferencia aplicada al espectro de la imagen. En esta representación, cada coeficiente de frecuencia se multiplica por un valor entre 0 y 1, que determina su atenuación o preservación. La forma de esta función define el efecto visual del filtro.

Corte abrupto y ringing. Los filtros ideales con transición instantánea en una frecuencia de corte \(D_0\) producen discontinuidades en el dominio de la frecuencia. Esta discontinuidad se refleja en el dominio espacial como oscilaciones cercanas a los bordes, conocidas como ringing. Este efecto está asociado a la convolución con funciones de soporte infinito en el espacio, como la función sinc, tal como se ilustra en la Figura 5.10.

Filtros con transición suave. Alternativas como los filtros gaussiano y Butterworth suavizan la transición entre las regiones de paso y de rechazo, reduciendo el ringing. En contrapartida, esta suavización implica una frontera de separación menos definida entre las frecuencias preservadas y las atenuadas.

%%writefile tmp/fig_05_conv_teorema_zoom.cpp
#define MM_OUT "tmp/fig_05_conv_teorema_zoom.png"
#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>


int main() {
    //| label: fig-05-conv-teorema-zoom
    //| fig-cap: "A Dualidade Perigosa: o corte abrupto na Frequência (cilindro Ideal) vira obrigatoriamente uma *sinc* no espaço. Suas ondulações causam o *ringing* fantasma nas bordas da imagem."
    //| echo: true
    //| output: true

    // Filtro Ideal na frequência (cilindro) e sua resposta espacial (sinc 2D).
    int N = 128;
    cv::Mat H_freq = mm::idealFilter(N, N, 20);      // 1 dentro do raio 20, 0 fora
    mm::Image h_space = mm::spatialKernel(H_freq);   // ifft2(H) -> ondulações da sinc (ringing)

    mm::show(
        std::vector<mm::Image>{H_freq, h_space},
        MM_OUT,
        std::vector<std::string>{
            "Frequencia: filtro Ideal (cilindro)",
            "Espaco: ondulacoes da sinc (causa do ringing)"
        },
        2
    );

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(H_freq, "tmp/fig_05_conv_teorema_zoom_0.png");
mm::write(h_space, "tmp/fig_05_conv_teorema_zoom_1.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_conv_teorema_zoom.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_conv_teorema_zoom.cpp -o tmp/fig_05_conv_teorema_zoom -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_conv_teorema_zoom \
  && test -f "tmp/fig_05_conv_teorema_zoom.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_conv_teorema_zoom.png"

[1] Frequencia: filtro Ideal (cilindro)
[2] Espaco: ondulacoes da sinc (causa do ringing)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_conv_teorema_zoom_0.png"),
            mm.read("tmp/fig_05_conv_teorema_zoom_1.png"),
        ],
        titles=[
            'Frequencia: filtro Ideal (cilindro)',
            'Espaco: ondulacoes da sinc (causa do ringing)',
        ],
        cols=2,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_conv_teorema_zoom_0.png (ver a versao Python)")

Figura 5.10: A Dualidade Perigosa: o corte abrupto na Frequência (cilindro Ideal) vira obrigatoriamente uma sinc no espaço. Suas ondulações causam o ringing fantasma nas bordas da imagem.

5.5.1 Filtros Pasa-Bajos

Los filtros pasa-bajos atenúan los componentes de alta frecuencia, lo que resulta en un suavizado de la imagen y una reducción del ruido. Tras la centralización del espectro (FFT Shift), la distancia de cada punto al centro viene dada por:

\[ D(u,v) = \sqrt{\left(u - \tfrac{M}{2}\right)^2 + \left(v - \tfrac{N}{2}\right)^2} \tag{5.5}\]

Filtro Ideal (LPFI): \[ H_{\text{ideal}}(u,v) = \begin{cases} 1, & D(u,v) \leq D_0 \\ 0, & D(u,v) > D_0 \end{cases} \tag{5.6}\]

El corte abrupto en \(D_0\) introduce discontinuidades en el dominio de la frecuencia, lo que da lugar a oscilaciones en el dominio espacial conocidas como ringing. Este efecto está asociado con la convolución con funciones de soporte infinito.

Filtro Gaussiano (LPFG): \[ H_{\text{gauss}}(u,v) = e^{-D^2(u,v)/(2\sigma^2)} \tag{5.7}\]

La suavidad de la función gaussiana en el dominio de la frecuencia evita discontinuidades, lo que elimina el ringing y produce una transición gradual entre las frecuencias preservadas y las atenuadas.

Filtro de Butterworth (LPFB) de orden \(n\): \[ H_{\text{BW}}(u,v) = \frac{1}{1 + \left[D(u,v)/D_0\right]^{2n}} \tag{5.8}\]

El parámetro \(n\) controla la suavidad de la transición entre el paso y el rechazo de frecuencias. Los valores pequeños producen transiciones suaves, mientras que los valores grandes aproximan el comportamiento del filtro ideal, con un mayor riesgo de ringing. En la Figura 5.11. se presenta un ejemplo comparativo.

Perfis de los Filtros Pasa-Baja — comparación visual (D₀ = 30)

A medida que aumenta la orden del Butterworth, el perfil se acerca al filtro Ideal — y el ringing aumenta.

Figura 5.11: Filtros pasa-baja.

5.5.2 Filtros Pasa-Alto y Pasa-Banda

Los filtros pasa-alto pueden obtenerse a partir de un filtro pasa-bajo complementario, definido como:

\[ H_{\text{HP}}(u,v) = 1 - H_{\text{LP}}(u,v) \]

Este tipo de filtro preserva componentes de alta frecuencia, resaltando bordes y detalles, mientras atenúa regiones de variación suave.

Los filtros pasa-banda preservan solo un rango intermedio de frecuencias, limitado por dos radios \(D_L\) y \(D_H\):

\[ H_{\text{BP}}(u,v) = H_{\text{LP}}^{(D_H)}(u,v)\cdot \left[1 - H_{\text{LP}}^{(D_L)}(u,v)\right] \]

Este tipo de filtrado es útil cuando se desea eliminar simultáneamente componentes de baja y alta frecuencia, preservando únicamente estructuras de escala intermedia.

Una aplicación importante es la eliminación de ruido periódico, en la cual patrones regulares aparecen como picos localizados en el espectro de magnitud. Estos picos pueden atenuarse mediante filtros notch (rechaza-banda), posicionados específicamente en las frecuencias no deseadas.

Ejemplos de filtros en el dominio de la frecuencia se presentan en el simulador de la Figura 5.12, Figura 5.13 y Figura 5.14..

🎛️ Simulador: Filtros en el Dominio de la Frecuencia Pasa-Baja / Pasa-Alta

Frecuencia de corte D₀ 30 Tipo de filtro

pasa-baja pasa-alta

Respuesta H(D)

Espectro filtrado |F · H|

Señal 1D — original vs filtrada

Energía retenida por banda (%)

Figura 5.12: Simulador interactivo de filtros en el dominio de la frecuencia.

%%writefile tmp/fig_05_filtros_freq.cpp
#define MM_OUT "tmp/fig_05_filtros_freq.png"
//| label: fig-05-filtros-freq
//| fig-cap: "Comparación entre filtros pasa-bajos: Ideal (D₀=30), Gaussiano (D₀=30) y Butterworth (D₀=30, n=2). Perfiles de H(u,v) a lo largo de una línea central e imágenes filtradas correspondientes."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include "morph.hpp"

// Helper function to visualize H(u,v) as an image
static cv::Mat H_vis(const cv::Mat& H) {
    cv::Mat scaled;
    cv::Mat H_uint8;
    // Scale H to 0-255 range and convert to uint8
    H.convertTo(H_uint8, CV_8U, 255.0);
    cv::normalize(H_uint8, scaled, 0, 255, cv::NORM_MINMAX);
    return scaled;
}

// Helper function to draw a profile curve
static void traca(const cv::Mat& row, const cv::Scalar& cor, cv::Mat& perfil, int M, int N, int Wp, int Hp, int linha) {
    cv::Point ant(-1, -1);
    for (int x = 0; x < N; x++) {
        int px = static_cast<int>(std::round(x * (Wp - 1) / (N - 1)));
        int py = static_cast<int>(std::round(10 + (1.0 - static_cast<float>(row.at<double>(0, x))) * (Hp - 20)));
        if (ant.x >= 0) {
            cv::line(perfil, ant, cv::Point(px, py), cor, 2, cv::LINE_AA);
        }
        ant = cv::Point(px, py);
    }
}

int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // Convert mm::Image to cv::Mat
    cv::Mat img_gray_cv = img_gray;

    int M = img_gray_cv.rows;
    int N = img_gray_cv.cols;
    double D0 = 30.0;
    int n_bw = 2;

    // ── Funciones de transferencia (H centrado) vía constructores de morph.hpp ───────
    //   Ideal:        1 si D<=D0, sino 0
    //   Gaussiano:    exp(-D^2/(2 D0^2))
    //   Butterworth:  1 / (1 + (D/D0)^(2n))
    cv::Mat H_ideal = mm::idealFilter(M, N, D0);
    cv::Mat H_gauss = mm::gaussFilter(M, N, D0);
    cv::Mat H_bw    = mm::butterFilter(M, N, D0, n_bw);

    // ── Imágenes filtradas vía FFT ────────────────────────────────────────────────
    mm::Image img_ideal = mm::freqFilter(img_gray, H_ideal);
    mm::Image img_gauss = mm::freqFilter(img_gray, H_gauss);
    mm::Image img_bw    = mm::freqFilter(img_gray, H_bw);

    // ── Perfiles de H(u,v) en la línea central, dibujados con cv::line ───────────────
    int linha = M / 2;
    int Wp = 512, Hp = 288;
    cv::Mat perfil(Hp, Wp, CV_8UC3, cv::Scalar(255, 255, 255));

    // Extraer la fila central de cada filtro H
    cv::Mat row_ideal = H_ideal.row(linha);
    cv::Mat row_gauss = H_gauss.row(linha);
    cv::Mat row_bw    = H_bw.row(linha);

    traca(row_ideal, cv::Scalar(48, 90, 216), perfil, M, N, Wp, Hp, linha);
    traca(row_gauss, cv::Scalar(117, 158, 29), perfil, M, N, Wp, Hp, linha);
    traca(row_bw,    cv::Scalar(183, 74, 83), perfil, M, N, Wp, Hp, linha);

    // Convert perfil from cv::Mat to mm::Image for display
    mm::Image perfil_img(perfil);

    mm::show(
        std::vector<mm::Image>{img_gray, img_ideal, img_gauss, img_bw,
                               H_vis(H_ideal), H_vis(H_gauss), H_vis(H_bw), perfil_img},
        MM_OUT,
        std::vector<std::string>{
            "Original", "LPF Ideal", "LPF Gaussiano", "LPF Butterworth (n=2)",
            "H Ideal", "H Gaussiano", "H Butterworth", "Perfiles H(u,v)",
        },
        4
    );

    return 0;
}

Overwriting tmp/fig_05_filtros_freq.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_filtros_freq.cpp -o tmp/fig_05_filtros_freq -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_filtros_freq \
  && test -f "tmp/fig_05_filtros_freq.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_filtros_freq.png"

[1] Original
[2] LPF Ideal
[3] LPF Gaussiano
[4] LPF Butterworth (n=2)
[5] H Ideal
[6] H Gaussiano
[7] H Butterworth
[8] Perfiles H(u,v)

try:
    mm.show(mm.read("tmp/fig_05_filtros_freq.png"), figsize=(16, 9))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_filtros_freq.png (ver a versao Python)")

Figura 5.13: Comparação entre filtros passa-baixa: Ideal (D₀=30), Gaussiano (D₀=30) e Butterworth (D₀=30, n=2). Perfis de H(u,v) ao longo de uma linha central e imagens filtradas correspondentes.

%%writefile tmp/fig_05_filtros_passa_alta.cpp
#define MM_OUT "tmp/fig_05_filtros_passa_alta.png"
//| label: fig-05-filtros-passa-alta
//| fig-cap: "Filtro passa-alta Gaussiano. (a) Original; (b) Filtro passa-alta (D₀=30) - as bordas das moedas e fundo texturizado são realçados."

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>

// Filtro passa-alta = complemento do passa-baixa Gaussiano (1 - H_gauss),
// construído com mm.gaussFilter(..., highpass=true) e aplicado via FFT.
int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // img_gray é fornecido automaticamente (mm::Image)

    int M = img_gray.h;
    int N = img_gray.w;
    double D0 = 30;
    cv::Mat H_alta = mm::gaussFilter(M, N, D0, true);
    mm::Image img_alta = mm::freqFilter(img_gray, H_alta);

    mm::show(
        std::vector<mm::Image>{img_gray, img_alta},
        MM_OUT,
        {"Original", "Passa-alta Gaussiano ($D_0=30$)"},
        2
    );

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(img_gray, "tmp/fig_05_filtros_passa_alta_0.png");
mm::write(img_alta, "tmp/fig_05_filtros_passa_alta_1.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_filtros_passa_alta.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_filtros_passa_alta.cpp -o tmp/fig_05_filtros_passa_alta -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_filtros_passa_alta \
  && test -f "tmp/fig_05_filtros_passa_alta.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_filtros_passa_alta.png"

[1] Original
[2] Passa-alta Gaussiano ($D_0=30$)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_filtros_passa_alta_0.png"),
            mm.read("tmp/fig_05_filtros_passa_alta_1.png"),
        ],
        titles=[
            'Original',
            'Passa-alta Gaussiano ($D_0=30$)',
        ],
        cols=2,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_filtros_passa_alta_0.png (ver a versao Python)")

Figura 5.14: Filtro passa-alta Gaussiano. (a) Original; (b) Filtro passa-alta (D₀=30) - as bordas das moedas e fundo texturizado são realçados.

5.5.3 Eliminación de Ruido Periódico

El ruido periódico — asociado a interferencias eléctricas, patrones regulares de sensores o artefactos de digitalización — aparece en el espectro de Fourier como picos puntuales simétricos alrededor del centro.

El filtro rechaza-banda (notch) atenúa selectivamente esas frecuencias, preservando las demás componentes de la imagen. Un ejemplo de aplicación se presenta en la Figura 5.15..

%%writefile tmp/fig_05_ruido_periodico.cpp
#define MM_OUT "tmp/fig_05_ruido_periodico.png"
//| label: fig-05-ruido-periodico
//| fig-cap: "Remoção de ruído periódico via filtro *notch* no domínio da frequência: (a) imagem com ruído senoidal, (b) espectro mostrando os picos do ruído, (c) máscara *notch* centrada nos picos, (d) imagem restaurada."
//| echo: true
//| output: true

// Ruído senoidal 2D -> espectro -> máscara notch nos 4 picos simétricos -> restauração.
#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include "morph.hpp"

int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // Variables pre-loaded: img_gray (mm::Image)

    int h_img = img_gray.h;
    int w_img = img_gray.w;

    // Criar meshgrid equivalente em 2D
    cv::Mat X(h_img, w_img, CV_64F), Y(h_img, w_img, CV_64F);
    for(int i = 0; i < h_img; i++) {
        for(int j = 0; j < w_img; j++) {
            X.at<double>(i,j) = j;
            Y.at<double>(i,j) = i;
        }
    }

    double u0 = 20, v0 = 20;  // frequências exatas do ruído
    cv::Mat ruido(h_img, w_img, CV_64F);
    for(int i = 0; i < h_img; i++) {
        for(int j = 0; j < w_img; j++) {
            ruido.at<double>(i,j) = 40 * std::sin(2 * M_PI * (u0 * j / w_img + v0 * i / h_img));
        }
    }

    // Converter imagem para double e adicionar ruído
    cv::Mat img_gray_cv = img_gray;  // conversão implícita
    cv::Mat img_float;
    img_gray_cv.convertTo(img_float, CV_64F);
    cv::Mat img_ruidosa_float = img_float + ruido;

    cv::Mat img_ruidosa_8u;
    cv::threshold(img_ruidosa_float, img_ruidosa_float, 0, 255, cv::THRESH_TRUNC);
    cv::threshold(img_ruidosa_float, img_ruidosa_float, 0, 255, cv::THRESH_TOZERO);
    img_ruidosa_float.convertTo(img_ruidosa_8u, CV_8U);
    mm::Image img_ruidosa(img_ruidosa_8u);

    // Espectro de magnitude (log) da imagem ruidosa
    mm::Image mag_vis = mm::spectrumMag(img_ruidosa);

    // Máscara notch: 1.0 em todo o plano, disco de 0.0 em cada um dos 4 picos
    cv::Mat mascara = cv::Mat::ones(h_img, w_img, CV_64F);
    int r_notch = 8;
    int cy = h_img / 2, cx = w_img / 2;

    // Gerar coordenadas de grade
    cv::Mat yy(h_img, w_img, CV_64F), xx(h_img, w_img, CV_64F);
    for(int i = 0; i < h_img; i++) {
        for(int j = 0; j < w_img; j++) {
            yy.at<double>(i,j) = i;
            xx.at<double>(i,j) = j;
        }
    }

    // Aplicar notches nos 4 picos
    int dys[] = {static_cast<int>(v0), -static_cast<int>(v0), static_cast<int>(v0), -static_cast<int>(v0)};
    int dxs[] = {static_cast<int>(u0), -static_cast<int>(u0), -static_cast<int>(u0), static_cast<int>(u0)};

    for(int k = 0; k < 4; k++) {
        int dy = dys[k], dx = dxs[k];
        cv::Mat dist;
        cv::Mat diff_y = yy - (cy + dy);
        cv::Mat diff_x = xx - (cx + dx);
        cv::magnitude(diff_x, diff_y, dist);

        for(int i = 0; i < h_img; i++) {
            for(int j = 0; j < w_img; j++) {
                if(dist.at<double>(i,j) <= r_notch) {
                    mascara.at<double>(i,j) = 0.0;
                }
            }
        }
    }

    cv::Mat mascara_8u;
    cv::Mat mascara_scaled = mascara * 255;
    mascara_scaled.convertTo(mascara_8u, CV_8U);
    mm::Image mascara_vis(mascara_8u);

    // Filtragem: aplica a máscara centrada como filtro no domínio da frequência
    mm::Image img_rest_vis = mm::freqFilter(img_ruidosa, mascara);

    // Calcular PSNR
    cv::Mat original_cv = img_gray;
    cv::Mat restaurada_cv = img_rest_vis;
    double psnr = cv::PSNR(original_cv, restaurada_cv);

    std::cout << "PSNR (original vs restaurada): " << psnr << " dB" << std::endl;

    // Construir títulos
    std::string titulo1 = "Com ruído periódico";
    std::string titulo2 = "Espectro (log)";
    std::string titulo3 = "Máscara notch";

    char buffer[100];
    snprintf(buffer, sizeof(buffer), "Restaurada (PSNR=%.1f dB)", psnr);
    std::string titulo4 = buffer;

    std::vector<mm::Image> imagens = {img_ruidosa, mag_vis, mascara_vis, img_rest_vis};
    std::vector<std::string> titulos = {titulo1, titulo2, titulo3, titulo4};

    mm::show(imagens, MM_OUT, titulos, 4);

    return 0;
}

Overwriting tmp/fig_05_ruido_periodico.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_ruido_periodico.cpp -o tmp/fig_05_ruido_periodico -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_ruido_periodico \
  && test -f "tmp/fig_05_ruido_periodico.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_ruido_periodico.png"

PSNR (original vs restaurada): 7.46158 dB
[1] Com ruído periódico
[2] Espectro (log)
[3] Máscara notch
[4] Restaurada (PSNR=7.5 dB)

try:
    mm.show(mm.read("tmp/fig_05_ruido_periodico.png"), figsize=(16, 4))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_ruido_periodico.png (ver a versao Python)")

Figura 5.15: Remoção de ruído periódico via filtro notch no domínio da frequência: (a) imagem com ruído senoidal, (b) espectro mostrando os picos do ruído, (c) máscara notch centrada nos picos, (d) imagem restaurada.

📌 Síntesis — Filtros Espectrales

Filtro	Efecto visual	Artefacto	Uso
Paso-bajo ideal	Suavizado intenso	Ringing	Ilustrativo
Paso-bajo gaussiano	Suavizado suave	No presenta ringing	Suavizado general
Paso-bajo Butterworth	Suavizado controlado	Ringing (órdenes altos)	Compromiso entre suavizado y selectividad
Paso-alto	Realce de bordes	Amplificación de ruido	Detección de contornos
Notch	Eliminación selectiva de frecuencias	Posibles distorsiones locales	Eliminación de ruido periódico

El diseño de filtros en el dominio de la frecuencia consiste en la definición de máscaras espectrales. Sin embargo, efectos en el dominio espacial, como ringing y desenfoque, emergen directamente de esas elecciones en el espectro.

5.6 Wavelets y Multirresolución

La Transformada de Fourier descompone la señal en frecuencias globales: cada coeficiente \(F(u,v)\) recibe contribuciones de toda la imagen, sin información explícita sobre la ubicación espacial de esas frecuencias. Así, las estructuras localizadas, como los bordes, se representan de forma distribuida en el espectro.

Las wavelets (ondículas) superan esta limitación al utilizar funciones base localizadas en el espacio, que pueden ser desplazadas y escaladas. Estas funciones poseen soporte compacto, es decir, son diferentes de cero solo en una región finita del dominio, permitiendo una representación simultánea en términos de frecuencia y ubicación espacial.

5.6.1 El Límite de la Transformada de Fourier: localización espacial

La Transformada de Fourier describe con precisión qué frecuencias están presentes en una señal, pero no representa explícitamente dónde ocurren esas frecuencias en el espacio.

En el experimento presentado en la Figura 5.16, dos imágenes con estructuras localizadas en posiciones diferentes producen espectros de magnitud prácticamente idénticos. Esto ocurre porque la representación de Fourier es global: cada coeficiente recibe contribución de toda la imagen.

Como consecuencia, el espectro de magnitud no representa explícitamente la localización de bordes u otras estructuras, solo la distribución de las frecuencias presentes. Esta limitación motivó el desarrollo de representaciones multirresolución, como la Transformada Wavelet Discreta (DWT), capaces de describir simultáneamente la frecuencia y la localización espacial de las estructuras de la imagen.

%%writefile tmp/fig_05_fracasso_fourier.cpp
#define MM_OUT "tmp/fig_05_fracasso_fourier.png"
//| label: fig-05-fracasso-fourier
//| fig-cap: "Fourier global es ciego a la posición. Los espectros no dicen dónde están los bordes."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>

int main() {
    // Crear img_sinal1: ceros de 128x128, columna 20-25 = 1, fila 100-105 = 1
    cv::Mat img_sinal1 = cv::Mat::zeros(128, 128, CV_8UC1);
    img_sinal1.colRange(20, 25).setTo(1);
    img_sinal1.rowRange(100, 105).setTo(1);

    // Crear img_sinal2: ceros de 128x128, columna 90-95 = 1, fila 30-35 = 1
    cv::Mat img_sinal2 = cv::Mat::zeros(128, 128, CV_8UC1);
    img_sinal2.colRange(90, 95).setTo(1);
    img_sinal2.rowRange(30, 35).setTo(1);

    // mag1 = log(1+|FFT shift|)
    cv::Mat img1_64f, img2_64f;
    img_sinal1.convertTo(img1_64f, CV_64F);
    img_sinal2.convertTo(img2_64f, CV_64F);

    cv::Mat fft1, fft2;
    cv::dft(img1_64f, fft1, cv::DFT_COMPLEX_OUTPUT);
    cv::dft(img2_64f, fft2, cv::DFT_COMPLEX_OUTPUT);

    // Función para fftshift de espectros
    auto fftshift = [](cv::Mat& spectrum) {
        cv::Mat shifted;
        int cx = spectrum.cols / 2;
        int cy = spectrum.rows / 2;

        // Cuadrantes para intercambiar
        cv::Mat q0(spectrum, cv::Rect(0, 0, cx, cy));
        cv::Mat q1(spectrum, cv::Rect(cx, 0, spectrum.cols - cx, cy));
        cv::Mat q2(spectrum, cv::Rect(0, cy, cx, spectrum.rows - cy));
        cv::Mat q3(spectrum, cv::Rect(cx, cy, spectrum.cols - cx, spectrum.rows - cy));

        cv::Mat tmp;
        q0.copyTo(tmp);
        q3.copyTo(q0);
        tmp.copyTo(q3);

        q1.copyTo(tmp);
        q2.copyTo(q1);
        tmp.copyTo(q2);
    };

    fftshift(fft1);
    fftshift(fft2);

    cv::Mat mag1, mag2;
    cv::Mat planes1[2], planes2[2];
    cv::split(fft1, planes1);
    cv::split(fft2, planes2);
    cv::magnitude(planes1[0], planes1[1], mag1);
    cv::magnitude(planes2[0], planes2[1], mag2);
    cv::log(1.0 + mag1, mag1);
    cv::log(1.0 + mag2, mag2);

    // Normalizar y convertir a 8-bit para visualización
    cv::Mat mag1_8u, mag2_8u;
    cv::normalize(mag1, mag1_8u, 0, 255, cv::NORM_MINMAX, CV_8U);
    cv::normalize(mag2, mag2_8u, 0, 255, cv::NORM_MINMAX, CV_8U);

    // Convertir imágenes para mostrar
    mm::Image img1_out(img_sinal1);
    mm::Image img2_out(img_sinal2);
    mm::Image mag1_out(mag1_8u);
    mm::Image mag2_out(mag2_8u);

    mm::show({img1_out, mag1_out, img2_out, mag2_out},
             MM_OUT,
             {"Sinal A", "Espectro A", "Sinal B (Desplazado)", "Espectro B"},
             4);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(img_sinal1, "tmp/fig_05_fracasso_fourier_0.png");
mm::write(mag1, "tmp/fig_05_fracasso_fourier_1.png");
mm::write(img_sinal2, "tmp/fig_05_fracasso_fourier_2.png");
mm::write(mag2, "tmp/fig_05_fracasso_fourier_3.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_fracasso_fourier.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_fracasso_fourier.cpp -o tmp/fig_05_fracasso_fourier -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_fracasso_fourier \
  && test -f "tmp/fig_05_fracasso_fourier.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_fracasso_fourier.png"

[1] Sinal A
[2] Espectro A
[3] Sinal B (Desplazado)
[4] Espectro B

try:
    mm.show(
        [
            mm.read("tmp/fig_05_fracasso_fourier_0.png"),
            mm.read("tmp/fig_05_fracasso_fourier_1.png"),
            mm.read("tmp/fig_05_fracasso_fourier_2.png"),
            mm.read("tmp/fig_05_fracasso_fourier_3.png"),
        ],
        titles=[
            'Sinal A',
            'Espectro A',
            'Sinal B (Deslocado)',
            'Espectro B',
        ],
        cols=4,
        figsize=(14, 4),
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_fracasso_fourier_0.png (ver a versao Python)")

Figura 5.16: Fourier global é cego para a posição. Os espectros não dizem onde as bordas estão.

5.6.2 Transformada Wavelet Discreta 2D

La Transformada Wavelet Discreta (DWT) aplica, de manera separada en las direcciones horizontal y vertical, dos filtros complementarios: un pasa-bajos \(h\) (aproximación) y un pasa-altos \(g\) (detalles), seguidos de un submuestreo por un factor de 2 en cada dimensión. Este proceso produce cuatro subbandas, cuyos nombres indican la combinación de los filtros aplicados en cada dirección (L = Low-pass, pasa-bajos; H = High-pass, pasa-altos). Las características de cada subbanda se resumen en la Tabla 5.3.

\[ \text{DWT}(f)=\{\underbrace{\text{LL}}_{\text{aprox.}},\; \underbrace{\text{LH}}_{\text{detalles horizontales}},\; \underbrace{\text{HL}}_{\text{detalles verticales}},\; \underbrace{\text{HH}}_{\text{detalles diagonales}}\}. \]

Tabla 5.3: Subbandas producidas por la Transformada Wavelet Discreta 2D (DWT), indicando los filtros aplicados en cada dirección y el contenido predominante de cada componente.

Subbanda	Filtros aplicados	Contenido visual
LL	baja × baja	Aproximación de la imagen (versión suavizada y reducida)
LH	baja × alta	Bordes horizontales y variaciones verticales
HL	alta × baja	Bordes verticales y variaciones horizontales
HH	alta × alta	Detalles diagonales y texturas

La descomposición puede aplicarse recursivamente sobre la subbanda LL, generando una representación multirresolución. Tras \(J\) niveles, se obtiene una estructura con \(3J+1\) subbandas, donde cada nuevo nivel reduce la resolución de la componente de aproximación.

NotaConexión con CNNs

La descomposición multirresolución de las wavelets posee una relación conceptual con las representaciones jerárquicas utilizadas en redes neuronales convolucionales (CNNs). En ambos casos, sucesivas etapas de filtrado y reducción de resolución producen descripciones cada vez más abstractas de la imagen. Sin embargo, las wavelets utilizan filtros matemáticamente definidos y reconstruibles, mientras que las CNNs aprenden sus filtros durante el entrenamiento.

5.6.3 Familias de Wavelets

Diferentes familias de wavelets presentan compromisos distintos entre soporte espacial, suavidad y capacidad de compresión. El soporte corresponde a la extensión de la función wavelet en el dominio espacial: cuanto menor es el soporte, más localizada está la función; cuanto mayor es, más suave tiende a ser su representación, aunque con un mayor costo computacional. La Tabla 5.4 compara algunas de las familias más utilizadas.

Tabla 5.4: Comparación entre familias de wavelets, destacando la longitud del soporte, el número de momentos nulos, la simetría y las aplicaciones típicas.

Wavelet	Longitud del soporte	Momentos nulos	Simetría	Uso típico
Haar	2	1	Asimétrica	Introducción y análisis básico
Daubechies db4	8	4	Asimétrica	Compresión y análisis general
Symlet sym4	8	4	Casi simétrica	Reconstrucción de señales
Biortogonal 5/3	5/3	2/2	Simétrica	JPEG 2000 sin pérdida
Biortogonal 9/7	9/7	4/4	Simétrica	JPEG 2000 con pérdida

Los momentos nulos miden la capacidad de la wavelet para representar regiones suaves de la imagen con pocos coeficientes distintos de cero. Una wavelet con \(p\) momentos nulos anula exactamente polinomios de grado hasta \(p-1\). En consecuencia, cuanto mayor es el número de momentos nulos, mayor tiende a ser la eficiencia de compresión en regiones homogéneas, aunque esto generalmente implica funciones con un soporte más largo.

La Figura 5.17 presenta las funciones de base (wavelets) \(\psi(t)\) en el dominio espacial. Estas funciones poseen soporte compacto, es decir, son distintas de cero solo en una región finita del dominio, a diferencia de las sinusoides de la Transformada de Fourier, que se extienden por todo el dominio.

%%writefile tmp/fig_05_wavelet_functions.cpp
#define MM_OUT "tmp/fig_05_wavelet_functions.png"
#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>

int main() {
    // psi via mm.wavefun (algoritmo em cascata) — suporte compacto: decai a zero.
    std::vector<double> x_h, phi_h, psi_h;
    mm::wavefun("haar", 6, x_h, phi_h, psi_h);
    std::vector<double> x_d, phi_d, psi_d;
    mm::wavefun("db4", 6, x_d, phi_d, psi_d);

    mm::Image chart_h = mm::lineChart(
        std::vector<std::vector<double>>{x_h},
        std::vector<std::vector<double>>{psi_h},
        std::vector<std::string>{"psi Haar"},
        std::vector<cv::Scalar>{cv::Scalar(180, 60, 40)},
        "Ondaleta Haar (psi)", "t");

    mm::Image chart_d = mm::lineChart(
        std::vector<std::vector<double>>{x_d},
        std::vector<std::vector<double>>{psi_d},
        std::vector<std::string>{"psi Daubechies 4"},
        std::vector<cv::Scalar>{cv::Scalar(60, 140, 40)},
        "Ondaleta Daubechies 4 (psi)", "t");

    mm::show(std::vector<mm::Image>{chart_h, chart_d}, MM_OUT,
             std::vector<std::string>{"Haar", "Daubechies 4"}, 2);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(chart_h, "tmp/fig_05_wavelet_functions_0.png");
mm::write(chart_d, "tmp/fig_05_wavelet_functions_1.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_wavelet_functions.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_wavelet_functions.cpp -o tmp/fig_05_wavelet_functions -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_wavelet_functions \
  && test -f "tmp/fig_05_wavelet_functions.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_wavelet_functions.png"

[1] Haar
[2] Daubechies 4

try:
    mm.show(
        [
            mm.read("tmp/fig_05_wavelet_functions_0.png"),
            mm.read("tmp/fig_05_wavelet_functions_1.png"),
        ],
        titles=[
            'Haar',
            'Daubechies 4',
        ],
        cols=2,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_wavelet_functions_0.png (ver a versao Python)")

Figura 5.17: Funções da Wavelet (ψ). Note como elas rapidamente decaem para zero (suporte compacto), ao contrário das senoides infinitas de Fourier.

El diagrama de la Figura 5.18 ilustra el análisis multirresolución realizado por la DWT, en el cual la subbanda de aproximación (LL) se descompone sucesivamente, formando una representación jerárquica con dos niveles.

Descomposición Wavelet 2D — Estructura Multirresolución (2 niveles)

Figura 5.18: Diagrama de la descomposición wavelet 2D en dos niveles.

El simulador de la Figura 5.19 permite explorar interactivamente la Transformada Wavelet Discreta 2D (DWT) utilizando la wavelet de Haar. La descomposición en subbandas evidencia la separación entre la componente de aproximación y las componentes de detalle de la imagen.

Los diferentes patrones de entrada permiten observar el comportamiento direccional de los filtros. En imágenes con bordes horizontales y verticales, las subbandas LH y HL destacan, respectivamente, las variaciones verticales y horizontales de la intensidad. En regiones de variación suave, la mayor parte de la energía se concentra en la subbanda de aproximación LL, mientras que las subbandas de detalle presentan coeficientes cercanos a cero.

El análisis multirresolución también se puede observar al aumentar el número de niveles de descomposición. En este caso, solo la subbanda \(\text{LL}_1\) se descompone nuevamente, dando origen a las subbandas \(\text{LL}_2\), \(\text{LH}_2\), \(\text{HL}_2\) y \(\text{HH}_2\), que forman el segundo nivel de la representación jerárquica.

En patrones formados por regiones homogéneas de gran extensión, como un degradado suave o un tablero compuesto por bloques grandes, la energía permanece predominantemente concentrada en la subbanda LL. En el degradado, esto ocurre porque las diferencias entre píxeles vecinos son pequeñas. En el tablero, por su parte, los píxeles poseen prácticamente la misma intensidad en el interior de cada bloque, de modo que solo las fronteras entre bloques producen coeficientes no nulos en las subbandas de detalle. Como estas fronteras ocupan solo una pequeña fracción de la imagen, su contribución a la energía total permanece reducida.

Para posibilitar el análisis visual de estas variaciones sutiles, el simulador incorpora un control de ganancia de contraste de los detalles (que varía de 1 a 8). Este parámetro funciona como un factor de amplificación lineal aplicado exclusivamente a los coeficientes de las subbandas de detalle (LH, HL y HH) antes de su renderización en pantalla. En escenarios de transición suave (como el gradiente) o de uniformidad local (como el interior de los bloques del tablero), las diferencias numéricas calculadas por el filtro paso-alto de Haar resultan en coeficientes muy cercanos a cero, lo que haría que los cuadrantes correspondientes fueran oscuros e imperceptibles a simple vista. Al multiplicar estos valores por la ganancia, el simulador recupera visualmente las estructuras de alta frecuencia ocultas y resalta la orientación de los bordes remanentes.

El gráfico de energía por subbanda cuantifica esta distribución entre la componente de aproximación y las componentes de detalle, demostrando que la ganancia visual no altera la métrica original de la energía. En imágenes naturales, la mayor parte de la energía se concentra en la subbanda LL, mientras que las subbandas LH, HL y HH representan principalmente bordes, texturas y otras variaciones locales de la intensidad.

🌊 Simulador: Descomposición Wavelet 2D Transformada de Haar

Patrón Sintético Combinado (formas + textura) Formas (bordes h/v) Textura (alta frecuencia) Gradiente suave

Niveles de Descomposición

1 nivel 2 niveles

Ganancia de Contraste: 3.0

Imagen Original

Descomposición Wavelet (Mosaico)

Energía por Subbanda (%) — Suma Preservada (Parseval)

LL — Aproximación

Versión suavizada y reducida de la imagen

LH — Detalle Horizontal

Realza bordes horizontales (variación vertical)

HL — Detalle Vertical

Realza bordes verticales (variación horizontal)

HH — Detalle Diagonal

Texturas y esquinas (variación en ambas direcciones)

Figura 5.19: Simulación de la descomposición wavelet 2D.

5.6.4 Análisis Multirresolución con la DWT 2D

La Transformada Wavelet Discreta 2D (DWT) descompone una imagen en componentes de aproximación y detalle, organizadas de forma jerárquica en diferentes escalas y orientaciones. Como las subbandas de detalle en imágenes naturales frecuentemente presentan coeficientes de bajo contraste, los ejemplos prácticos a continuación utilizan un patrón geométrico sintético generado en Python. Este enfoque replica el comportamiento del simulador de la Figura 5.19, haciendo visualmente explícitos los efectos del filtrado espacial y de la descomposición multirresolución.

5.6.4.1 Descomposición en Mosaico de Múltiples Niveles

La Figura 5.20 ilustra la estructura jerárquica de la DWT en dos niveles utilizando la wavelet de Haar. El proceso se basa en la aplicación combinada de filtros pasa-baja y pasa-alta en las direcciones horizontal y vertical, seguidos por submuestreo por un factor de 2.

En el primer nivel, la imagen original origina la subbanda de aproximación (\(LL_1\)) y las componentes de detalle horizontal (\(LH_1\)), vertical (\(HL_1\)) y diagonal (\(HH_1\)). En el análisis multirresolución, la subbanda \(LL_1\) se filtra y submuestrea nuevamente, generando el segundo nivel de descomposición (\(LL_2\), \(LH_2\), \(HL_2\) y \(HH_2\)).

Para viabilizar la interpretación visual de las componentes de detalle, el código extrae el valor absoluto de sus coeficientes y aplica una normalización lineal (min-max) para ocupar toda la gama dinámica de tonos de gris [0, 255]. Esta operación transforma regiones homogéneas (coeficientes nulos) en negro y destaca en blanco los bordes y texturas extraídas en cada escala y orientación.

%%writefile tmp/fig_05_dwt_subbandas.cpp
#define MM_OUT "tmp/fig_05_dwt_subbandas.png"
// Compile: g++ -std=c++17 -o snippet snippet.cpp $(pkg-config --cflags --libs opencv4)
#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include "morph.hpp"

//| label: fig-05-dwt-subbandas
//| fig-cap: "Decomposição *wavelet* 2D de 2 níveis com *wavelet* Haar: subbandas LL, LH, HL, HH em cada nível. As subbandas de detalhe revelam estruturas orientadas em diferentes escalas utilizando um padrão sintético."
//| echo: true
//| output: true

// ── Geração da Imagem Sintética (Mesmo padrão 'combined' do simulador) ────────
cv::Mat gerar_imagem_sintetica(int N = 256) {
    cv::Mat img = cv::Mat::zeros(N, N, CV_64F);
    for (int y = 0; y < N; ++y) {
        for (int x = 0; x < N; ++x) {
            double v = 55 + 35 * (static_cast<double>(x) / N) + 15 * std::sin(static_cast<double>(y) / 24);
            // Quadrado
            if (24 < x && x < 100 && 24 < y && y < 100) {
                v = 225;
            }
            // Círculo
            int cx = 190, cy = 76, r = 34;
            if ((x - cx) * (x - cx) + (y - cy) * (y - cy) < r * r) {
                v = 205;
            }
            // Textura periódica (inferior)
            if (y > 164 && y < 244) {
                int p = 12;
                v = (((x / p) + (y / p)) % 2 == 0) ? 185.0 : 65.0;
            }
            // Linha diagonal
            if (std::abs(x - y) < 4) {
                v = 240;
            }
            v = std::max(0.0, std::min(255.0, v));
            img.at<double>(y, x) = v;
        }
    }
    cv::Mat img8;
    img.convertTo(img8, CV_8U);
    return img8;
}

// Normaliza subbanda para visualização [0,255]
cv::Mat sb_vis(const cv::Mat& sb) {
    cv::Mat abs_sb;
    cv::absdiff(sb, cv::Scalar(0), abs_sb);
    cv::Mat normalized;
    cv::normalize(abs_sb, normalized, 0, 255, cv::NORM_MINMAX, CV_8U);
    return normalized;
}

int main() {
    // Substitui a imagem escura de moedas pelo padrão sintético claro
    cv::Mat img_gray = gerar_imagem_sintetica(256);

    // ── Decomposição wavelet 2 níveis ─────────────────────────────────────────────
    std::string wavelet = "haar";
    cv::Mat img_float;
    img_gray.convertTo(img_float, CV_64F);

    // Nível 1
    mm::Subbands coefs1 = mm::dwt2(img_float, wavelet);
    cv::Mat LL1 = coefs1.LL;
    cv::Mat LH1 = coefs1.LH;
    cv::Mat HL1 = coefs1.HL;
    cv::Mat HH1 = coefs1.HH;

    // Nível 2 (aplicado sobre LL1)
    mm::Subbands coefs2 = mm::dwt2(LL1, wavelet);
    cv::Mat LL2 = coefs2.LL;
    cv::Mat LH2 = coefs2.LH;
    cv::Mat HL2 = coefs2.HL;
    cv::Mat HH2 = coefs2.HH;

    std::cout << "Forma original     : " << img_gray.rows << "x" << img_gray.cols << std::endl;
    std::cout << "LL1 (nível 1)      : " << LL1.rows << "x" << LL1.cols << "  |  LH1/HL1/HH1: " << LH1.rows << "x" << LH1.cols << std::endl;
    std::cout << "LL2 (nível 2)      : " << LL2.rows << "x" << LL2.cols << "    |  LH2/HL2/HH2: " << LH2.rows << "x" << LH2.cols << std::endl;

    std::vector<mm::Image> imgs_dwt = {mm::Image(img_gray), 
                                        mm::Image(sb_vis(LL1)), mm::Image(sb_vis(LH1)), 
                                        mm::Image(sb_vis(HL1)), mm::Image(sb_vis(HH1)),
                                        mm::Image(sb_vis(LL2)), mm::Image(sb_vis(LH2)), 
                                        mm::Image(sb_vis(HL2)), mm::Image(sb_vis(HH2))};
    std::vector<std::string> titles_dwt = {"Original",
                                            "LL₁ (aprox.)", "LH₁ (horiz.)", "HL₁ (vert.)", "HH₁ (diag.)",
                                            "LL₂ (aprox.)", "LH₂ (horiz.)", "HL₂ (vert.)", "HH₂ (diag.)"};

    mm::show(imgs_dwt, MM_OUT, titles_dwt, 5);

    return 0;
}

Overwriting tmp/fig_05_dwt_subbandas.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dwt_subbandas.cpp -o tmp/fig_05_dwt_subbandas -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dwt_subbandas \
  && test -f "tmp/fig_05_dwt_subbandas.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dwt_subbandas.png"

Forma original     : 256x256
LL1 (nível 1)      : 128x128  |  LH1/HL1/HH1: 128x128
LL2 (nível 2)      : 64x64    |  LH2/HL2/HH2: 64x64
[1] Original
[2] LL₁ (aprox.)
[3] LH₁ (horiz.)
[4] HL₁ (vert.)
[5] HH₁ (diag.)
[6] LL₂ (aprox.)
[7] LH₂ (horiz.)
[8] HL₂ (vert.)
[9] HH₂ (diag.)

try:
    mm.show(mm.read("tmp/fig_05_dwt_subbandas.png"), figsize=(16, 7))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dwt_subbandas.png (ver a versao Python)")

Figura 5.20: Decomposição wavelet 2D de 2 níveis com wavelet Haar: subbandas LL, LH, HL, HH em cada nível. As subbandas de detalhe revelam estruturas orientadas em diferentes escalas utilizando um padrão sintético.

5.6.4.2 El Compromiso entre Localización y Suavidad

La elección de la función de base (wavelet) influye directamente en la forma en que las características de la imagen se distribuyen y codifican mediante los coeficientes de la DWT. La Figura 5.21 compara los resultados prácticos obtenidos al aplicar cuatro familias distintas sobre el patrón geométrico sintético: haar, db4, sym4 y bior2.2.

Por poseer soporte corto y formato de función escalón, la wavelet de Haar produce coeficientes altamente localizados en las discontinuidades espaciales, generando bordes finos y nítidos en las subbandas de detalle. En contrapartida, familias como Daubechies (db4) y Symlets (sym4), que presentan mayor soporte (filtros más largos) y mayor número de momentos nulos, generan respuestas más suaves y distribuidas alrededor de las transiciones, lo que puede introducir ligeras oscilaciones o desenfoques en las fronteras abruptas.

Este comportamiento evidencia el clásico compromiso (trade-off) del análisis de multirresolución: soportes menores favorecen la localización espacial exacta de los bordes, mientras que soportes mayores y un mayor número de momentos nulos tienden a producir representaciones más dispersas y suaves. Esa suavidad y capacidad de atenuación de altas frecuencias garantizan una mayor eficiencia en la compactación de la energía, características fundamentales para aplicaciones de compresión de datos y eliminación de ruido (denoising).

%%writefile tmp/fig_05_dwt_wavelets.cpp
#define MM_OUT "tmp/fig_05_dwt_wavelets.png"
//| label: fig-05-dwt-wavelets
//| fig-cap: "Comparação entre famílias de *wavelets*: Haar, db4, sym4 e bior2.2. Subbanda LL₁ (aproximação) e HH₁ (diagonal) para cada escolha, ilustrando o compromisso entre compactação e suavidade com base no padrão sintético."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include "morph.hpp"

// Garante que img_gray e img_float utilizem o mesmo padrão sintético claro
// (Fallback caso o bloco anterior não tenha sido executado na mesma sessão)
static cv::Mat gerar_imagem_sintetica(int N = 256) {
    cv::Mat img(N, N, CV_64F);
    for (int y = 0; y < N; ++y) {
        for (int x = 0; x < N; ++x) {
            double v = 55 + 35 * ((double)x / N) + 15 * std::sin(y / 24.0);
            if (24 < x && x < 100 && 24 < y && y < 100) v = 225;
            int cx = 190, cy = 76, r = 34;
            if ((x - cx)*(x - cx) + (y - cy)*(y - cy) < r*r) v = 205;
            if (y > 164 && y < 244) {
                int p = 12;
                v = (((x / p + y / p) % 2) == 0) ? 185 : 65;
            }
            if (std::abs(x - y) < 4) v = 240;
            img.at<double>(y, x) = std::max(0.0, std::min(255.0, v));
        }
    }
    cv::Mat img8;
    img.convertTo(img8, CV_8U);
    return img8;
}

// Visor de subbanda (assume que existe em morph.hpp ou é implementado aqui)
static mm::Image sb_vis(const cv::Mat& sb) {
    cv::Mat normalized;
    cv::normalize(sb, normalized, 0, 255, cv::NORM_MINMAX, CV_8U);
    return mm::Image(normalized);
}

int main() {
    cv::Mat img_gray;
    // Se a variável global estiver disponível (injetada por célula anterior), use-a;
    // caso contrário, gere o fallback.
    // (Na prática, como não há persistência real entre células neste contexto, 
    //  sempre geramos o fallback.)
    img_gray = gerar_imagem_sintetica(256);

    cv::Mat img_float;
    img_gray.convertTo(img_float, CV_64F);

    std::vector<std::string> wavelets_comp = {"haar", "db4", "sym4", "bior2.2"};
    std::vector<mm::Image> imgs_comp;
    std::vector<std::string> titles_comp;

    for (const auto& wname : wavelets_comp) {
        // pywt.dwt2 -> mm::Subbands via morph.hpp
        mm::Subbands s = mm::dwt2(img_float, wname);
        cv::Mat LL = s.LL;
        cv::Mat HH = s.HH;

        imgs_comp.push_back(sb_vis(LL));
        imgs_comp.push_back(sb_vis(HH));
        titles_comp.push_back(wname + " — LL₁");
        titles_comp.push_back(wname + " — HH₁");
    }

    mm::show(imgs_comp, MM_OUT, titles_comp, 4);
    return 0;
}

Overwriting tmp/fig_05_dwt_wavelets.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dwt_wavelets.cpp -o tmp/fig_05_dwt_wavelets -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dwt_wavelets \
  && test -f "tmp/fig_05_dwt_wavelets.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dwt_wavelets.png"

[1] haar — LL₁
[2] haar — HH₁
[3] db4 — LL₁
[4] db4 — HH₁
[5] sym4 — LL₁
[6] sym4 — HH₁
[7] bior2.2 — LL₁
[8] bior2.2 — HH₁

try:
    mm.show(mm.read("tmp/fig_05_dwt_wavelets.png"), figsize=(14, 8))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dwt_wavelets.png (ver a versao Python)")

Figura 5.21: Comparação entre famílias de wavelets: Haar, db4, sym4 e bior2.2. Subbanda LL₁ (aproximação) e HH₁ (diagonal) para cada escolha, ilustrando o compromisso entre compactação e suavidade com base no padrão sintético.

5.6.4.3 Limiarización de Coeficientes y Compresión

Una de las principales aplicaciones de la Transformada Wavelet Discreta (DWT) es la compresión de datos, impulsada por la capacidad de representación dispersa de los coeficientes. La Figura 5.22 ilustra el efecto de la limiarización abrupta (hard thresholding), técnica en la cual los coeficientes de detalle con magnitud inferior a un umbral \(T\) se anulan por completo antes del proceso de síntesis realizado por la Transformada Wavelet Discreta Inversa (IDWT).

A medida que se eleva el umbral \(T\), un volumen creciente de coeficientes de alta frecuencia se pone a cero. Al concentrar menor energía, la eliminación de estos componentes reduce considerablemente la cantidad de información necesaria para representar la imagen, manteniendo intacta la componente de aproximación global (la subbanda \(LL\) más profunda) para preservar la estructura macro. Visualmente, este descarte de coeficientes se manifiesta a través de la desaparición progresiva de texturas finas y el suavizado de transiciones abruptas de intensidad.

La fidelidad de la imagen reconstruida frente a la original se cuantifica mediante la métrica de Pico de Relación Señal-Ruido (PSNR, Peak Signal-to-Noise Ratio), expresada en decibelios (dB). Valores más altos de PSNR indican menor distorsión y mayor proximidad matemática con la señal original. El experimento práctico evidencia la disminución gradual del PSNR conforme aumenta la agresividad de la limiarización, lo que permite evaluar numéricamente el umbral óptimo para el equilibrio entre compresión y degradación visual.

%%writefile tmp/fig_05_dwt_reconstrucao.cpp
#define MM_OUT "tmp/fig_05_dwt_reconstrucao.png"
// Compile: g++ -std=c++17 -O2 snippet.cpp -o snippet $(pkg-config --cflags --libs opencv4)
#include <opencv2/opencv.hpp>
#include <string>
#include <vector>
#include <cmath>
#include <iostream>
#include "morph.hpp"

//| label: fig-05-dwt-reconstrucao
//| fig-cap: "Reconstrução *wavelet* com limiarização de coeficientes (*hard thresholding*): à medida que o limiar aumenta, mais detalhes são zerados, produzindo imagens progressivamente mais suaves. Métrica PSNR quantifica a perda de qualidade sobre o padrão sintético."
//| echo: true
//| output: true

// Garante que img_gray utilize o mesmo padrão sintético claro
// (Função auxiliar para gerar o padrão sintético)
cv::Mat gerar_imagem_sintetica(int N = 256) {
    cv::Mat img = cv::Mat::zeros(N, N, CV_32F);
    for (int y = 0; y < N; ++y) {
        for (int x = 0; x < N; ++x) {
            double v = 55 + 35 * (double(x) / N) + 15 * std::sin(y / 24.0);
            if (24 < x && x < 100 && 24 < y && y < 100) v = 225;
            double cx = 190, cy = 76, r = 34;
            if ((x - cx) * (x - cx) + (y - cy) * (y - cy) < r * r) v = 205;
            if (y > 164 && y < 244) {
                int p = 12;
                v = (((x / p + y / p) % 2) == 0) ? 185 : 65;
            }
            if (std::abs(x - y) < 4) v = 240;
            img.at<float>(y, x) = static_cast<float>(std::max(0.0, std::min(v, 255.0)));
        }
    }
    cv::Mat out;
    img.convertTo(out, CV_8U);
    return out;
}

cv::Mat dwt_threshold_reconstruct(const cv::Mat& img, const std::string& wavelet = "db4", int nivel = 2, double threshold = 0.0) {
    // Decompõe, aplica limiar e reconstrói via IDWT.
    cv::Mat img64;
    img.convertTo(img64, CV_64F);
    mm::WaveDec2 c = mm::wavedec2(img64, wavelet, nivel);

    // Copia e aplica hard thresholding em todos os detalhes
    // LL final não é limiarizado
    std::vector<cv::Mat> detail_t;
    // Reconstruir a estrutura de coeficientes
    // c_rec será construído no formato esperado por waverec2
    mm::WaveDec2 c_rec = c; // copia inicial

    // Aplicar limiar aos detalhes (nível mais grosseiro ao mais fino)
    for (int j = 0; j < nivel; ++j) {
        int idx = nivel - 1 - j; // ordem: from coarse to fine
        auto& detalhe = c.detail[idx];
        for (int k = 0; k < 3; ++k) {
            c_rec.detail[idx][k] = mm::wave_threshold(detalhe[k], threshold, "hard");
        }
    }

    cv::Mat rec = mm::waverec2(c_rec, wavelet);

    // Recorte para dimensão original
    cv::Mat rec_crop = rec(cv::Rect(0, 0, img.cols, img.rows)).clone();

    // Clip e converter para uint8
    cv::Mat rec8;
    cv::normalize(rec_crop, rec_crop, 0, 255, cv::NORM_MINMAX);
    rec_crop.convertTo(rec8, CV_8U);
    return rec8;
}

int main() {
    cv::Mat img_gray = gerar_imagem_sintetica(256);

    std::vector<int> thresholds = {0, 10, 30, 60, 100};
    std::vector<cv::Mat> imgs_thr = {img_gray};
    std::vector<std::string> titles_thr = {"Original"};

    for (int t : thresholds) {
        cv::Mat rec = dwt_threshold_reconstruct(img_gray, "db4", 2, static_cast<double>(t));
        double psnr = cv::PSNR(img_gray, rec);
        imgs_thr.push_back(rec);
        titles_thr.push_back("T=" + std::to_string(t) + "  PSNR=" + std::to_string(psnr).substr(0, 5) + " dB");
    }

    // Convert vector<cv::Mat> to vector<mm::Image> for display
    std::vector<mm::Image> imgs_display;
    for (const auto& m : imgs_thr) {
        imgs_display.push_back(mm::Image(m));
    }
    mm::show(imgs_display, MM_OUT, titles_thr, 3);

    return 0;
}

Overwriting tmp/fig_05_dwt_reconstrucao.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dwt_reconstrucao.cpp -o tmp/fig_05_dwt_reconstrucao -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dwt_reconstrucao \
  && test -f "tmp/fig_05_dwt_reconstrucao.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dwt_reconstrucao.png"

[1] Original
[2] T=0  PSNR=19.54 dB
[3] T=10  PSNR=19.57 dB
[4] T=30  PSNR=22.94 dB
[5] T=60  PSNR=21.42 dB
[6] T=100  PSNR=18.46 dB

try:
    mm.show(mm.read("tmp/fig_05_dwt_reconstrucao.png"), figsize=(14, 10))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dwt_reconstrucao.png (ver a versao Python)")

Figura 5.22: Reconstrução wavelet com limiarização de coeficientes (hard thresholding): à medida que o limiar aumenta, mais detalhes são zerados, produzindo imagens progressivamente mais suaves. Métrica PSNR quantifica a perda de qualidade sobre o padrão sintético.

Síntesis — Fourier vs. Wavelets: ¿cuándo utilizar cada enfoque?

La Tabla 5.5 sintetiza las principales diferencias estructurales y operativas entre la Transformada Discreta de Fourier (DFT) y la Transformada Wavelet Discreta (DWT).

Tabla 5.5: Comparación entre la Transformada Discreta de Fourier (DFT) y la Transformada Wavelet Discreta (DWT), destacando sus principales características y aplicaciones.

Criterio	Fourier (DFT)	Wavelet (DWT)
Funciones de base	Senoides de soporte infinito	Funciones de soporte compacto
Localización espacial	No explícita (global)	Explícita (local)
Filtrado espectral	Excelente para el control fino de frecuencias	Basado en subbandas (escalas)
Compresión de imágenes	Base de la DCT (JPEG tradicional)	Base de la DWT (JPEG 2000)
Análisis multiescala	No	Sí
Eliminación de ruido periódico	Altamente eficiente	Poco indicada
Señales no estacionarias	Limitada	Altamente eficiente

En términos prácticos, la DFT se consolida como la herramienta ideal para el análisis espectral puro, el diseño de filtros selectivos en el dominio de la frecuencia y la atenuación de ruidos periódicos y armónicos. Por otro lado, la DWT sobresale en escenarios que exigen la preservación rigurosa de la localización espacial de las características asociada a su contenido frecuencial, destacándose en la compresión de datos, el análisis multirresolución y el procesamiento de transiciones abruptas. De este modo, ambas transformadas deben entenderse como técnicas perfectamente complementarias, que mapean caminos distintos y específicos para la resolución de problemas en PDI-VC.

NotaAnalogías con Audio: Limitaciones y Precauciones

Al establecer analogías entre el procesamiento de imágenes y el audio, es importante considerar las diferencias fundamentales:

• En los sistemas de audio estéreo/multicanal, la fase entre canales es crucial para la percepción de la localización espacial (diferencias interaurales de fase y tiempo).

• En los sistemas monoaurales, la fase tiene una influencia perceptual limitada: el oído humano es relativamente insensible a la fase absoluta de componentes sinusoidales aislados.

• En las imágenes, la fase de la DFT siempre es fundamental para la localización espacial de las estructuras, independientemente de que se trate de una imagen monocromática o en color.

La analogía entre la fase en audio y la fase en imágenes debe emplearse con cautela, destacando que, aunque ambas transportan información sobre la organización espacial/temporal de la señal, los mecanismos perceptuales son fundamentalmente diferentes.

5.7 Compresión de Imágenes

Mientras que las wavelets establecen la base teórica del estándar JPEG 2000, el estándar JPEG tradicional se basa en la Transformada Discreta de Cosenos (DCT, Discrete Cosine Transform). A pesar de las diferencias estructurales, ambos enfoques comparten el mismo principio fundamental: compactar la energía de la imagen en un número reducido de coeficientes y descartar los componentes de menor relevancia con un impacto visual mínimo.

El objetivo central de la compresión es reducir el volumen de datos necesario para el almacenamiento o transmisión de una imagen. Este proceso es posible gracias a la identificación y eliminación de redundancias estructurales y perceptuales.

5.7.1 Taxonomía de las Redundancias

El desarrollo de algoritmos de compresión se fundamenta en la identificación y eliminación de tres categorías principales de redundancia, sintetizadas en la Tabla 5.6.

Tabla 5.6: Categorías de redundancia en imágenes digitales y sus respectivos mecanismos de explotación.

Tipo	Definición	Enfoque de Explotación
Espacial (interpíxel)	Alta correlación y dependencia estadística entre píxeles vecinos.	DCT, DWT y codificación predictiva.
Espectral (intercanal)	Correlación estadística entre los canales de color de una misma imagen.	Transformaciones de espacio de color (ej: RGB a \(YC_bC_r\)).
Psicovisual	Insensibilidad del sistema visual humano (SVH) a variaciones de alta frecuencia y bajo contraste.	Procesos de cuantización selectiva de coeficientes.

Dependiendo de la preservación de la información original tras el proceso de decodificación, los métodos de compresión se dividen en dos clases fundamentales:

• Sin pérdida (lossless): Garantiza una reconstrucción bit a bit idéntica a la imagen original. Se emplea en escenarios donde la integridad de los datos es estrictamente crítica, como en imágenes médicas, diagnósticos por imagen y almacenamiento de documentos textuales.
• Con pérdida (lossy): Admite la introducción de una distorsión controlada en la señal a cambio de tasas de compresión sustancialmente más elevadas. Es el enfoque estándar para fotografías de consumo y streaming de vídeo, ecosistemas en los cuales el SVH tolera pequeñas atenuaciones de alta frecuencia sin percepción de degradación de la calidad visual.

5.7.2 Transformada de Cosenos Discreta (DCT-II 2D)

La Transformada de Cosenos Discreta (DCT) constituye la operación central del estándar JPEG. A diferencia de la DFT, que utiliza una base compleja, la DCT se basa en funciones trigonométricas puramente reales. Para un bloque de imagen \(f(x,y)\) de dimensiones \(N \times N\), la DCT-II 2D mapea la señal espacial al dominio de las frecuencias espaciales, generando la matriz de coeficientes \(C(u,v)\) mediante:

\[ C(u,v) = \alpha(u)\,\alpha(v) \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y)\, \cos\!\left[\frac{\pi(2x+1)u}{2N}\right] \cos\!\left[\frac{\pi(2y+1)v}{2N}\right] \tag{5.9}\]

donde los factores de normalización ortogonal están dados por \(\alpha(0) = \sqrt{1/N}\) y \(\alpha(k) = \sqrt{2/N}\) para \(k > 0\).

Cada coeficiente \(C(u,v)\) cuantifica la contribución —o “peso”— de una frecuencia espacial específica dentro de ese bloque. El término \(C(0,0)\) se denomina componente DC y representa la intensidad media del bloque (frecuencia nula). Los demás coeficientes, llamados componentes AC (Alternating Current), corresponden a las frecuencias espaciales progresivamente mayores.

5.7.3 Las Funciones de Base de la DCT

Desde una perspectiva geométrica, la Ecuación 5.9 realiza la proyección del bloque de píxeles sobre un conjunto de funciones ortogonales. Para el caso estándar de JPEG (\(N=8\)), el bloque espacial se descompone en una combinación lineal de 64 funciones de base bidimensionales, denotadas por \(B_{u,v}(x,y)\) y generadas por el producto de funciones cosenoidales:

\[B_{u,v}(x,y) = \cos\left[ \frac{\pi (2x+1)u}{16} \right] \cos\left[ \frac{\pi (2y+1)v}{16} \right]\]

De esta manera, la operación inversa puede interpretarse como la reconstrucción exacta del bloque original mediante la suma ponderada de estas 64 matrices de base, donde cada coeficiente \(C(u,v)\) actúa como el peso analítico de su respectiva componente armónica.

La frecuencia espacial indicada por los índices \((u,v)\) determina el número de ciclos de oscilación a lo largo de las dimensiones horizontales y verticales del bloque. Como se ilustra en la Figura 5.23 —cuyo código aísla cada base aplicando la transformación inversa sobre impulsos unitarios—, estas 64 funciones se organizan en una matriz \(8 \times 8\). La esquina superior izquierda (\(u=0, v=0\)) muestra el patrón uniforme de frecuencia nula (DC), mientras que el avance hacia la derecha (eje \(u\)) o hacia abajo (eje \(v\)) mapea variaciones armónicas progresivamente mayores, representando transiciones rápidas, bordes y texturas en las orientaciones horizontales, verticales y diagonales.

NotaDCT vs DFT: Ventaja de la Compactación de Energía

Tanto la DCT como la DFT mapean un bloque espacial \(N \times N\) en una matriz de coeficientes de la misma dimensión. Sin embargo, para imágenes naturales, la DCT presenta una mayor eficiencia en la compactación de energía en las bajas frecuencias. Esto ocurre porque la DCT asume implícitamente una simetría par de la señal en las fronteras del bloque, lo que equivale a una extensión periódica continua, minimizando el efecto de dispersión espectral (ringing). Como consecuencia, la mayoría de los coeficientes AC decae rápidamente hacia valores cercanos a cero, optimizando el pipeline de compresión sin introducir degradación visual perceptible.

%%writefile tmp/fig_05_dct_basis.cpp
#define MM_OUT "tmp/fig_05_dct_basis.png"
//| label: fig-05-dct-basis
//| fig-cap: "O Alfabeto Visual do JPEG: As 64 funções de base da DCT-II. O coeficiente DC fica no topo esquerdo (suave). Ao descer e avançar à direita, a oscilação espacial aumenta drasticamente."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>

int main() {
    // Cada base é a IDCT de um único coeficiente unitário — montadas num mosaico 8×8.
    int tile = 32;
    cv::Mat montagem = cv::Mat::zeros(8 * tile, 8 * tile, CV_8UC1);

    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            cv::Mat coef = cv::Mat::zeros(8, 8, CV_64F);
            coef.at<double>(i, j) = 1.0;
            cv::Mat base = mm::idct2(coef);
            cv::Mat base_norm;
            cv::normalize(base, base_norm, 0, 255, cv::NORM_MINMAX);
            base_norm.convertTo(base_norm, CV_8U);
            cv::Mat base_resized;
            cv::resize(base_norm, base_resized, cv::Size(tile, tile), 0, 0, cv::INTER_NEAREST);

            cv::Rect roi(j * tile, i * tile, tile, tile);
            base_resized.copyTo(montagem(roi));
        }
    }

    mm::Image montagem_img(montagem);
    mm::show(std::vector<mm::Image>{montagem_img}, MM_OUT,
             std::vector<std::string>{"As 64 bases da DCT-II 8×8 (DC no topo-esquerdo)"}, 1);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(montagem, "tmp/fig_05_dct_basis_0.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_dct_basis.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dct_basis.cpp -o tmp/fig_05_dct_basis -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dct_basis \
  && test -f "tmp/fig_05_dct_basis.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dct_basis.png"

[1] As 64 bases da DCT-II 8×8 (DC no topo-esquerdo)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_dct_basis_0.png"),
        ],
        titles=[
            'As 64 bases da DCT-II 8×8 (DC no topo-esquerdo)',
        ],
        cols=1,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dct_basis_0.png (ver a versao Python)")

Figura 5.23: O Alfabeto Visual do JPEG: As 64 funções de base da DCT-II. O coeficiente DC fica no topo esquerdo (suave). Ao descer e avançar à direita, a oscilação espacial aumenta drasticamente.

5.7.4 Concentración de Energía y Reconstrucción Progresiva

Antes de la aplicación de la DCT, los píxeles del bloque de intensidad se trasladan rutinariamente (restando \(128\) para imágenes de 8 bits) con el fin de centrar la señal en torno a cero, eliminando componentes continuas innecesarias. Al calcular la DCT sobre el bloque resultante, la propiedad de compactación de energía se hace evidente: casi la totalidad de la varianza y de la información de la imagen original se concentra en el coeficiente DC (\(C(0,0)\)) y en los primeros armónicos AC de baja frecuencia.

La Figura 5.24 demuestra este fenómeno mediante una reconstrucción progresiva por truncamiento abrupto. En lugar de utilizar los 64 coeficientes, el algoritmo preserva únicamente los \(k\) primeros componentes —seleccionados con base en un barrido que prioriza las bajas frecuencias espaciales— y anula los restantes.

La síntesis inversa (IDCT) realizada con solo una fracción de los coeficientes (como 15% o 30%) ya es capaz de recuperar las estructuras y la iluminación macro del bloque original de píxeles. A medida que los armónicos de frecuencias más altas se reincorporan progresivamente, los detalles finos y las transiciones rápidas se restauran. Este comportamiento valida el principio de la compresión perceptual: las altas frecuencias descartadas poseen poca energía y su ausencia, en condiciones normales, genera un impacto visual secundario en la percepción del observador.

%%writefile tmp/fig_05_dct_bloco.cpp
#define MM_OUT "tmp/fig_05_dct_bloco.png"
#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include <iostream>
#include <algorithm>
#include "morph.hpp"

//| label: fig-05-dct-bloco
//| fig-cap: "DCT 2D em bloco 8×8: coeficientes e reconstrução progressiva."
//| echo: true
//| output: true

// DCT-II 2D ortogonal (separável) — usar mm::dct2 do morph.hpp
// IDCT-II 2D ortogonal — usar mm::idct2 do morph.hpp

int main() {
// [pdi:state-io] auto-generated — do not edit by hand
mm::Image img_gray = mm::_read_state("tmp/state/img_gray_20.png");
// [pdi:state-io:end]

    // ── Bloco 8×8 centralizado da imagem ─────────────────────────────────────────
    int cy = img_gray.h / 2;
    int cx = img_gray.w / 2;

    // Extrair bloco 8×8
    cv::Mat img_mat = img_gray;
    cv::Mat bloco_mat = img_mat(cv::Rect(cx, cy, 8, 8)).clone();

    // Converter para float64 e subtrair 128
    cv::Mat bloco_f;
    bloco_mat.convertTo(bloco_f, CV_64F);
    bloco_f -= 128.0;

    // Aplicar DCT2
    cv::Mat C = mm::dct2(bloco_f);

    std::cout << "Coeficientes DCT do bloco 8×8:" << std::endl;

    // Arredondar e imprimir
    cv::Mat C_rounded;
    cv::Mat C_abs, C_int;
    cv::normalize(C, C_abs, 0, 255, cv::NORM_MINMAX);
    C.convertTo(C_int, CV_32S, 1.0, 0.5); // arredondamento via conversão
    std::cout << C_int << std::endl;

    double dc_energy = C.at<double>(0, 0) * C.at<double>(0, 0);
    double total_energy = 0.0;
    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            total_energy += C.at<double>(i, j) * C.at<double>(i, j);
        }
    }

    std::cout << "\nEnergia DC     : " << dc_energy << std::endl;
    std::cout << "Energia total  : " << total_energy << std::endl;
    std::cout << "Fração no DC   : " << (dc_energy / total_energy) * 100.0 << "% ← concentração de energia" << std::endl;

    // ── Reconstrução progressiva ──────────────────────────────────────────────────
    std::vector<mm::Image> imgs_rec;
    std::vector<std::string> titles_rec;

    // Bloco original
    cv::Mat bloco_original;
    bloco_mat.convertTo(bloco_original, CV_8U);
    imgs_rec.push_back(mm::Image(bloco_original));
    titles_rec.push_back("Bloco original\n(8×8 pixels)");

    // Diferentes quantidades de coeficientes
    int keeps[] = {1, 4, 10, 20, 40, 64};

    for (int keep : keeps) {
        cv::Mat C_trunc = cv::Mat::zeros(8, 8, CV_64F);

        // Zerar coeficientes mantendo apenas os "keep" primeiros por ordem zigzag
        // (ordem por soma u+v, que é uma aproximação da ordem zigzag)
        std::vector<std::pair<int, int>> indices;
        for (int u = 0; u < 8; u++) {
            for (int v = 0; v < 8; v++) {
                indices.push_back({u, v});
            }
        }

        // Ordenar por soma u+v (frequência crescente)
        std::sort(indices.begin(), indices.end(), 
                  [](const std::pair<int, int>& a, const std::pair<int, int>& b) {
                      return (a.first + a.second) < (b.first + b.second);
                  });

        for (int k = 0; k < keep && k < 64; k++) {
            int u = indices[k].first;
            int v = indices[k].second;
            C_trunc.at<double>(u, v) = C.at<double>(u, v);
        }

        // Reconstruir
        cv::Mat rec_f = mm::idct2(C_trunc);
        rec_f += 128.0;

        // Recortar para [0, 255] e converter para uint8
        cv::Mat rec;
        cv::threshold(rec_f, rec_f, 255, 255, cv::THRESH_TRUNC);
        cv::threshold(rec_f, rec_f, 0, 0, cv::THRESH_TOZERO);
        rec_f.convertTo(rec, CV_8U);

        imgs_rec.push_back(mm::Image(rec));

        char buf[64];
        snprintf(buf, sizeof(buf), "%d coef.\n(%.0f%% do total)", 
                 keep, (keep / 64.0) * 100.0);
        titles_rec.push_back(buf);
    }

    mm::show(imgs_rec, MM_OUT, titles_rec, 4);

    return 0;
}

Overwriting tmp/fig_05_dct_bloco.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_dct_bloco.cpp -o tmp/fig_05_dct_bloco -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_dct_bloco \
  && test -f "tmp/fig_05_dct_bloco.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_dct_bloco.png"

Coeficientes DCT do bloco 8×8:
[192, -47, 1, 9, 0, 0, 0, 0;
 -95, 55, 7, -10, 0, 1, 1, 0;
 14, -14, 9, 1, 1, 0, 1, 0;
 0, 1, -11, 1, 1, 1, 1, 0;
 10, -10, 0, 1, 1, 0, 1, 0;
 0, 1, 0, 1, 0, 0, 1, 1;
 0, 1, 0, 0, 0, 1, 1, 0;
 1, 0, 0, 0, 1, 1, 1, 1]

Energia DC     : 36816
Energia total  : 52253
Fração no DC   : 70.4572% ← concentração de energia
[1] Bloco original
(8×8 pixels)
[2] 1 coef.
(2% do total)
[3] 4 coef.
(6% do total)
[4] 10 coef.
(16% do total)
[5] 20 coef.
(31% do total)
[6] 40 coef.
(62% do total)
[7] 64 coef.
(100% do total)

try:
    mm.show(mm.read("tmp/fig_05_dct_bloco.png"), figsize=(12, 7))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_dct_bloco.png (ver a versao Python)")

Figura 5.24: DCT 2D em bloco 8×8: coeficientes e reconstrução progressiva.

5.7.5 El Pipeline de Compresión JPEG

El estándar JPEG opera dividiendo la imagen en bloques disjuntos de \(8 \times 8\) píxeles, procesados mediante una secuencia de transformaciones espaciales, perceptuales y estadísticas. El pipeline completo de codificación se estructura en seis etapas principales:

\[ \text{RGB} \xrightarrow{\text{(1) } YC_bC_r} \xrightarrow{\text{(2) Submuestreo}} \xrightarrow{\text{(3) Bloques } 8 \times 8} \xrightarrow{\text{(4) DCT}} \xrightarrow{\text{(5) Cuantificación}} \xrightarrow{\text{(6) Codificación Entrópica}} \]

La Tabla 5.7 detalla la función analítica y el fundamento perceptual que justifica cada una de estas etapas.

Tabla 5.7: Etapas del pipeline de compresión JPEG y sus respectivos fundamentos de diseño.

Etapa	Operación	Fundamento Perceptual y Estadístico
1	Conversión \(RGB \rightarrow YC_bC_r\)	Separa la luminancia (\(Y\)) de la crominancia (\(C_b, C_r\)). El sistema visual humano (SVH) presenta mayor sensibilidad a variaciones de brillo que de color.
2	Submuestreo de crominancia (ej: 4:2:0)	Reduce la resolución espacial de los canales de color a la mitad, descartando datos redundantes con impacto visual despreciable.
3–4	Centralización y aplicación de la DCT \(8 \times 8\)	Traslada los píxeles al intervalo \([-128, 127]\) y compacta la energía espectral del bloque en los coeficientes de baja frecuencia.
5	Cuantificación lineal selectiva	Divide cada coeficiente \(C(u,v)\) por el elemento correspondiente de la matriz \(Q(u,v)\), aplicando redondeo entero. Constituye la principal fuente de compresión con pérdida.
6	Barrido en zigzag y codificación	Ordena los coeficientes cuantificados para maximizar secuencias nulas consecutivas, optimizando la codificación por longitud de corrida (RLE) y la codificación de Huffman.

La matriz de cuantificación \(Q(u,v)\) es el mecanismo central de control del compromiso entre tasa de compresión y calidad visual. En el algoritmo práctico de la Figura 5.25, el factor de calidad estipulado por el usuario (escala de 1 a 100) se convierte en un escalar que parametriza la severidad de la matriz \(Q\). Valores reducidos de calidad expanden los divisores de \(Q(u,v)\), forzando el truncamiento masivo de los coeficientes AC a cero. Cuando esta eliminación es excesiva, la discontinuidad en las fronteras de los bloques adyacentes no se atenúa en la reconstrucción, generando los denominados artefactos de bloque (blocking artifacts).

La Lógica del Barrido en Zigzag

La eficiencia del codificador entrópico posterior a la cuantificación depende directamente de la ordenación de los datos. Como la DCT concentra la energía vital en el vértice superior izquierdo de la matriz (bajas frecuencias) y empuja los coeficientes nulos hacia las extremidades opuestas, la lectura lineal por filas o columnas fragmentaría las secuencias de ceros.

La ordenación en zigzag soluciona esta limitación al recorrer la matriz diagonalmente en orden creciente de frecuencia espacial. Este mapeo agrupa los coeficientes significativos al inicio del vector y concentra los coeficientes nulos en una única secuencia continua al final del arreglo, permitiendo que el algoritmo RLE codifique grandes bloques de datos de forma compacta y eficiente.

Nota¿Qué es RLE?

RLE (Run-Length Encoding) es una técnica de compresión sin pérdidas que codifica secuencias consecutivas de valores idénticos — especialmente ceros — como un par (recuento, valor). En JPEG, después del barrido en zigzag, los coeficientes cuantificados se organizan de modo que los ceros se concentren al final del vector. El RLE entonces comprime esa larga corrida de ceros con extrema eficiencia, optimizando el almacenamiento y la transmisión de la imagen comprimida.

%%writefile tmp/fig_05_jpeg_pipeline.cpp
#define MM_OUT "tmp/fig_05_jpeg_pipeline.png"
#include <opencv2/opencv.hpp>
#include <string>
#include <vector>
#include <cstdio>
#include "morph.hpp"

//| label: fig-05-jpeg-pipeline
//| fig-cap: "*Pipeline* JPEG simplificado aplicado à imagem clássica do *Cameraman*: DCT em blocos 8×8, quantização com diferentes fatores de qualidade e reconstrução via IDCT. Os artefatos de bloco (*blocking artifacts*) tornam-se visualmente evidentes em fatores de qualidade reduzidos ($Q=10$ e $Q=25$)."
//| echo: true
//| output: true

int main() {
    // Pipeline JPEG (DCT 8x8 -> quantizacao -> IDCT) via mm.jpegCompress,
    // na imagem classica do Cameraman (asset do capitulo) reduzida a 256x256.
    cv::Mat img_src_mat = cv::imread("imagens/cameraman.png", cv::IMREAD_GRAYSCALE);
    if (img_src_mat.empty()) {
        return 1;
    }
    cv::Mat img_resized;
    cv::resize(img_src_mat, img_resized, cv::Size(256, 256));
    mm::Image img_src = mm::gray(img_resized);

    std::vector<mm::Image> imgs;
    imgs.push_back(img_src);
    std::vector<std::string> titles;
    titles.push_back("Original (Cameraman)");

    // Store PSNR values before pushing to titles
    std::vector<double> psnr_vals;
    for (int q : {10, 25, 50, 75, 90}) {
        mm::Image rec = mm::jpegCompress(img_src, q);
        imgs.push_back(rec);
        double psnr_val = mm::psnr(img_src, rec);
        psnr_vals.push_back(psnr_val);
        char title[128];
        std::snprintf(title, sizeof(title), "Q=%d (PSNR=%.1f dB)", q, psnr_val);
        titles.push_back(std::string(title));
    }

    mm::show(imgs, MM_OUT, titles, 3);

    return 0;
}

Overwriting tmp/fig_05_jpeg_pipeline.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_jpeg_pipeline.cpp -o tmp/fig_05_jpeg_pipeline -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_jpeg_pipeline \
  && test -f "tmp/fig_05_jpeg_pipeline.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_jpeg_pipeline.png"

[1] Original (Cameraman)
[2] Q=10 (PSNR=28.0 dB)
[3] Q=25 (PSNR=30.7 dB)
[4] Q=50 (PSNR=32.8 dB)
[5] Q=75 (PSNR=35.2 dB)
[6] Q=90 (PSNR=40.0 dB)

try:
    mm.show(mm.read("tmp/fig_05_jpeg_pipeline.png"), figsize=(14, 10))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_jpeg_pipeline.png (ver a versao Python)")

Figura 5.25: Pipeline JPEG simplificado aplicado à imagem clássica do Cameraman: DCT em blocos 8×8, quantização com diferentes fatores de qualidade e reconstrução via IDCT. Os artefatos de bloco (blocking artifacts) tornam-se visualmente evidentes em fatores de qualidade reduzidos (\(Q=10\) e \(Q=25\)).

5.7.6 Simulador Interactivo: Cuantización DCT

El simulador de la Figura 5.26 permite explorar el impacto del proceso de cuantización sobre un bloque \(8 \times 8\) extraído de una imagen real, sintetizando en tiempo real los siguientes componentes:

• Bloque original y reconstruido: Representación directa de los píxeles en el dominio espacial en escala de grises [0, 255].
• Coeficientes DCT: Distribución de la energía mapeada de forma logarítmica en un gradiente cromático, evidenciando la concentración de intensidad en el vértice superior izquierdo (bajas frecuencias).
• Coeficientes cuantizados: Exhibición de los valores enteros resultantes de la división por la matriz \(Q(u,v)\), haciendo visualmente explícita la aparición masiva de coeficientes nulos (en tonos oscuros) conforme el factor de calidad se reduce.
• Métricas de compresión: Panel de monitoreo que cuantifica el Error Cuadrático Medio (MSE), el número de coeficientes preservados y el volumen de ceros generados para la codificación entrópica.

⊞ Simulador: Cuantización DCT-JPEG (bloque 8×8) bloques 8×8

Calidad

50

Coef. ≠ 0

–

Ceros

–

Error MSE

–

Bloque Original (8×8)

Coef. DCT (abs, log)

Coef. Cuantizados

Bloque Reconstruido

Calidad JPEG: 50

Figura 5.26: Simulador interactivo de compresión DCT-JPEG: ajuste el factor de calidad y visualice en tiempo real los coeficientes anulados, el bloque reconstruido y el error de cuantización.

5.8 Comparación de Formatos de Imagen

La elección de un formato de almacenamiento digital impacta directamente en el compromiso entre calidad visual, tamaño de archivo y costo computacional de decodificación. Los tres formatos de mayor relevancia para arquitecturas web y sistemas de computación visual son JPEG, PNG y WebP.

5.8.1 Características de los Formatos

La Tabla 5.8 sintetiza las propiedades estructurales de los principales formatos de imagen rasterizados.

Tabla 5.8: Comparación estructural entre los principales formatos de imagen rasterizados.

Característica	JPEG	PNG	WebP
Compresión	Con pérdida	Sin pérdida	Con y sin pérdida.
Transparencia (canal alfa)	No	Sí	Sí.
Soporte de animación	No	Limitado (APNG)	Sí.
Algoritmo base	DCT + Huffman	DEFLATE (LZ77 + Huffman)	VP8 / VP8L.
Mejor para	Fotografía	Gráficos, texto e iconos	Uso universal en entorno Web.
Peor para	Texto y bordes nítidos	Imágenes fotográficas complejas	Compatibilidad heredada.

5.8.2 Métricas de Evaluación de Calidad

Dos métricas objetivas son ampliamente adoptadas para cuantificar la distorsión introducida por procesos de compresión:

Pico de la Relación Señal-Ruido (PSNR, Peak Signal-to-Noise Ratio): \[ \text{PSNR} = 10\,\log_{10}\!\left(\frac{L^2}{\text{MSE}}\right) \quad [\text{dB}] \tag{5.10}\]

donde \(L = 255\) para imágenes cuantizadas en 8 bits y \(\text{MSE}\) representa el Error Cuadrático Medio (Mean Squared Error). Valores de PSNR superiores a 40 dB indican excelente fidelidad; entre 30 dB y 40 dB representan buena calidad; y valores inferiores a 30 dB corresponden a degradaciones visuales fácilmente perceptibles.

Índice de Similitud Estructural (SSIM, Structural Similarity Index): \[ \text{SSIM}(f,g) = \frac{(2\mu_f\mu_g + c_1)(2\sigma_{fg} + c_2)}{(\mu_f^2+\mu_g^2+c_1)(\sigma_f^2+\sigma_g^2+c_2)} \tag{5.11}\]

El SSIM evalúa ventanas locales de la imagen basándose en tres componentes complementarios: luminancia (\(\mu_f, \mu_g\)), contraste (\(\sigma_f, \sigma_g\)) y estructura (\(\sigma_{fg}\)), ponderados por constantes de estabilidad \(c_1\) y \(c_2\). El índice varía en el intervalo \([-1, 1]\), donde la unidad representa la identidad perfecta. A diferencia del PSNR, el SSIM considera la organización espacial de los errores, alineándose con la percepción del sistema visual humano (SVH).

NotaPSNR vs SSIM: Aplicación de Métricas Perceptuales

El PSNR posee una formulación matemática simple y un bajo costo computacional; sin embargo, tiende a sobreestimar la calidad en imágenes con distorsiones localizadas o subestimarla en variaciones globales de brillo toleradas por el observador. El SSIM modela con mayor fidelidad la percepción biológica, pero exige un mayor esfuerzo de procesamiento. Para análisis rigurosos de codificadores, se recomienda reportar ambas métricas estadísticas en carácter complementario.

5.8.3 Inspección Visual: Naturaleza de los Artefactos de Compresión

La naturaleza matemática del codificador determina el tipo de degradación introducida en tasas de bits reducidas. Como se ilustra en la Figura 5.27, la compresión agresiva mediante DCT en el estándar JPEG segmenta la imagen en mallas rígidas, generando los artefactos de bloque (blocking artifacts). En contrapartida, los algoritmos basados en codificación predictiva o representaciones sometidas a transformadas espaciales avanzadas (como WebP y JPEG 2000) eliminan las discontinuidades de bloque, pero introducen pérdida de textura fina y desenfoques característicos alrededor de bordes de alto contraste.

%%writefile tmp/fig_05_zoom_artefatos.cpp
#define MM_OUT "tmp/fig_05_zoom_artefatos.png"
//| label: fig-05-zoom-artefatos
//| fig-cap: "Análise comparativa de artefatos de compressão sob fator de qualidade reduzido ($Q=10$). À esquerda, observa-se o artefato de bloco característico da discretização por DCT no JPEG. À direita, evidencia-se o efeito de atenuação e suavização de bordas intrínseco ao padrão WebP."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include "morph.hpp"
#include <filesystem>

cv::Mat zoom(const cv::Mat& img) {
    // Recorta a região 120:200, 150:230 e redimensiona para 320x320 com interpolação mais próxima
    cv::Mat crop = img(cv::Rect(150, 120, 80, 80));
    cv::Mat resized;
    cv::resize(crop, resized, cv::Size(320, 320), 0, 0, cv::INTER_NEAREST);
    return resized;
}

int main() {
    // Lê a imagem, converte para escala de cinza e redimensiona para 256x256
    cv::Mat img_src_m = cv::imread("imagens/cameraman.png", cv::IMREAD_GRAYSCALE);
    cv::Mat img_src_resized;
    cv::resize(img_src_m, img_src_resized, cv::Size(256, 256));
    mm::Image img_src = img_src_resized;  // converte para mm::Image para usar mm::gray se necessário, mas já é gray

    // Salva com qualidade JPEG 10 e WebP 10
    std::filesystem::create_directories("tmp");
    cv::imwrite("tmp/zoom_q10.jpg", img_src_resized, {cv::IMWRITE_JPEG_QUALITY, 10});
    cv::imwrite("tmp/zoom_q10.webp", img_src_resized, {cv::IMWRITE_WEBP_QUALITY, 10});

    // Aplica zoom na imagem original, JPEG e WebP
    cv::Mat z_orig = zoom(img_src_resized);
    cv::Mat z_jpeg = zoom(cv::imread("tmp/zoom_q10.jpg", cv::IMREAD_GRAYSCALE));
    cv::Mat z_webp = zoom(cv::imread("tmp/zoom_q10.webp", cv::IMREAD_GRAYSCALE));

    // Exibe as três imagens em uma grade
    mm::show(
        std::vector<mm::Image>{mm::Image(z_orig), mm::Image(z_jpeg), mm::Image(z_webp)},
        MM_OUT,
        {"Zoom Original", "JPEG Q=10 (Artefato de Bloco)", "WebP Q=10 (Suavizacao)"},
        3
    );

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(z_orig, "tmp/fig_05_zoom_artefatos_0.png");
mm::write(z_jpeg, "tmp/fig_05_zoom_artefatos_1.png");
mm::write(z_webp, "tmp/fig_05_zoom_artefatos_2.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_zoom_artefatos.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_zoom_artefatos.cpp -o tmp/fig_05_zoom_artefatos -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_zoom_artefatos \
  && test -f "tmp/fig_05_zoom_artefatos.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_zoom_artefatos.png"

[1] Zoom Original
[2] JPEG Q=10 (Artefato de Bloco)
[3] WebP Q=10 (Suavizacao)

try:
    mm.show(
        [
            mm.read("tmp/fig_05_zoom_artefatos_0.png"),
            mm.read("tmp/fig_05_zoom_artefatos_1.png"),
            mm.read("tmp/fig_05_zoom_artefatos_2.png"),
        ],
        titles=[
            'Zoom Original',
            'JPEG Q=10 (Artefato de Bloco)',
            'WebP Q=10 (Suavizacao)',
        ],
        cols=3,
        figsize=(14, 5),
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_zoom_artefatos_0.png (ver a versao Python)")

Figura 5.27: Análise comparativa de artefatos de compressão sob fator de qualidade reduzido (\(Q=10\)). À esquerda, observa-se o artefato de bloco característico da discretização por DCT no JPEG. À direita, evidencia-se o efeito de atenuação e suavização de bordas intrínseco ao padrão WebP.

5.8.4 Evaluación Cuantitativa y Espacial de la Compresión

La validación de los algoritmos de compresión con pérdida exige un análisis que correlacione el costo de almacenamiento con la fidelidad de la señal reconstruida. Esta evaluación se realiza de manera complementaria mediante curvas de rendimiento global y mediante el mapeo local de las distorsiones inducidas por los codificadores.

5.8.4.1 Curvas de Tasa-Distorsión

La Figura 5.28 presenta la evaluación empírica del pipeline JPEG y WebP mediante curvas de tasa-distorsión, que monitorean la ganancia de compresión (tamaño del archivo en KB) en función del PSNR. El formato PNG actúa como línea de base ideal (\(\text{PSNR} = \infty\)), ya que su naturaleza lossless impide cualquier degradación, aunque requiere un volumen de datos sustancialmente mayor.

El análisis de las curvas demuestra la superioridad y la eficiencia del estándar WebP sobre el JPEG tradicional: para alcanzar un mismo nivel de fidelidad matemática (como el rango de calidad excelente, donde \(\text{PSNR} > 40\text{ dB}\)), el codificador WebP genera archivos significativamente más pequeños. Este comportamiento refleja el impacto práctico de la evolución de los algoritmos en la optimización de los sistemas de transmisión y almacenamiento digital.

%%writefile tmp/fig_05_formatos_comparacao.cpp
#define MM_OUT "tmp/fig_05_formatos_comparacao.png"
//| label: fig-05-formatos-comparacao
//| fig-cap: "Curva taxa-distorção: PSNR vs tamanho de arquivo para JPEG, WebP e PNG aplicada à imagem do *Cameraman*."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <string>
#include <vector>
#include <filesystem>
#include <cstdio>
#include "morph.hpp"

int main() {
    // Criar diretório temporário se não existir
    std::filesystem::create_directories("tmp");

    // Carregar e redimensionar a imagem do Cameraman para 256x256
    cv::Mat img_full = cv::imread("imagens/cameraman.png", cv::IMREAD_GRAYSCALE);
    cv::Mat img_resized;
    cv::resize(img_full, img_resized, cv::Size(256, 256));
    mm::Image src = mm::gray(img_resized);

    // Loop para JPEG com diferentes qualidades
    std::vector<double> jpeg_kb, jpeg_psnr;
    std::vector<int> jpeg_qualities = {10, 20, 30, 40, 50, 60, 70, 80, 90, 95};
    for (int q : jpeg_qualities) {
        std::string p = "tmp/fmt_q" + std::to_string(q) + ".jpg";
        cv::Mat src_mat = src;  // converter mm::Image para cv::Mat
        std::vector<int> params = {cv::IMWRITE_JPEG_QUALITY, q};
        cv::imwrite(p, src_mat, params);
        cv::Mat rec = cv::imread(p, cv::IMREAD_GRAYSCALE);
        mm::Image rec_img = rec;
        jpeg_kb.push_back((double)std::filesystem::file_size(p) / 1024.0);
        jpeg_psnr.push_back(mm::psnr(src, rec_img));
    }

    // Loop para WebP com diferentes qualidades
    std::vector<double> webp_kb, webp_psnr;
    std::vector<int> webp_qualities = {30, 50, 70, 85, 95};
    for (int q : webp_qualities) {
        std::string p = "tmp/fmt_w" + std::to_string(q) + ".webp";
        cv::Mat src_mat = src;  // converter mm::Image para cv::Mat
        std::vector<int> params = {cv::IMWRITE_WEBP_QUALITY, q};
        cv::imwrite(p, src_mat, params);
        cv::Mat rec = cv::imread(p, cv::IMREAD_GRAYSCALE);
        mm::Image rec_img = rec;
        webp_kb.push_back((double)std::filesystem::file_size(p) / 1024.0);
        webp_psnr.push_back(mm::psnr(src, rec_img));
    }

    // PNG sem perda
    std::string p_png = "tmp/fmt.png";
    cv::Mat src_mat_png = src;  // converter mm::Image para cv::Mat
    std::vector<int> params_png = {cv::IMWRITE_PNG_COMPRESSION, 9};
    cv::imwrite(p_png, src_mat_png, params_png);
    double png_kb = (double)std::filesystem::file_size(p_png) / 1024.0;

    // Criar o gráfico de linha comparando JPEG e WebP
    std::string title = "Curva Taxa-Distorcao (PNG sem perda: " + 
                        std::to_string(png_kb).substr(0, std::to_string(png_kb).find(".") + 2) + 
                        " KB)";

    // Preparar vetores como vetores de vetores (uma curva por elemento)
    std::vector<std::vector<double>> xs = {jpeg_kb, webp_kb};
    std::vector<std::vector<double>> ys = {jpeg_psnr, webp_psnr};
    std::vector<std::string> labels = {"JPEG", "WebP"};

    mm::Image chart = mm::lineChart(xs, ys, labels, {}, title, 
                                    "Tamanho do arquivo (KB)", "PSNR (dB)");

    // Exibir o resultado
    std::vector<mm::Image> charts = {chart};
    std::vector<std::string> chart_titles = {"JPEG x WebP x PNG"};
    mm::show(charts, MM_OUT, chart_titles, 1);

    
// [pdi:panel-io] auto-generated — do not edit by hand
std::filesystem::create_directories("tmp");
mm::write(chart, "tmp/fig_05_formatos_comparacao_0.png");
// [pdi:panel-io:end]
return 0;
}

Overwriting tmp/fig_05_formatos_comparacao.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_formatos_comparacao.cpp -o tmp/fig_05_formatos_comparacao -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_formatos_comparacao \
  && test -f "tmp/fig_05_formatos_comparacao.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_formatos_comparacao.png"

[1] JPEG x WebP x PNG

try:
    mm.show(
        [
            mm.read("tmp/fig_05_formatos_comparacao_0.png"),
        ],
        titles=[
            'JPEG x WebP x PNG',
        ],
        cols=1,
    )
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_formatos_comparacao_0.png (ver a versao Python)")

Figura 5.28: Curva taxa-distorção: PSNR vs tamanho de arquivo para JPEG, WebP e PNG aplicada à imagem do Cameraman.

NotaTamaño original de la imagen

La imagen Cameraman (\(256 \times 256\) píxeles en escala de grises) ocupa 64 KB en formato bruto (sin compresión). Como referencia, el PNG lossless comprime ese volumen a 36,2 KB — evidenciando que la compresión sin pérdidas ya reduce significativamente el almacenamiento para imágenes con regiones homogéneas. En contrapartida, los formatos con pérdida (JPEG y WebP) alcanzan tamaños aún menores: el JPEG con calidad 95 ocupa 22,3 KB (PSNR ≈ 45 dB), mientras que el WebP con calidad 90 alcanza 12,5 KB con PSNR equivalente, demostrando su superioridad en eficiencia de compresión.

5.8.4.2 Mapeo Espacial de Errores y Correlación Perceptual

Aunque el PSNR ofrece un indicativo numérico rápido, las métricas globales no logran discriminar cómo se distribuye geométricamente la pérdida de información sobre la imagen. La Figura 5.29 soluciona esta limitación al asociar las reconstrucciones en diferentes calidades con sus respectivos mapas de error absoluto y con el SSIM.

Los mapas residuales — obtenidos mediante la diferencia absoluta normalizada entre la imagen original y la comprimida — revelan la firma espacial intrínseca de cada arquitectura de codificación:

• En calidades altas (\(Q=95\) a \(Q=75\)): Las distorsiones se concentran predominantemente alrededor de transiciones abruptas de intensidad (bordes), como resultado del reflejo espectral derivado del descarte de altas frecuencias. El índice SSIM permanece cercano a la unidad, atestiguando la integridad de las estructuras originales.
• En calidades agresivas (\(Q=50\) a \(Q=25\)): El error adopta una estructura de malla ortogonal regularizada. Este patrón geométrico evidencia la aparición de los artefactos de bloque (blocking artifacts), indicando que la cuantización severa ha corrompido la correlación espacial entre bloques adyacentes de \(8 \times 8\) píxeles.

El SSIM captura esta degradación morfológica de manera mucho más sensible que el PSNR, penalizando la puntuación final a medida que la organización estructural y las texturas finas — a las cuales el sistema visual humano es altamente receptivo — son eliminadas por el codificador.

%%writefile tmp/fig_05_ssim_artefatos.cpp
#define MM_OUT "tmp/fig_05_ssim_artefatos.png"
//| label: fig-05-ssim-artefatos
//| fig-cap: "Análise espacial de degradação: imagens reconstruídas e respectivos mapas de erro absoluto normalizados para diferentes fatores de qualidade JPEG."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>
#include "morph.hpp"

// Implementación de SSIM (Wang et al.) con ventana gaussiana
double compute_ssim(const cv::Mat& img1, const cv::Mat& img2) {
    const double C1 = 6.5025, C2 = 58.5225;
    cv::Mat I1, I2;
    img1.convertTo(I1, CV_64F);
    img2.convertTo(I2, CV_64F);

    cv::Mat mu1, mu2;
    cv::GaussianBlur(I1, mu1, cv::Size(11, 11), 1.5);
    cv::GaussianBlur(I2, mu2, cv::Size(11, 11), 1.5);

    cv::Mat mu1_sq = mu1.mul(mu1);
    cv::Mat mu2_sq = mu2.mul(mu2);
    cv::Mat mu1_mu2 = mu1.mul(mu2);

    cv::Mat sigma1_sq, sigma2_sq, sigma12;
    cv::GaussianBlur(I1.mul(I1), sigma1_sq, cv::Size(11, 11), 1.5);
    sigma1_sq -= mu1_sq;
    cv::GaussianBlur(I2.mul(I2), sigma2_sq, cv::Size(11, 11), 1.5);
    sigma2_sq -= mu2_sq;
    cv::GaussianBlur(I1.mul(I2), sigma12, cv::Size(11, 11), 1.5);
    sigma12 -= mu1_mu2;

    cv::Mat ssim_map;
    cv::Mat cs_map;
    cv::Mat temp1 = 2 * mu1_mu2 + C1;
    cv::Mat temp2 = 2 * sigma12 + C2;
    cv::Mat temp3 = mu1_sq + mu2_sq + C1;
    cv::Mat temp4 = sigma1_sq + sigma2_sq + C2;

    cv::divide(temp1.mul(temp2), temp3.mul(temp4), ssim_map);
    cv::divide(temp2, temp4, cs_map);

    cv::Scalar mssim = cv::mean(ssim_map);
    return mssim[0];
}

int main() {
    // Cameraman (activo del capítulo), 256x256 — misma imagen de la pista py.
    cv::Mat img_raw = cv::imread("imagens/cameraman.png", cv::IMREAD_GRAYSCALE);
    if (img_raw.empty()) {
        return 1;
    }
    cv::Mat img_resized;
    cv::resize(img_raw, img_resized, cv::Size(256, 256));
    mm::Image img_gray(img_resized);

    std::vector<mm::Image> imgs_ssim;
    std::vector<std::string> titles_ssim;
    imgs_ssim.push_back(img_gray);
    titles_ssim.push_back("Original");

    int qs[] = {25, 50, 75, 95};
    for (int idx = 0; idx < 4; ++idx) {
        int q = qs[idx];
        mm::Image rec = mm::jpegCompress(img_gray, q);  // DCT 8x8 -> quant -> IDCT
        double psnr_v = mm::psnr(img_gray, rec);

        cv::Mat gray_mat = img_gray;
        cv::Mat rec_mat = rec;

        double ssim_v = compute_ssim(gray_mat, rec_mat);

        cv::Mat diff = cv::abs(gray_mat - rec_mat);
        cv::Mat diff_vis;
        cv::normalize(diff, diff_vis, 0, 255, cv::NORM_MINMAX, CV_8U);

        imgs_ssim.push_back(rec);
        imgs_ssim.push_back(mm::Image(diff_vis));

        char title1[100];
        snprintf(title1, sizeof(title1), "Q=%d (PSNR=%.1fdB | SSIM=%.3f)", q, psnr_v, ssim_v);
        titles_ssim.push_back(title1);

        char title2[100];
        snprintf(title2, sizeof(title2), "Mapa de erro (Q=%d) - bordas e blocagem", q);
        titles_ssim.push_back(title2);
    }

    mm::show(imgs_ssim, MM_OUT, titles_ssim, 3);

    return 0;
}

Overwriting tmp/fig_05_ssim_artefatos.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_ssim_artefatos.cpp -o tmp/fig_05_ssim_artefatos -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_ssim_artefatos \
  && test -f "tmp/fig_05_ssim_artefatos.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_ssim_artefatos.png"

[1] Original
[2] Q=25 (PSNR=30.7dB | SSIM=0.860)
[3] Mapa de erro (Q=25) - bordas e blocagem
[4] Q=50 (PSNR=32.8dB | SSIM=0.905)
[5] Mapa de erro (Q=50) - bordas e blocagem
[6] Q=75 (PSNR=35.2dB | SSIM=0.938)
[7] Mapa de erro (Q=75) - bordas e blocagem
[8] Q=95 (PSNR=44.8dB | SSIM=0.990)
[9] Mapa de erro (Q=95) - bordas e blocagem

try:
    mm.show(mm.read("tmp/fig_05_ssim_artefatos.png"), figsize=(14, 14))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_ssim_artefatos.png (ver a versao Python)")

Figura 5.29: Análise espacial de degradação: imagens reconstruídas e respectivos mapas de erro absoluto normalizados para diferentes fatores de qualidade JPEG.

NotaInterpretando los mapas de error

Los mapas de error presentados fueron normalizados individualmente (cv2.NORM_MINMAX) para maximizar el contraste visual y revelar la estructura espacial de las distorsiones. Esto significa que:

• En Q=95, el error absoluto es del orden de 0.5–1.5 niveles de gris (imperceptible visualmente), pero la normalización lo amplifica a blanco y negro para evidenciar su ubicación en bordes y transiciones.
• En Q=25, el error absoluto es 10–20 veces mayor (5–15 niveles de gris), pero la normalización también lo lleva al mismo intervalo [0, 255].

Por lo tanto, la intensidad del blanco en los mapas NO es comparable entre diferentes calidades — los mapas sirven únicamente para revelar la firma espacial del error (bordes vs bloques), no su magnitud. La magnitud correcta está dada por los valores de PSNR y SSIM, que muestran claramente que Q=95 tiene un error mucho menor que Q=25.

Síntesis — Compresión JPEG

El proceso de compresión en el estándar JPEG se basa en la aplicación combinada de transformaciones espaciales, perceptuales y estadísticas para reducir las redundancias de una imagen. La Tabla 5.9 resume el papel de cada etapa en el pipeline y su respectivo impacto en la reducción de datos.

Tabla 5.9: Síntesis de las etapas del pipeline de compresión JPEG y sus respectivos impactos.

Etapa	Operación Analítica	Mecanismo de Ganancia / Compresión
Conversión \(YC_bC_r\)	Aislamiento de los canales de luminancia y crominancia.	Modela la percepción del SVH, permitiendo tratar el color y el brillo de forma independiente.
Submuestreo 4:2:0	Reducción de la resolución espacial de los canales de color (\(C_b\) y \(C_r\)).	Elimina aproximadamente el 50% de los datos brutos con un impacto visual mínimo.
DCT \(8 \times 8\)	Mapeo del dominio espacial al dominio de frecuencias espaciales.	Compactación de energía, concentrando la información vital en los primeros coeficientes.
Cuantización Lineal	División entera de los coeficientes por una matriz de ponderación \(Q(u,v)\).	Principal fuente de compresión con pérdida; elimina altas frecuencias imperceptibles.
Codificación Entrópica	Aplicación de algoritmos RLE y codificación de Huffman.	Compresión estadística sin pérdida, optimizada por las largas series de coeficientes nulos.

Artefactos de Degradación Característicos

La aplicación de tasas de compresión excesivamente agresivas (factores de calidad reducidos) introduce distorsiones predecibles en la imagen reconstruida, derivadas de las limitaciones matemáticas del modelo:

• Artefactos de bloqueo (blocking artifacts): Discontinuidades geométricas visibles en los límites de los bloques de \(8 \times 8\) píxeles, causadas por la pérdida de correlación espacial tras la cuantización severa de los componentes de CA.
• Efecto de dispersión (ringing): Oscilaciones fantasma o distorsiones de “humo” alrededor de bordes nítidos y de alto contraste, provocadas por la eliminación abrupta de armónicos de alta frecuencia necesarios para reconstruir funciones escalón.
• Pérdida de textura fina: Atenuación de detalles de alta frecuencia y bajo contraste (como céspedes, tejidos o porosidad), lo que hace que regiones originalmente texturizadas adopten un aspecto excesivamente liso u homogeneizado.

5.9 Aplicación Práctica: Eliminación de Ruido mediante Filtrado Híbrido

Reuniendo las técnicas consolidadas a lo largo de este capítulo, se presenta un pipeline completo de restauración de imágenes que combina el análisis espectral en el dominio de la frecuencia con el filtrado adaptativo en el dominio espacial. El objetivo es atenuar un ruido mixto (compuesto por degradación gaussiana e interferencia periódica) preservando al máximo los detalles estructurales de la imagen original.

\[ \text{Imagen Ruidosa} \xrightarrow{\text{FFT2}} \xrightarrow{\text{Filtro Notch Gaussiano}} \xrightarrow{\text{IFFT2}} \xrightarrow{\text{Filtro Bilateral}} \text{Imagen Restaurada} \]

NotaEvaluación Complementaria: PSNR vs. SSIM

El par de métricas estadísticas PSNR y SSIM proporciona una evaluación cualitativa y morfológica complementaria del proceso de restauración:

• PSNR: Penaliza uniformemente la desviación cuadrática media píxel a píxel.
• SSIM: Evalúa la preservación de estructuras locales perceptualmente relevantes (luminancia, contraste y contornos).

En la práctica, existe un compromiso analítico (trade-off) entre reducción de ruido y preservación de detalles: los filtros espaciales excesivamente agresivos atenúan bien el ruido de alta frecuencia, pero degradan texturas finas y suavizan bordes nítidos — lo que reduce simultáneamente tanto el PSNR como el SSIM en relación con la imagen original. El desafío del diseño de filtros es encontrar el punto de equilibrio que maximice ambas métricas, garantizando una restauración fiel y visualmente agradable.

5.9.1 Análisis de Rendimiento y Conclusión del Capítulo

Los resultados numéricos y visuales generados por la Figura 5.30 demuestran la relevancia práctica de asociar diferentes dominios de procesamiento. La inserción simultánea de ruido periódico y estocástico corrompe las propiedades morfológicas de la señal, reduciendo severamente los índices de similitud y la relación señal-ruido de la imagen de referencia.

El aislamiento y la supresión de los picos armónicos en el dominio de la frecuencia mediante la máscara notch eliminan las franjas de interferencia senoidales dispersas sobre el espacio bidimensional. Como se evidencia en los datos impresos de la Figura 5.30, este filtrado quirúrgico promueve un salto inmediato y sustancial en la métrica PSNR. No obstante, el ruido Gaussiano de alta frecuencia permanece activo de forma homogénea en el espectro, lo que exige un enfoque complementario.

La restauración final se consolida en el dominio espacial con la introducción del filtro bilateral. A diferencia de los operadores de paso bajo convencionales (como el Gaussiano o el de media), que suavizarían indiscriminadamente el ruido y los contornos estructurales, el filtrado bilateral calcula pesos ponderados por la proximidad geométrica y por la diferencia de intensidad radiométrica. Este comportamiento adaptativo atenúa las fluctuaciones estocásticas remanentes en las regiones de transición suave y preserva la nitidez de los bordes espaciales.

La convergencia de ambos enfoques resulta en una mejora sustancial y simultánea del PSNR y del SSIM en relación con la imagen ruidosa — aunque los valores finales permanecen inferiores a los de la imagen original (PSNR = \(\infty\), SSIM = 1,0), debido a la pérdida inevitable de información espectral y de textura durante los procesos de filtrado. La atenuación suave (gaussiana) de los picos en el espectro evita los artefactos de ringing, mientras que el filtro bilateral elimina el ruido estocástico residual sin comprometer la nitidez de los bordes. Los resultados confirman la eficacia y la complementariedad práctica de las herramientas de análisis de frecuencia presentadas en este capítulo, demostrando que el filtrado híbrido (frecuencia + espacial) es superior a cualquier enfoque aislado para la restauración de imágenes degradadas por ruido mixto.

%%writefile tmp/fig_05_pipeline_denoising.cpp
#define MM_OUT "tmp/fig_05_pipeline_denoising.png"
//| label: fig-05-pipeline-denoising
//| fig-cap: "*Pipeline* completo de remoção de ruído misto: (1) adição de ruído gaussiano e periódico; (2) identificação de picos de interferência no espectro de frequências; (3) aplicação de máscara *notch* com atenuação gaussiana suave; (4) pós-processamento via filtro bilateral para eliminação do ruído estocástico residual."
//| echo: true
//| output: true

#include <opencv2/opencv.hpp>
#include <vector>
#include <string>
#include <cmath>
#include <random>
#include <iostream>
#include "morph.hpp"

// Función auxiliar: suprimir pico gaussiano en la máscara
static void suprimir_pico_gaussiano(cv::Mat& mask, int cy, int cx, double sigma = 3.0) {
    int H = mask.rows, W = mask.cols;
    double inv2sigma2 = 1.0 / (2.0 * sigma * sigma);
    for (int y = 0; y < H; ++y) {
        for (int x = 0; x < W; ++x) {
            double dy = y - cy;
            double dx = x - cx;
            double notch = std::exp(-(dy * dy + dx * dx) * inv2sigma2);
            mask.at<double>(y, x) *= (1.0 - notch);
        }
    }
}

int main() {
    // Cameraman (asset del capítulo), 256x256 — misma imagen de la trilha py.
    mm::Image img_orig = mm::read("imagens/cameraman.png");
    cv::Mat img_resized;
    cv::resize((cv::Mat)img_orig, img_resized, cv::Size(256, 256));
    mm::Image img_gray = mm::gray(img_resized);
    int h_img = img_gray.h;
    int w_img = img_gray.w;

    // ── 1. Ruido mixto: gaussiano + periódico ───────────────────────────────────
    std::mt19937 gen(42);
    std::normal_distribution<double> dist(0.0, 15.0);

    int u0 = 15, v0 = 10;

    // Ruido gaussiano
    cv::Mat ruido_gauss(h_img, w_img, CV_64F);
    for (int y = 0; y < h_img; ++y) {
        for (int x = 0; x < w_img; ++x) {
            ruido_gauss.at<double>(y, x) = dist(gen);
        }
    }

    // Ruido periódico
    cv::Mat ruido_period(h_img, w_img, CV_64F);
    for (int y = 0; y < h_img; ++y) {
        for (int x = 0; x < w_img; ++x) {
            ruido_period.at<double>(y, x) = 30.0 * std::sin(2.0 * M_PI * (u0 * (double)x / w_img + v0 * (double)y / h_img));
        }
    }

    // Imagen ruidosa
    cv::Mat img_gray64f(h_img, w_img, CV_64F);
    for (int y = 0; y < h_img; ++y) {
        for (int x = 0; x < w_img; ++x) {
            double val = (double)img_gray.data[y * img_gray.w + x] + ruido_gauss.at<double>(y, x) + ruido_period.at<double>(y, x);
            val = std::min(255.0, std::max(0.0, val));
            img_gray64f.at<double>(y, x) = val;
        }
    }
    cv::Mat img_noisy8u;
    img_gray64f.convertTo(img_noisy8u, CV_8U);
    mm::Image img_noisy(img_noisy8u);

    // ── 2. Espectro (log-magnitud) de la imagen ruidosa ──────────────────────────
    mm::Image mag_n = mm::spectrumMag(img_noisy);

    // ── 3. Máscara notch gaussiana en los 4 picos periódicos ──────────────────────
    int cy0 = h_img / 2, cx0 = w_img / 2;
    cv::Mat mascara_notch(h_img, w_img, CV_64F, cv::Scalar(1.0));

    // Picos: (v0,u0), (-v0,-u0), (v0,-u0), (-v0,u0)
    int picos[4][2] = {{v0, u0}, {-v0, -u0}, {v0, -u0}, {-v0, u0}};
    for (int i = 0; i < 4; ++i) {
        suprimir_pico_gaussiano(mascara_notch, cy0 + picos[i][0], cx0 + picos[i][1], 3.0);
    }

    // ── 4. Filtrado: notch en el dominio de la frecuencia + bilateral ───────────────
    mm::Image img_notch = mm::freqFilter(img_noisy, mascara_notch);
    cv::Mat img_notch_cv = img_notch; // convertir a cv::Mat
    cv::Mat img_den_cv;
    cv::bilateralFilter(img_notch_cv, img_den_cv, 7, 25, 7);
    mm::Image img_den(img_den_cv);

    // Máscara visible (0-255)
    cv::Mat mascara_vis64f;
    cv::multiply(mascara_notch, 255.0, mascara_vis64f);
    cv::Mat mascara_vis8u;
    mascara_vis64f.convertTo(mascara_vis8u, CV_8U);
    mm::Image mascara_vis(mascara_vis8u);

    // Calcular PSNR
    double psnr_n = mm::psnr(img_gray, img_noisy);
    double psnr_no = mm::psnr(img_gray, img_notch);
    double psnr_d = mm::psnr(img_gray, img_den);

    std::cout << "Ruidosa: PSNR=" << psnr_n << " dB | Apos notch: " << psnr_no << " dB | Notch+bilateral: " << psnr_d << " dB" << std::endl;

    // Mostrar resultados
    mm::show(
        std::vector<mm::Image>{img_gray, img_noisy, mag_n, mascara_vis, img_notch, img_den},
        MM_OUT,
        std::vector<std::string>{
            "Original",
            "Ruidosa (PSNR=" + std::to_string(psnr_n).substr(0, 4) + " dB)",
            "Espectro (picos visiveis)",
            "Mascara notch (gaussiana)",
            "Apos notch (PSNR=" + std::to_string(psnr_no).substr(0, 4) + " dB)",
            "Notch + bilateral (PSNR=" + std::to_string(psnr_d).substr(0, 4) + " dB)"
        },
        6
    );

    return 0;
}

Overwriting tmp/fig_05_pipeline_denoising.cpp

!g++ -I. -std=c++17 -DMM_USE_OPENCV -I/usr/include/opencv4 tmp/fig_05_pipeline_denoising.cpp -o tmp/fig_05_pipeline_denoising -lopencv_stitching -lopencv_alphamat -lopencv_aruco -lopencv_barcode -lopencv_bgsegm -lopencv_bioinspired -lopencv_ccalib -lopencv_dnn_objdetect -lopencv_dnn_superres -lopencv_dpm -lopencv_face -lopencv_freetype -lopencv_fuzzy -lopencv_hdf -lopencv_hfs -lopencv_img_hash -lopencv_intensity_transform -lopencv_line_descriptor -lopencv_mcc -lopencv_quality -lopencv_rapid -lopencv_reg -lopencv_rgbd -lopencv_saliency -lopencv_shape -lopencv_stereo -lopencv_structured_light -lopencv_phase_unwrapping -lopencv_superres -lopencv_optflow -lopencv_surface_matching -lopencv_tracking -lopencv_highgui -lopencv_datasets -lopencv_text -lopencv_plot -lopencv_ml -lopencv_videostab -lopencv_videoio -lopencv_viz -lopencv_wechat_qrcode -lopencv_ximgproc -lopencv_video -lopencv_xobjdetect -lopencv_objdetect -lopencv_calib3d -lopencv_imgcodecs -lopencv_features2d -lopencv_dnn -lopencv_flann -lopencv_xphoto -lopencv_photo -lopencv_imgproc -lopencv_core \
  && ./tmp/fig_05_pipeline_denoising \
  && test -f "tmp/fig_05_pipeline_denoising.png" \
  || echo "⚠ mm::show não gravou tmp/fig_05_pipeline_denoising.png"

Ruidosa: PSNR=20.1968 dB | Apos notch: 22.3552 dB | Notch+bilateral: 23.5442 dB
[1] Original
[2] Ruidosa (PSNR=20.1 dB)
[3] Espectro (picos visiveis)
[4] Mascara notch (gaussiana)
[5] Apos notch (PSNR=22.3 dB)
[6] Notch + bilateral (PSNR=23.5 dB)

try:
    mm.show(mm.read("tmp/fig_05_pipeline_denoising.png"), figsize=(20, 4))
except Exception as _e:
    print("figura indisponivel nesta trilha (C++): " + repr(_e) + " tmp/fig_05_pipeline_denoising.png (ver a versao Python)")

Figura 5.30: Pipeline completo de remoção de ruído misto: (1) adição de ruído gaussiano e periódico; (2) identificação de picos de interferência no espectro de frequências; (3) aplicação de máscara notch com atenuação gaussiana suave; (4) pós-processamento via filtro bilateral para eliminação do ruído estocástico residual.

5.10 Resumen del Capítulo

La transición del dominio espacial al dominio de la frecuencia revela la distribución espectral de energía de la imagen, estableciendo la base analítica para el filtrado avanzado, la restauración y la compresión de datos. La articulación estructural de estos conceptos se sintetiza en el mapa conceptual de la Figura 5.31.

Figura 5.31: Mapa conceptual de las transformaciones y propiedades en el dominio de la frecuencia.

Fundamentos Esenciales

• DFT y Percepción Visual: El espectro descompone la imagen en componentes armónicas. La fase retiene la inteligibilidad geométrica de la escena y la localización de contornos, mientras que la magnitud dicta la distribución de contraste y las amplitudes globales.
• Eficiencia Algorítmica: El Teorema de la Convolución hace viable el procesamiento de máscaras de gran escala en el dominio de la frecuencia mediante FFT, reduciendo la complejidad computacional asintótica de \(O(N^2 K^2)\) en el espacio a \(O(N^2 \log N)\).
• Fenómeno de Ringing: Los cortes abruptos en el espectro (Filtros Ideales) generan oscilaciones espaciales no deseadas (fenómeno de Gibbs). La atenuación suave mediante filtros de Butterworth o Gaussianos elimina estas discontinuidades.
• Análisis Multirresolución mediante Wavelets: Superando el carácter puramente global de Fourier, la DWT captura la frecuencia y la localización espacial simultáneamente, fundamentando el estándar JPEG 2000 y subsidiando representaciones jerárquicas análogas a las extracciones de características en Redes Neuronales Convolucionales (CNNs).
• Compresión Perceptual (DCT): El pipeline JPEG explora las limitaciones de contraste del sistema visual humano en altas frecuencias espaciales. La DCT aísla la energía de bloques \(8 \times 8\), permitiendo que la cuantización descarte coeficientes AC de detalles finos sin perjuicio perceptual severo.

Próximos Pasos: El Capítulo 6 inaugura la Parte II de la obra, aplicando las herramientas de procesamiento de imágenes en la resolución de problemas reales de inspección industrial. Se explorarán técnicas de segmentación y análisis de formas para la detección automática de fallas en líneas de producción — desde la identificación de defectos superficiales en piezas hasta la lectura QRCode en pruebas, consolidando el puente entre la teoría presentada en la Parte I y las demandas prácticas de la visión computacional.

5.11 🤖 Uso de Gemini Notebook como Tutor Complementario

En esta edición, se incentiva el uso de la plataforma Gemini Notebook como herramienta complementaria de aprendizaje — no como sustituta de la lectura atenta, de la resolución de ejercicios o de la experimentación práctica. Basado en arquitecturas de inteligencia artificial, el sistema utiliza exclusivamente el material didáctico y los documentos proporcionados por el autor como base de conocimiento, asegurando que las respuestas generadas estén conceptualmente alineadas con el contenido programático y con el enfoque pedagógico adoptado a lo largo de esta obra.

ImportanteAcceso al Tutor Inteligente

🚀 ACCEDER A Gemini Notebook: CAPÍTULO 05

🌐 Idioma y Lenguaje de Programación

El proyecto de este capítulo en Gemini Notebook fue construido únicamente con el texto en portugués y los ejemplos de código en Python. Si estás estudiando con la edición en inglés o francés, o siguiendo la ruta en C++, las respuestas del tutor pueden no corresponder exactamente con la versión que estás leyendo.

Directrices sobre el Contenido Generado por Inteligencia Artificial

Aunque las herramientas de inteligencia artificial constituyen aliados eficientes en el proceso de aprendizaje y revisión, el contenido generado está sujeto a inconsistencias o imprecisiones técnicas. Por ello, es indispensable la consulta sistemática de libros de texto, artículos científicos y fuentes académicas indexadas para la validación rigurosa de la información. Se recomienda encarecidamente la ejecución y la modificación de los ejemplos prácticos en Python proporcionados en este capítulo como método primario de verificación experimental de los resultados.

5.12 Lista de Ejercicios

1. (10%) Implementación Directa de la DFT 2D: Implemente analíticamente la Transformada Discreta de Fourier 2D (DFT) sin la ayuda de funciones nativas de bibliotecas (como np.fft.fft2), utilizando estrictamente la formulación matemática definida en la Ecuación 5.1 para una matriz de dimensiones \(16 \times 16\). Realice la validación numérica comparando los coeficientes generados con los resultados de la función np.fft.fft2, asegurándose de que la desviación absoluta máxima sea inferior a \(10^{-8}\). Mida los tiempos de ejecución de ambos métodos y presente una justificación teórica para la disparidad observada en términos de complejidad asintótica.

2. (15%) Supresión de Ruido Periódico: Agregue interferencias sinusoidales con frecuencias espaciales \((u_0, v_0) \in \{(5,10), (20,5), (30,30)\}\) a la imagen de prueba del Cameraman. Para cada escenario de degradación, diseñe una máscara de filtrado notch específica en el dominio de la frecuencia para aislar y atenuar los picos armónicos no deseados. Evalúe cuantitativamente la eficacia del proceso de restauración mediante el cálculo de las métricas de PSNR y SSIM. Discuta analíticamente el compromiso (trade-off) entre la atenuación del ruido sinusoidal y la indeseada atenuación de características estructurales legítimas de la imagen.

3. (15%) Análisis Comparativo de Operadores Pasa-Bajas: Realice un estudio comparativo entre los filtros pasa-bajas Ideal, Gaussiano y Butterworth (con órdenes armónicas \(n = 1, 2, 4\)), parametrizados con frecuencias de corte \(D_0 = 20, 40, 60\) píxeles. Para cada combinación estructural, calcule los índices PSNR y SSIM de la imagen resultante frente a la señal original de referencia. Organice los datos cuantitativos en una tabla estructurada y grafique las curvas unidimensionales de las funciones de transferencia correspondientes a lo largo del perfil horizontal \(H(u, 0)\).

4. (15%) Banco de Filtros Multirresolución de Haar: Desarrolle un script para ejecutar manualmente la descomposición wavelet discreta 2D de primer nivel utilizando la familia Haar. El algoritmo debe calcular los coeficientes de los filtros correspondientes pasa-bajas (\(h\)) y pasa-altas (\(g\)), aplicándolos de forma separable sobre las filas y columnas de la matriz, seguidos por la operación de diezmado (submuestreo espacial por un factor de 2). Valide numéricamente la exactitud de su implementación contrastando las subbandas obtenidas con la salida de la función pywt.dwt2(img, 'haar').

5. (15%) Compresión Dispersa por Umbralización Wavelet: Aplique la técnica de filtrado por umbralización abrupta (hard thresholding) sobre los coeficientes de detalle de la descomposición wavelet, adoptando los umbrales numéricos \(T \in \{5, 10, 20, 40, 80\}\) para las familias Haar, Daubechies (db4) y Symlets (sym4). Después de realizar el proceso de síntesis mediante la transformada inversa (pywt.waverec2), calcule los valores de PSNR y SSIM de cada imagen reconstruida. Identifique y justifique qué combinación de familia wavelet y umbral \(T\) maximiza la similitud estructural.

6. (15%) Construcción de Codificador JPEG Simplificado: Implemente el flujo completo de compresión de datos simulando el estándar JPEG. El flujo debe abarcar: conversión espacial \(RGB \rightarrow YC_bC_r\), submuestreo cromático en la proporción 4:2:0, segmentación de la luminancia en bloques disjuntos de \(8 \times 8\) píxeles, aplicación de la DCT-II 2D ortogonal y cuantización lineal basada en la matriz normalizada de luminancia escalada por factores de calidad deseados. Realice la decodificación inversa y compare cuantitativamente las reconstrucciones con los archivos generados por la función cv2.imencode para los factores de calidad de 20, 50 y 80.

7. (15%) Análisis Perceptual en Contenidos Heterogéneos: Desarrolle una imagen sintética compuesta por tres regiones distintas y de características espectrales contrastantes: una textura fotográfica compleja (representando altas frecuencias estocásticas), un área de texto vectorizado con bordes nítidos (representando transiciones escalón puras) y un gradiente lineal continuo (representando bajas frecuencias homogéneas). Someta esta imagen mixta a los procesos de compresión bajo los formatos JPEG, PNG y WebP. Evalúe e interprete los resultados correlacionando el tamaño final del archivo en disco con las métricas PSNR y SSIM obtenidas, justificando qué formato exhibe el mejor desempeño para señales de naturaleza heterogénea y por qué ocurre esa ventaja en términos de compactación de energía y preservación perceptual.

Referencias del Capítulo

La fundamentación teórica y el desarrollo analítico de los conceptos tratados en este capítulo se basan en las siguientes obras de referencia:

• Gonzalez (2018) — Formulaciones clásicas de Transformadas Discretas de Fourier 2D (DFT), diseño de filtros analíticos en el dominio de la frecuencia, Transformada Discreta de Cosenos (DCT) y principios fundamentales de sistemas de compresión de imágenes.
• Oppenheim (2010) — Teoría formal de señales y sistemas aplicados en el dominio discreto, cubriendo las propiedades matemáticas de la DFT y el modelado analítico del Teorema de la Convolución.
• Mallat (1999) — Fundamentación matemática de la teoría de wavelets, formalización del análisis multirresolución (MRA) y arquitectura de bancos de filtros diádicos.
• Wallace (1991) — Especificación original y aspectos de ingeniería del estándar de compresión ISO/IEC JPEG, con énfasis en los criterios psicovisuales para el diseño de matrices de cuantización DCT.
• Szeliski (2022) — Modelado computacional y caracterización de métricas modernas de fidelidad y calidad perceptual (PSNR y SSIM), así como el análisis comparativo de formatos de imagen rasterizados de alto rendimiento.

---

5.13 💻 Parte Práctica con Ejercicios de Programación

La presente lista de ejercicios de programación (EP) consolida las formulaciones teóricas presentadas a lo largo del Capítulo 5 — Transformadas y Compresión — mediante una ruta práctica aplicada. Los ejercicios se estructuran a partir de matrices de dimensiones reducidas, lo que permite la validación analítica y la inspección manual de cada coeficiente, manteniendo la consistencia metodológica adoptada en los capítulos anteriores.

El encadenamiento de los ejercicios reproduce rigurosamente el flujo conceptual del capítulo: se comienza con la implementación explícita de la Transformada Discreta de Fourier (DFT) a partir de su definición matemática fundamental; se avanza hacia el diseño de filtros pasa-baja y máscaras notch en el dominio de la frecuencia; se aplica la cuantización de coeficientes (núcleo de la compresión con pérdida); y se concluye con la integración de estas etapas en la construcción de un pipeline de compresión JPEG simplificado y en el análisis perceptual de formatos de imagen.

ImportanteDirectrices para la Resolución de los Ejercicios de Programación

En todos los ejercicios de este capítulo, las coordenadas del centro del espectro (origen de las frecuencias espaciales tras la aplicación del desplazamiento fftshift) deben determinarse mediante división entera. Para una matriz con \(L\) filas y \(C\) columnas, la componente de frecuencia nula se localiza en la posición:

\[ (c_y, c_x) = \left( \left\lfloor \frac{L}{2} \right\rfloor, \left\lfloor \frac{C}{2} \right\rfloor \right) \]

Esta convención es rigurosamente idéntica a la adoptada por la función np.fft.fftshift. Además, en todas las etapas que requieran discretización o redondeo numérico (ya sea en la cuantización de coeficientes AC o en la reconstrucción final de píxeles), debe emplearse el redondeo estándar al entero más cercano (round half away from zero), mitigando ambigüedades en valores con fracción exactamente igual a \(0.5\).

🎯 Objetivo de este Cuaderno

El cuaderno permite desarrollar, validar, organizar y probar soluciones de Ejercicios de Programación (EPs) en entornos interactivos, como Colab, con los mismos casos de prueba de Moodle, copiándolos allí solo en el momento de registrar la nota oficial.

Download

Descargue morph.py y testsuite.py ejecutando la celda a continuación:

import os, urllib.request

os.makedirs("tmp/state", exist_ok=True)

url = "https://raw.githubusercontent.com/fzampirolli/pdi-vc/master/morph/config.py"
if not os.path.exists("config.py"):
    urllib.request.urlretrieve(url, "config.py")

import config
config.setup(testsuite=True, cpp=True)
from morph import mm
from testsuite import TestSuite

✅ Entorno listo. Morph: 1.1.9 | OpenCV: 5.0.0 | TestSuite: 1.1.2

Ejecutando las Pruebas

Para evaluar las pruebas, ejecuta TestSuite("EP05_01.extensión").run() en una nueva celda, reemplazando la extensión por la del lenguaje utilizado (.py, .java, .c, .cpp, .js o .r). El sistema descarga los casos de prueba de GitHub, ejecuta el programa y calcula la nota automáticamente.

Para probar código Python directamente, sin guardar archivo, usa run_code(codigo) pasando el código como cadena en una variable codigo:

codigo = """
from morph import mm
# ... tu código aquí ...
"""
TestSuite("EP05_01").run_code(codigo)

5.13.1 EP05_01 🟢 Filtro Pasa-Bajas Ideal por Distancia en el Espectro

En un escáner de documentos antiguo, el sensor capta papel arrugado y textura de fibra junto con el texto — ruido de alta frecuencia que “contamina” el espectro en los bordes. El técnico de mantenimiento no tiene acceso a la imagen original, solo al espectro de magnitud ya calculado por el software del escáner. Su trabajo es simple y quirúrgico: mantener únicamente el círculo central de bajas frecuencias (la estructura global del documento) y borrar todo lo que esté fuera del radio \(D_0\), eliminando la textura fina sin siquiera tocar la imagen espacial.

Este es el Filtro Pasa-Bajas Ideal (LPFI): la operación espectral más directa del capítulo, pero también la que mejor revela la anatomía de un espectro centrado.

5.13.1.1 📋 Directrices de Implementación

1. Dimensiones: Leer los enteros \(L\) (filas) y \(C\) (columnas) del espectro de magnitud — ya proporcionado centrado (equivalente a la salida de np.fft.fftshift).
2. Frecuencia de corte: Leer el entero \(D_0\).
3. Datos: Leer los valores enteros de la matriz de magnitud, fila por fila.
4. Centro del espectro: Calcular \((c_y, c_x) = (L \mathbin{//} 2,\; C \mathbin{//} 2)\).
5. Distancia: Para cada posición \((u,v)\), calcular \[ D(u,v) = \sqrt{(u-c_y)^2 + (v-c_x)^2} \]
6. Máscara ideal: Aplicar \[ H(u,v) = \begin{cases} 1, & D(u,v) \le D_0 \\ 0, & D(u,v) > D_0 \end{cases} \]
7. Filtrado: El valor de salida es \(\text{mag}'(u,v) = \text{mag}(u,v) \cdot H(u,v)\).
8. Salida: Mostrar la matriz filtrada con dimensiones \(L \times C\).

5.13.1.2 📌 Restricciones Computacionales

• Comparación no estricta: el criterio usa \(D(u,v) \le D_0\) (la frontera pertenece al filtro, es decir, se mantiene).
• Tipo: todos los valores de entrada y salida son enteros; la distancia se calcula en punto flotante solo internamente.
• Sin redondeo de magnitud: como la entrada ya es entera y la máscara es binaria (0 o 1), la salida nunca necesita redondeo.

5.13.1.3 🧠 Fundamentación Teórica

Región	Distancia al centro	Efecto del filtro
Centro (\(D \le D_0\))	Bajas frecuencias	Preservadas — estructura global mantenida
Bordes (\(D > D_0\))	Altas frecuencias	Puestas a cero — textura y ruido eliminados
\(D_0\) pequeño	—	Imagen reconstruida quedaría muy borrosa
\(D_0\) grande	—	Poca filtración; casi toda la energía preservada

5.13.1.4 📦 Especificación de Entrada y Salida (VPL)

Entrada:

• Línea 1: Entero \(L\).
• Línea 2: Entero \(C\).
• Línea 3: Entero \(D_0\).
• Líneas siguientes: Elementos enteros de la matriz de magnitud (centrada).

Salida:

• Matriz filtrada en \(L\) filas y \(C\) columnas, separados por espacio.

5.13.1.5 📌 Ejemplos

Entrada	Salida	Observación
3 3 1 10 20 30 40 50 60 70 80 90	0 20 0 40 50 60 0 80 0	Centro \((1,1)\). Las esquinas tienen \(D=\sqrt{2}\approx1.41 > 1\), por lo que se ponen a cero; los vecinos ortogonales tienen \(D=1 \le 1\) y se mantienen.
1 3 0 5 9 7	0 9 0	\(L=1, C=3\): centro en \((0,1)\). Solo la propia posición central (\(D=0\)) sobrevive a \(D_0=0\).

🎮 Simulador EP05_01: Filtro Pasa-Baja Ideal H = (D ≤ D₀) ? 1 : 0

Radio de corte (D₀): 1

Ajusta D₀ y observa qué posiciones del espectro 5×5 sobreviven al filtro.

Espectro Original (Magnitud)

Resultado Filtrado

–

Figura 5.32: Simulador EP05_01: Filtro Pasabajas Ideal en el Espectro

%%writefile EP05_01.cpp
// your solution

Overwriting EP05_01.cpp

TestSuite("EP05_01.cpp").run()

✔️ EP05_01.cases ya existe en casos/
📋 5 caso(s) cargado(s) de casos/EP05_01.cases

🔍 Probando C++: EP05_01.cpp
⚠️ EP05_01.cpp: archivo vacío (menos de 3 líneas). Pruebas omitidas.

5.13.2 EP05_02 🟡 Filtro Notch: Eliminando Picos Periódicos

Una cámara de inspección industrial captura imágenes de placas de circuito, pero la fuente de alimentación de la línea de producción introduce una interferencia eléctrica periódica — un patrón de franjas casi imperceptible a simple vista, pero que aparece en el espectro de Fourier como pares de picos brillantes simétricamente posicionados alrededor del centro. El equipo de visión por computadora no puede reprocesar la captura: necesita localizar y borrar quirúrgicamente esos pares de picos en el espectro, preservando todo el resto de la información útil de la imagen.

Ese es el papel del filtro rechaza-banda notch: a diferencia del pasa-bajas (que afecta una región continua), ataca puntos específicos y sus simétricos, dejando el resto del espectro intacto.

5.13.2.1 📋 Directrices de Implementación

1. Dimensiones: Leer los enteros \(L\) (filas) y \(C\) (columnas) del espectro de magnitud centrado.
2. Datos: Leer los valores enteros de la matriz de magnitud, fila por fila.
3. Picos: Leer el entero \(K\) (cantidad de pares de picos a eliminar).
4. Para cada uno de los \(K\) picos: leer tres enteros \(\Delta v\), \(\Delta u\), \(r\) — desplazamiento vertical, desplazamiento horizontal y radio del notch.
5. Centro del espectro: \((c_y, c_x) = (L \mathbin{//} 2,\; C \mathbin{//} 2)\).
6. Supresión simétrica: para cada pico, poner a cero todas las posiciones \((u,v)\) tales que la distancia al punto \((c_y+\Delta v,\, c_x+\Delta u)\) sea \(\le r\), y también todas las posiciones con distancia \(\le r\) al punto simétrico \((c_y-\Delta v,\, c_x-\Delta u)\).
7. Salida: Mostrar la matriz resultante con dimensiones \(L \times C\).

5.13.2.2 📌 Restricciones Computacionales

• Simetría obligatoria: cada pico informado genera dos discos puestos a cero (el punto y su simétrico respecto al centro) — olvidar el simétrico es el error más común.
• Superposición: si dos discos se superponen, la posición permanece en cero (no hay “suma” ni restauración).
• Comparación no estricta: una posición se pone a cero si \(\text{distancia} \le r\).
• Orden de lectura: los \(K\) picos deben procesarse en el orden en que aparecen en la entrada, pero el resultado final no depende del orden (las operaciones de poner a cero son conmutativas).

5.13.2.3 🧠 Fundamentación Teórica

Concepto	Papel en el filtro notch
Pico en \((\Delta v, \Delta u)\)	Frecuencia de la interferencia periódica detectada visualmente en el espectro
Punto simétrico \((-\Delta v,-\Delta u)\)	Toda DFT de señal real es hermítica: los picos siempre aparecen en pares simétricos al centro
Radio \(r\)	Controla la “anchura” del rechazo — \(r\) grande elimina más energía alrededor del pico, pero también información útil

5.13.2.4 📦 Especificación de Entrada y Salida (VPL)

Entrada:

• Línea 1: Entero \(L\).
• Línea 2: Entero \(C\).
• Líneas siguientes: Elementos enteros de la matriz de magnitud (centrada), \(L\) filas.
• Siguiente línea: Entero \(K\).
• \(K\) líneas siguientes: tres enteros \(\Delta v\), \(\Delta u\), \(r\) (separados por espacios).

Salida:

• Matriz resultante en \(L\) filas y \(C\) columnas, separadas por espacios.

5.13.2.5 📌 Ejemplos

Entrada	Salida	Observación
5 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 1 1 0	1 2 3 4 5 6 0 8 9 10 11 12 13 14 15 16 17 18 0 20 21 22 23 24 25	Centro \((c_y, c_x) = (2, 2)\). El pico informado \((\Delta v, \Delta u) = (1, 1)\) genera el punto \((3, 3)\) (valor 19) y su simétrico \((1, 1)\) (valor 7), ambos puestos a cero con \(r=0\) (solo los puntos exactos).

🎮 Simulador EP05_02: Filtro Notch Par Simétrico

Δv 1

Δu 1

Radio (r) 0

Mueve Δv e Δu para elegir el pico — observa que el par simétrico también se filtra.

Espectro 5×5 (Rojo = Eliminado por el Filtro)

–

Figura 5.33: Simulador EP05_02: Filtro Notch

%%writefile EP05_02.cpp
// your solution

Overwriting EP05_02.cpp

TestSuite("EP05_02.cpp").run()

✔️ EP05_02.cases ya existe en casos/
📋 5 caso(s) cargado(s) de casos/EP05_02.cases

🔍 Probando C++: EP05_02.cpp
⚠️ EP05_02.cpp: archivo vacío (menos de 3 líneas). Pruebas omitidas.

5.13.3 EP05_03 🟠 Cuantización DCT: la Verdadera Fuente de Compresión

Una aplicación de galería de fotos necesita reducir el tamaño de miles de imágenes antes de hacer upload a la nube, sin recodificar todo desde cero. El ingeniero responsable ya tiene los coeficientes DCT de cada bloque \(4\times4\) calculados (la etapa costosa computacionalmente ya se ha realizado) — solo falta aplicar la tabla de cuantización, la etapa que realmente descarta información y genera compresión. Los coeficientes de alta frecuencia, menos perceptibles al ojo humano, reciben divisores grandes y tienden a convertirse en cero; los coeficientes de baja frecuencia, más perceptibles, reciben divisores pequeños y sobreviven casi intactos.

Vas a implementar exactamente esta etapa: cuantizar y descuantizar (dividir, redondear, multiplicar de vuelta) — el corazón de la compresión lossy del JPEG.

5.13.3.1 📋 Directrices de Implementación

1. Dimensión del bloque: Leer el entero \(N\) (bloque \(N \times N\)).
2. Coeficientes: Leer la matriz \(C\) de coeficientes DCT, \(N\) filas con \(N\) enteros cada una (pueden ser negativos).
3. Tabla de cuantización: Leer la matriz \(Q\), \(N\) filas con \(N\) enteros positivos cada una.
4. Cuantización: Para cada posición \((u,v)\), calcular el índice cuantizado \[ \tilde{C}(u,v) = \text{round}\!\left(\frac{C(u,v)}{Q(u,v)}\right) \] usando redondeo estándar al entero más cercano (los valores intermedios .5 nunca ocurren en los casos de prueba).
5. Descuantización (reconstrucción): Calcular \[ C'(u,v) = \tilde{C}(u,v) \times Q(u,v) \]
6. Salida: Mostrar la matriz reconstruida \(C'\), \(N \times N\), enteros.

5.13.3.2 📌 Restricciones Computacionales

• Round-trip completo: la salida es el coeficiente reconstruido (\(\tilde{C} \times Q\)), no el índice cuantizado aislado.
• División en punto flotante: la división \(C(u,v)/Q(u,v)\) debe realizarse en punto flotante antes del redondeo — la división entera truncada producirá un resultado incorrecto.
• Signo preservado: los coeficientes negativos mantienen el signo después de la cuantización y la reconstrucción.
• \(Q(u,v) > 0\) siempre: no hay necesidad de tratar la división por cero.

5.13.3.3 🧠 Fundamentación Teórica

Coeficiente	Frecuencia	Valor típico de \(Q\)	Efecto de la cuantización
\(C(0,0)\)	DC (promedio del bloque)	Pequeño	Casi siempre sobrevive — domina la energía
\(C(u,v)\) bajo \(u+v\)	Baja frecuencia	Pequeño/medio	Parcialmente preservado
\(C(u,v)\) alto \(u+v\)	Alta frecuencia	Grande	Frecuentemente se convierte en cero — fuente de la compresión

5.13.3.4 📦 Especificación de Entrada y Salida (VPL)

Entrada:

• Línea 1: Entero \(N\).
• \(N\) líneas siguientes: matriz \(C\) (coeficientes DCT, enteros, pueden ser negativos).
• \(N\) líneas siguientes: matriz \(Q\) (tabla de cuantización, enteros positivos).

Salida:

• Matriz reconstruida \(C'\), \(N \times N\), enteros separados por espacio.

5.13.3.5 📌 Ejemplos

Entrada	Salida	Observación
4 50 10 -5 0 8 -3 2 1 0 1 0 0 2 0 0 -1 2 5 7 8 4 7 8 11 6 8 11 12 9 11 12 14	50 10 -7 0 8 0 0 0 0 0 0 0 0 0 0 0	\(C(0,0)=50/2=25 \to 25\times2=50\) (preservado). \(C(0,2)=-5/7\approx-0.71\to-1\to-1\times7=-7\). Ya \(C(1,1)=-3/7\approx-0.43\to0\): anulado por la cuantización — la mayor parte del bloque se convierte en cero, ilustrando la compactación de energía en la esquina superior izquierda.

🎮 Simulador EP05_03: Cuantización DCT round(C / Q) × Q

Escala de Q (Agresividad): 1.00×

Ajuste la escala de Q y vea cuántos coeficientes sobreviven (no cero) tras el round-trip.

Coeficientes DCT (C)

Reconstruido (round(C / Q) · Q)

–

Figura 5.34: Simulador EP05_03: Cuantización DCT (round-trip)

%%writefile EP05_03.cpp
// your solution

Overwriting EP05_03.cpp

TestSuite("EP05_03.cpp").run()

✔️ EP05_03.cases ya existe en casos/
📋 5 caso(s) cargado(s) de casos/EP05_03.cases

🔍 Probando C++: EP05_03.cpp
⚠️ EP05_03.cpp: archivo vacío (menos de 3 líneas). Pruebas omitidas.

5.13.4 EP05_04 🔴 Implementando la DFT 2D a partir de la definición

Un laboratorio de investigación en astronomía computacional recibió, de una misión antigua, un pequeño sensor experimental cuyos datos brutos no pueden ser procesados por bibliotecas modernas de FFT — el entorno de validación está aislado y solo permite operaciones aritméticas básicas. El equipo necesita reimplementar la Transformada de Fourier Discreta 2D a partir de la propia definición matemática, célula por célula, para después comparar bit a bit con np.fft.fft2 en otro entorno.

Este es el ejercicio más conceptual de la lista: no hay atajos. Vas a implementar el doble sumatorio de la Ecuación 5.1 directamente, evidenciando por qué existe la FFT — y el costo computacional que evita.

5.13.4.1 📋 Directrices de implementación

1. Dimensiones: Leer los enteros \(M\) (filas) y \(N\) (columnas) de la imagen \(f(x,y)\).
2. Datos: Leer los valores enteros de \(f(x,y)\), fila a fila.
3. DFT 2D: Para cada par de frecuencias \((u,v)\) con \(u=0,\ldots,M-1\) y \(v=0,\ldots,N-1\), calcular \[ F(u,v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y)\, e^{-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)} \] usando la identidad de Euler \(e^{-j\theta} = \cos(\theta) - j\sin(\theta)\) para separar parte real e imaginaria — no se debe utilizar ninguna función de FFT predefinida.
4. Magnitud: Calcular \(|F(u,v)| = \sqrt{\text{Re}(F)^2 + \text{Im}(F)^2}\) y redondear al entero más cercano.
5. Salida: Mostrar la matriz de magnitudes redondeadas, \(M \times N\), en el mismo orden (sin fftshift — el DC permanece en \((0,0)\)).

5.13.4.2 📌 Restricciones computacionales

• Prohibido usar bibliotecas de FFT: la implementación debe calcular los sumatorios dobles explícitamente (bucles anidados), aunque sea más lenta.
• Sin fftshift: la salida mantiene la convención cruda de la DFT, con el componente DC en \(F(0,0)\) (esquina superior izquierda).
• Redondeo: la magnitud final debe redondearse al entero más cercano; en los casos de prueba no hay ambigüedad .5.
• Precisión: pequeños errores de punto flotante (del orden de \(10^{-6}\)) antes del redondeo son esperados y no afectan al resultado entero final.

5.13.4.3 🧠 Fundamentación teórica

Elemento	Significado
\(F(0,0)\)	Componente DC — suma de todos los píxeles, \(F(0,0) = \sum f(x,y)\)
Parte real \(\text{Re}(F)\)	Proyección de la señal sobre cosenos
Parte imaginaria \(\text{Im}(F)\)	Proyección de la señal sobre senos
Complejidad de esta implementación	\(\mathcal{O}((MN)^2)\) — por eso la FFT, con \(\mathcal{O}(MN\log(MN))\), resulta indispensable en imágenes reales

5.13.4.4 📦 Especificación de entrada y salida (VPL)

Entrada:

• Línea 1: Entero \(M\).
• Línea 2: Entero \(N\).
• Líneas siguientes: Elementos enteros de \(f(x,y)\), \(M\) líneas.

Salida:

• Matriz de magnitudes \(|F(u,v)|\) redondeadas, \(M \times N\), separadas por espacios.

5.13.4.5 📌 Ejemplos

Entrada	Salida	Observación
2 2 1 2 3 4	10 2 4 0	\(F(0,0)=1+2+3+4=10\) (DC = suma total). \(F(0,1)=(1-2)+(3-4)=-2 \to |F|=2\). \(F(1,0)=(1+2)-(3+4)=-4\to|F|=4\). \(F(1,1)=(1-2)-(3-4)=0\).

🎮 Simulador EP05_04: DFT 2D — Definición Directa ΣΣ f(x,y) e^-j2π(…)

Haz clic en las celdas de f(x,y) para cambiar los valores (incrementa +1; Shift + clic decrementa -1) y observa |F(u,v)| recalculado en vivo.

f(x,y) — Dominio Espacial

|F(u,v)| — Magnitud (Sin Shift)

–

Figura 5.35: Simulador EP05_04: DFT 2D manual

%%writefile EP05_04.cpp
// your solution

Overwriting EP05_04.cpp

TestSuite("EP05_04.cpp").run()

✔️ EP05_04.cases ya existe en casos/
📋 5 caso(s) cargado(s) de casos/EP05_04.cases

🔍 Probando C++: EP05_04.cpp
⚠️ EP05_04.cpp: archivo vacío (menos de 3 líneas). Pruebas omitidas.

5.13.5 EP05_05 🏆 Pipeline JPEG Completo: DCT, Cuantización y Reconstrucción

Usted ha sido contratado para crear, desde cero, un códec JPEG didáctico en un entorno embebido, sin ninguna biblioteca de imágenes disponible — solo operaciones matemáticas básicas. El cliente quiere entender exactamente dónde se pierde la calidad y dónde se recupera, bloque por bloque. Este es el desafío final del capítulo: integrar todo lo estudiado — la DCT-II ortonormal, la cuantización perceptual y la reconstrucción vía IDCT — en un único pipeline de extremo a extremo, procesando un bloque \(N \times N\) desde el inicio hasta el final, exactamente como el estándar JPEG lo hace internamente, \(8\times8\) píxeles a la vez.

5.13.5.1 📋 Directrices de Implementación

1. Dimensión del bloque: Leer el entero \(N\).
2. Bloque original: Leer la matriz de píxeles \(f(x,y)\), \(N\) líneas con \(N\) enteros en \([0,255]\).
3. Tabla de cuantización: Leer la matriz \(Q\), \(N \times N\) enteros positivos.
4. Centralización: Restar 128 de cada píxel: \(g(x,y) = f(x,y) - 128\).
5. DCT-II 2D ortonormal: Calcular \[ C(u,v) = \alpha(u)\,\alpha(v)\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} g(x,y)\,\cos\!\left[\frac{\pi(2x+1)u}{2N}\right]\cos\!\left[\frac{\pi(2y+1)v}{2N}\right] \] con \(\alpha(0)=\sqrt{1/N}\) y \(\alpha(k)=\sqrt{2/N}\) para \(k>0\).
6. Cuantización: \(\tilde{C}(u,v) = \text{round}(C(u,v)/Q(u,v))\).
7. Descuantización: \(C'(u,v) = \tilde{C}(u,v)\times Q(u,v)\).
8. IDCT-II 2D (inversa ortonormal): Calcular \(g'(x,y)\) a partir de \(C'(u,v)\) usando la transformada inversa correspondiente (misma base, sumatorio sobre \(u,v\)).
9. Reversión de la centralización y redondeo: \(f'(x,y) = \text{round}(g'(x,y) + 128)\), restringido al intervalo \([0,255]\) (clipping).
10. Salida: Mostrar el bloque reconstruido \(f'\), \(N \times N\), enteros.

5.13.5.2 📌 Restricciones Computacionales

• Pipeline completo obligatorio: todas las seis etapas (centralizar, DCT, cuantizar, descuantizar, IDCT, revertir) deben implementarse — omitir la cuantización no pasa las pruebas, pues el resultado sería idéntico al original.
• Clipping: los valores reconstruidos fuera de \([0,255]\) deben truncarse (0 si es negativo, 255 si es mayor que 255).
• Redondeo: tanto en la cuantización como en la reconstrucción final de los píxeles, use redondeo estándar; los casos de prueba evitan ambigüedad .5.
• Base ortonormal: la normalización \(\alpha(u)\) y \(\alpha(v)\) debe aplicarse exactamente como se especifica — sin ella, la IDCT no reconstruye correctamente.

5.13.5.3 🧠 Fundamentación Teórica

Etapa	Análoga en el estándar JPEG real	Dónde se pierde la calidad
Centralización	Misma — la DCT asume señal centrada en cero	Sin pérdida
DCT-II	Etapas 3–4 del pipeline (Tabla 5.7)	Sin pérdida (transformación exacta y reversible)
Cuantización	Etapa 5 — división por \(Q(u,v)\)	Principal fuente de pérdida — los coeficientes de alta frecuencia se vuelven cero
IDCT	Reconstrucción final	Reconstruye exactamente los coeficientes cuantizados, no los originales

5.13.5.4 📦 Especificación de Entrada y Salida (VPL)

Entrada:

• Línea 1: Entero \(N\).
• \(N\) líneas siguientes: bloque original \(f(x,y)\), enteros en \([0,255]\).
• \(N\) líneas siguientes: tabla de cuantización \(Q\), enteros positivos.

Salida:

• Bloque reconstruido \(f'(x,y)\), \(N \times N\), enteros en \([0,255]\), separados por espacios.

5.13.5.5 📌 Ejemplos

Entrada	Salida	Observación
4 120 130 125 128 115 140 135 122 118 150 160 130 110 120 145 138 4 6 8 10 6 8 10 12 8 10 12 16 10 12 16 20	118 126 119 131 114 143 140 119 117 149 159 130 107 121 146 139	Tras la DCT, cuantización agresiva en las altas frecuencias (valores grandes de \(Q\) en la esquina inferior derecha) y reconstrucción vía IDCT, el bloque queda cercano al original, pero no idéntico — la diferencia es el costo de la compresión lossy.

5.13.5.6 💡 Consejo de Depuración

Si el resultado no coincide, verifique en este orden: (1) los coeficientes DCT brutos (antes de la cuantización) — deben reconstruir el original exactamente vía IDCT si omite las etapas 6–7; (2) la tabla \(\alpha(u)\) — error común es aplicar \(\sqrt{2/N}\) también para \(u=0\); (3) el redondeo de la cuantización, que debe ocurrir antes de multiplicar de vuelta por \(Q\).

🎮 Simulador EP05_05: Pipeline JPEG (Bloque 4×4) DCT → Q → IDCT

Escala de Q (1 = Tabla Base, Mayor = Más Pérdida): 1.00×

Ajusta el factor de escala de cuantización y observa el bloque reconstruido alejarse (o acercarse) del original.

Bloque Original

Reconstruido (DCT → Q → IDCT)

–

Figura 5.36: Simulador EP05_05: Pipeline JPEG completo en bloques

%%writefile EP05_05.cpp
// your solution

Overwriting EP05_05.cpp

TestSuite("EP05_05.cpp").run()

✔️ EP05_05.cases ya existe en casos/
📋 5 caso(s) cargado(s) de casos/EP05_05.cases

🔍 Probando C++: EP05_05.cpp
⚠️ EP05_05.cpp: archivo vacío (menos de 3 líneas). Pruebas omitidas.