Filtrado de Ruido Periódico en el Dominio de Frecuencia

1. Introducción

1.1. Ruido en imágenes: conceptos básicos

En el procesamiento digital de imágenes, el ruido es cualquier elemento no deseado que altera la fidelidad visual y dificulta el análisis posterior. Sus causas pueden ir desde fluctuaciones aleatorias en los sensores hasta interferencias generadas por el propio equipo o el entorno.

En la práctica, el ruido modifica la representación original de la escena: reduce el contraste, introduce patrones que no existen en la realidad o esconde detalles importantes. Por eso, eliminarlo o al menos reducirlo es un paso clave en áreas donde la calidad de la imagen es crítica, como el diagnóstico médico, la observación satelital o la visión por computadora.

1.2. Ruido periódico: características y causas

El ruido periódico se distingue del aleatorio porque forma patrones repetitivos y predecibles: franjas, bandas o texturas regulares que, a simple vista, parecen parte de la imagen, pero en realidad no lo son.

En el dominio espacial, estos patrones suelen mantener una orientación y amplitud constantes. En el dominio de Fourier, se hacen evidentes como picos bien definidos, situados en posiciones específicas lejos del centro de la transformada.

Entre sus causas más comunes se encuentran:

Interferencias eléctricas durante la captura, como en sistemas de escaneo o transmisión de datos.
Defectos de calibración en sensores, que generan líneas o bandas fijas.
Problemas de muestreo como el aliasing, que replican patrones indeseados.

Un ejemplo claro es el banding en imágenes satelitales causado por fallos en detectores CCD, o las resonancias magnéticas con artefactos producidos por vibraciones mecánicas.

2. Filtrado tradicional: el filtro notch

2.1. Fundamentos del filtro notch en el dominio de Fourier

El filtro notch es una herramienta diseñada para atenuar frecuencias específicas sin alterar el resto del espectro. Cuando el ruido es periódico, en la transformada de Fourier aparece como pares de picos simétricos respecto al centro. Si se eliminan o reducen esas frecuencias, gran parte del patrón indeseado desaparece, mientras que el resto de la imagen se conserva.

2.2. Implementación paso a paso

El proceso clásico consiste en:

Calcular la transformada de Fourier de la imagen para obtener su espectro.
Localizar los picos asociados al ruido —a menudo por inspección visual—.
Diseñar el filtro notch, que puede ser de corte abrupto (ideal) o con transición suave (gaussiano).
Aplicar el filtro multiplicándolo por el espectro.
Obtener la imagen filtrada mediante la transformada inversa.

def NotchFiltering(img, d0, notch_coords, n=2):
    """
    Applies a Butterworth notch filter to an image at specified frequency coordinates.

    Parameters:
        img : np.ndarray
            Input grayscale image to be filtered.
        d0 : float
            The cutoff radius of the notch filter. Frequencies within this radius will be attenuated.
        notch_coords : list of tuples
            A list of (u_k, v_k) coordinates in the frequency domain where periodic noise is present.
        n : int, optional
            The order of the Butterworth filter. Higher values result in a sharper transition. Default is 2.

    Returns:
        img_filtered : np.ndarray
            The filtered image in the spatial domain.
        magnitude_spectrum : np.ndarray
            The magnitude spectrum of the original image in the frequency domain.
        H_total : np.ndarray
            The combined notch filter applied in the frequency domain.

    Notes:
    -----
    - The function computes the Fourier Transform of the input image, applies the notch filter, 
      and then performs the inverse Fourier Transform to return the filtered image.
    - The filter is applied at all specified coordinates in [notch_coords] and their symmetric counterparts.
    """
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)

    magnitude_spectrum = 20 * np.log(np.abs(fshift) + 1)

    # Construct combined filter
    H_total = np.ones_like(img, dtype=np.float32)
    for u_k, v_k in notch_coords:
        H = ButterworthNotchFilter(img.shape, d0, u_k, v_k, n)
        H_total *= H

    # Apply the notch filter
    filtered_spectrum = fshift * H_total
    f_ishift = np.fft.ifftshift(filtered_spectrum)
    img_filtered = np.fft.ifft2(f_ishift)
    img_filtered = np.abs(img_filtered)

    return img_filtered, magnitude_spectrum, H_total

2.3. Ejemplo práctico

Supongamos que tenemos una imagen satelital con bandas horizontales. En el espectro, estas se manifiestan como picos sobre el eje vertical.

Al aplicar un filtro notch Butterworth en esas posiciones, el patrón de bandas se reduce considerablemente, conservando la mayoría de los detalles originales.

3. Limitaciones del filtro notch

hay que identificar manualmente las frecuencias a eliminar y ajustar el ancho del filtro. Esto resulta poco práctico cuando el espectro está plagado de múltiples picos dispersos, o cuando el ruido se encuentra muy cerca de la frecuencia cero (DC), ya que puede confundirse con información real de la imagen.

4. Filtrado automático de ruido cuasiperiódico

4.1. Idea general del método

El enfoque propuesto por Sur y Grédiac parte de una observación clave: cuando una imagen está contaminada con ruido periódico que se extiende por toda su superficie, ese patrón es el único elemento que se repite de manera consistente en cualquier región que seleccionemos.

Aprovechando esto, se calcula un espectro de potencia promedio a partir de múltiples parches extraídos de la imagen. Al promediar, las texturas y detalles propios de la escena tienden a cancelarse, mientras que el patrón repetitivo del ruido permanece visible en el dominio de Fourier como picos bien localizados.

4.2. Modelado y detección de picos espurios

El espectro promedio obtenido se compara con el comportamiento esperado para imágenes “naturales”, cuyo espectro de potencia sigue una ley de decaimiento en frecuencia (aproximadamente proporcional a 1/f^α).
Cualquier desviación significativa respecto a esta tendencia —en especial, picos muy por encima del nivel esperado— se interpreta como un indicio de ruido cuasiperiódico.

En vez de seleccionar manualmente las frecuencias problemáticas, el algoritmo identifica automáticamente estas “anomalías” como valores atípicos en el espectro.

4.3. Construcción del filtro notch adaptativo

Una vez localizadas las frecuencias sospechosas, se genera un mapa de picos que actúa como plantilla para construir un filtro notch. Este filtro no se limita a eliminar un par de frecuencias concretas, sino que atenúa de forma controlada todas las regiones del espectro asociadas al ruido detectado.

La imagen filtrada se reconstruye aplicando la transformada inversa de Fourier, resultando en una versión con el patrón periódico reducido o eliminado.

5. Descripción del algoritmo paso a paso

El procedimiento para eliminar ruido cuasiperiódico de manera automática se puede dividir en varias fases bien definidas, inspiradas en el método de Sur y Grédiac y adaptadas a una implementación práctica en Python:

Conversión y normalización de la imagen
La imagen se transforma a escala de grises y se normaliza a valores entre 0 y 1. Esto unifica el formato de entrada y evita que las diferencias de iluminación o codificación de color influyan en el análisis espectral.
División en parches con ventana de Hann
Se fragmenta la imagen en parches cuadrados de tamaño fijo (patch size) usando un solapamiento moderado. Cada parche se multiplica por una ventana de Hann bidimensional, lo que suaviza los bordes y reduce artefactos en el espectro causados por discontinuidades.
Cálculo del espectro promedio
Para cada parche se obtiene su espectro de potencia (módulo al cuadrado de la Transformada de Fourier). Estos espectros se combinan mediante media geométrica, que atenúa valores extremos y realza la estructura común: el patrón del ruido.
Ajuste de un modelo estadístico
Se calcula la frecuencia radial de cada punto del espectro promedio y se ajusta una ley de potencia $1/f^\alpha$ mediante regresión robusta (Huber). Este modelo describe cómo debería decaer la energía en una imagen natural sin ruido periódico.
Detección de picos anómalos
Se calculan los residuos entre el espectro real y el modelo ajustado. Aquellos puntos que superan un umbral estadístico (kσ) y que están por encima de una frecuencia mínima f₂ se marcan como outliers, es decir, posibles frecuencias de ruido.
Construcción del mapa de picos
Los outliers se organizan en una máscara binaria que respeta la simetría del espectro. Esta máscara se interpola al tamaño original de la imagen y se suaviza con un filtro gaussiano para evitar cortes bruscos.
Filtrado notch adaptativo
Se multiplica el espectro original de la imagen por el complemento de la máscara suavizada. Este paso atenúa o elimina las frecuencias asociadas al ruido detectado, dejando intactas las demás.
Reconstrucción y extracción del ruido
Mediante la transformada inversa de Fourier se obtiene la imagen filtrada. El componente de ruido se calcula restando la imagen filtrada de la original, lo que permite analizar qué ha sido eliminado.
Visualización y validación
Se muestran los espectros antes y después del filtrado, junto con el espectro del ruido, lo que permite verificar visualmente la efectividad del proceso.


def RemoveQuasiperiodicNoise(image, patch_size=128, threshold=3.0, fmax=0.61):
    """
    Removes quasiperiodic noise from images using adaptive notch filtering.
    Based on the method by Sur & Grédiac (2015) with practical adjustments.

    Parameters:
        image (np.ndarray): Input image in grayscale or color (BGR).
        patch_size (int): Size of the square patch for spectral analysis.
        threshold (float): Standard deviation factor for detecting noise peaks.
        fmax (float): Maximum frequency for noise detection.

    Returns:
        denoised_image (uint8): Filtered image (values in range 0-255).
        noise_component (uint8): Extracted noise component (values in range 0-255).
    """
    # Convert to grayscale and normalize to [0, 1]
    if image.ndim == 3:
        image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
    image = image.astype(np.float32) / 255.0
    height, width = image.shape

    # Adjust parameters according to image dimensions
    patch_size = min(patch_size, height, width)
    step = max(1, patch_size // 8)  # Overlap L/8
    f2 = 8 / patch_size             # Minimum frequency

    # Precompute Hann window
    hann_window = np.outer(np.hanning(patch_size), np.hanning(patch_size))

    # Extract patches and compute power spectra
    patches = [
        image[y:y+patch_size, x:x+patch_size] * hann_window
        for y in range(0, height - patch_size, step)
        for x in range(0, width - patch_size, step)
    ]
    power_spectra = np.array([np.abs(fftshift(fft2(p)))**2 for p in patches])

    # Average power spectrum (geometric mean)
    avg_power_spectrum = np.exp(np.mean(np.log(power_spectra + 1e-10), axis=0))

    # Radial frequencies
    fy = np.fft.fftfreq(patch_size)[:, np.newaxis]
    fx = np.fft.fftfreq(patch_size)
    f = np.sqrt(fx**2 + fy**2)
    valid_mask = (f > f2 / 4) & (f < fmax)

    # Robust fit of the power law
    log_f = np.log(f[valid_mask]).reshape(-1, 1)
    log_P = np.log(avg_power_spectrum[valid_mask]).ravel()
    model = HuberRegressor().fit(log_f, log_P)

    log_P_pred = model.predict(log_f)
    residuals = log_P - log_P_pred
    std_res = np.std(residuals)
    upper_bound = log_P_pred + threshold * std_res

    # Noise peak detection
    outliers = (log_P > upper_bound) & (f[valid_mask].ravel() >= f2)

    # Outlier map with symmetry
    outlier_mask = np.zeros_like(avg_power_spectrum, dtype=bool)
    outlier_mask[valid_mask] = outliers
    outlier_mask |= np.flip(outlier_mask, axis=0)
    outlier_mask |= np.flip(outlier_mask, axis=1)

    # Resize and smooth the mask
    outlier_map = cv.resize(outlier_mask.astype(np.float32), (width, height), interpolation=cv.INTER_LINEAR)
    outlier_map = gaussian_filter(outlier_map, sigma=2.0)

    # Protect the DC component
    cy, cx = height // 2, width // 2
    outlier_map[cy-1:cy+2, cx-1:cx+2] = 0.0

    # Notch filtering
    fft_image = fftshift(fft2(image))
    fft_filtered = fft_image * (1 - outlier_map)
    denoised_image = np.real(ifft2(ifftshift(fft_filtered)))
    noise_component = image - denoised_image

    # Normalize and convert to uint8
    denoised_image = np.clip(denoised_image * 255, 0, 255).astype(np.uint8)
    noise_component = ((noise_component - noise_component.min()) / 
                       (noise_component.max() - noise_component.min()) * 255).astype(np.uint8)

    return denoised_image, noise_component

6. Limitaciones y consideraciones prácticas

Suposición de naturalidad: el modelo espectral está pensado para imágenes naturales. En datos sintéticos o experimentales, puede no ajustarse bien.
Sensibilidad en frecuencias bajas: cuando el ruido está muy cerca de la frecuencia cero (DC), su separación respecto a la información útil se complica, pudiendo introducir pérdidas de detalle.
Ruido de alta frecuencia: no aborda el ruido aleatorio fino ni patrones no periódicos.
Confusión con detalles reales: si el patrón de ruido se parece a estructuras de la imagen (p. ej., líneas finas en la misma dirección), puede eliminar información válida o generar artefactos.

El algoritmo ofrece parámetros ajustables, como el tamaño de parche, el umbral de detección y el ancho del suavizado, que permiten optimizar su rendimiento según el tipo de imagen. Sin embargo, estos ajustes deben hacerse con cuidado, especialmente cuando se intenta eliminar frecuencias cercanas a DC, donde existe

7. Resultados

Derecha: imagen original, Centro: image resultante, Izquierda: ruido filtrado.

8. Conclusiones

El filtrado automático de ruido cuasiperiódico mediante análisis estadístico del espectro promedio ofrece una alternativa potente a los métodos manuales tradicionales. Su capacidad para identificar picos espurios sin intervención humana lo convierte en una herramienta ideal para flujos de trabajo donde se procesan grandes volúmenes de imágenes.

Si bien no es infalible —especialmente en escenas no naturales o cuando el ruido se confunde con detalles reales—, en la mayoría de los casos logra un equilibrio eficaz entre limpieza y preservación de la información visual. Además, su naturaleza adaptativa le permite enfrentarse a patrones complejos y distribuciones irregulares de ruido que serían tediosas de manejar de forma manual.

bibliografía

Frédéric Sur and Michel Grediac "Automated removal of quasiperiodic noise using frequency domain statistics," Journal of Electronic Imaging 24(1), 013003 (11 February 2015). https://doi.org/10.1117/1.JEI.24.1.013003

https://github.com/Nobody-1321/Learn-Image-Processing