Eliminar un artefacto sinusoidal de un conjunto de fotogramas de película

7

Estoy haciendo un análisis post-hoc de un conjunto de datos que consiste en una serie de cuadros de película que están contaminados por un artefacto muy periódico. Me gustaría eliminar este artefacto de mis marcos.

Para facilitar el trazado, acabo de cambiar la forma de mi matriz Mde valores de píxeles [nframes, npixels], luego promedié todos los valores de píxeles para darme un vector 1D m. Así es como se ve esta señal en el dominio del tiempo. Puede ver la oscilación con bastante claridad en el recuadro ampliado.

ingrese la descripción de la imagen aquí

Luego hice un periodograma tomando Fm = rfft(m), y tracé abs(Fm)**2contra la frecuencia. Veo un pico muy fuerte a ~ 1.5Hz:

ingrese la descripción de la imagen aquí

Además de la periodicidad temporal, también parece haber un componente espacial más débil para este artefacto, ya que en el valor exacto de la frecuencia pico parece haber una variación suave en la fase a través del eje x de mis cuadros, de modo que los píxeles en el derecha tienden a retrasar píxeles a la izquierda:

ingrese la descripción de la imagen aquí

Como enfoque de fuerza bruta, he intentado simplemente filtrar cada píxel en el dominio del tiempo usando un filtro de muesca centrado en 1.5Hz. Utilicé un filtro Butterworth de orden 4 con frecuencias críticas de 1.46 y 1.52Hz (no estoy muy versado en el diseño de filtros, así que estoy seguro de que puede haber opciones más apropiadas).

Así es como se ve la señal media de píxeles después del filtrado: ingrese la descripción de la imagen aquí

Y el periodograma correspondiente: ingrese la descripción de la imagen aquí

El filtro de muesca hace un trabajo razonablemente bueno para reducir el artefacto, pero dado que básicamente parece una sinusoide estacionaria pura, no puedo evitar pensar que podría hacerlo mejor que solo atenuar esa parte del espacio de frecuencia.

Mi idea inicial (muy ingenua) fue hacer algo como:

  1. Obtenga la frecuencia, fase y amplitud de la oscilación del espectro de Fourier para cada píxel de la película.
  2. Reconstruir la oscilación en el dominio del tiempo.
  3. Restarlo de los fotogramas de la película.

Me doy cuenta de que esto no es algo que la gente suele hacer, ya que la interferencia generalmente no es tan puramente espectral y temporalmente estacionaria, pero me pregunto si podría tener sentido en mi caso.

Datos

Pila TIFF completa de 16 bits (~ 2 GB sin comprimir)

Versión de 8 bits diezmada espacialmente (~ 35 MB sin comprimir)

ali_m
fuente
A partir de una serie de fotogramas de películas, ¿puede aclarar más claramente cómo está generando el PSD exactamente?
Tarin Ziyaee
@ user4619 de manera muy cruda: para cada cuadro acabo de calcular el valor promedio de píxeles para generar un vector x, luego tomo Fx = rfft(x)y obtengo la potencia comoabs(Fx)**2
ali_m
Tienes un marco 2D y luego generas un vector 1D promedio. A lo largo de x? A lo largo de y?
Tarin Ziyaee
@ user4619 a lo largo de x e y: modifico mi película en una matriz de nframes por npixels, luego promedio en todos los píxeles
ali_m
Ok, gracias por ese detalle, es importante en el análisis. Agrega esta información a tu publicación.
Tarin Ziyaee

Respuestas:

1

Su solución propuesta (calcular una sinusoide en el dominio del tiempo basada en el pico en la FFT, luego restarla) debería funcionar, pero hay una manera más fácil de hacer esencialmente lo mismo: modificar ese valor pico en la FFT, luego tomar el inverso transformar.

Entonces, para su video rasterizado M[nframes, npixels], encuentra el contenedor de frecuencia que contiene el artefacto, luego lo aplana sistemáticamente (por ejemplo, establece su magnitud al promedio de sus vecinos) para cada píxel:

import numpy as np
nframes, npixels = np.shape(M)
# Identify the bin containing the sinusoidal artifact
# Use the average intensity for each image
m = np.mean(M, axis=1)
# Calculate the FFT
Fm = np.fft.rfft(m)
# Find the largest bin away from the low-frequency region
lowfreq = 100  # or something
badbin = lowfreq + np.argmax(Fm[lowfreq:]**2)

# Now adjust the amplitude of that bin in the FFT of each pixel
for pixel in range(npixels):
   Fpix = np.fft.rfft(M[:, pixel])
   # Scale magnitude of artifact bin to be the mean of its neighbors
   Fpix[badbin] *= np.mean(np.absolute(Fpix[[badbin-1, badbin+1]]))/np.absolute(Fpix[badbin])
   # Rewrite the time sequence of that pixel
   M[:,pixel] = np.fft.irfft(Fpix)

Esto debería funcionar si el artefacto es exactamente de amplitud y frecuencia constantes, y su frecuencia cae directamente en un submúltiplo de la longitud de la secuencia (es decir, los sinusoides representados por la FFT). En general, es posible que desee aplanar uno o dos contenedores a cada lado badbinpara lidiar con un conjunto un poco más amplio de corrupciones de banda estrecha, por ejemplo

# ...

   # Scale magnitude of artifact binS to be the mean of neighbors
   spread = 3  # flatten bins from (badbin - (spread-1)) to (badbin + (spread-1))
   # target value for new bins
   targetmag = np.mean(np.absolute(Fpix[[badbin-spread, badbin+spread]]))
   bins = range(badbin - (spread-1), badbin + spread)
   Fpix[bins] *= targetmag/np.abs(Fpix[bins])
   # ...

Si desea restringir el componente eliminado de cada píxel para que tenga la misma frecuencia y fase del artefacto detectado en la intensidad media, puede eliminar solo la proyección de la badbinmagnitud en esa fase, por ejemplo

badbinphase = np.angle(Fm[badbin])
# ...

   Ncomponent = np.abs(Fpix[badbin])*np.cos(np.angle(Fpix[badbin]) - badbinphase)
   Fpix[badbin] -= Ncomponent * np.exp(0+1j * badbinphase)
   # ...

Tenga en cuenta que el componente resultante en badbinahora siempre será desplazado de fase 90 ° (ortogonal) del global badbinphaseen cada píxel: cualquier componente de señal exactamente a esa frecuencia y fase no puede separarse del artefacto.

dpwe
fuente
¿No es eso esencialmente lo que ya estoy haciendo con el filtro de muesca? Sin embargo, sigo pensando que no es lo ideal, ya que este enfoque no tiene en cuenta la información sobre la fase de la oscilación que estoy tratando de eliminar, que es constante en el tiempo. Me parece que debería ser posible eliminar selectivamente el artefacto sin afectar la señal 'genuina' que cae en la misma banda de frecuencia.
ali_m
No es exactamente lo que hiciste: en primer lugar, es una muesca mucho más estrecha; En segundo lugar, a diferencia del filtro de Butterworth, no causa distorsión de fase. Si resta las FFT antes y después de la modificación, obtendrá un único componente distinto de cero con cierta amplitud y la fase de la espiga original. Esta es la señal que estamos eliminando en el dominio del tiempo, es decir, una sinusoide de amplitud constante exactamente a la frecuencia y fase del pico espectral. Si la señal subyacente tiene energía en esta región, aparece como ruido en la estimación, pero debería desaparecer.
dpwe
Me doy cuenta de que puede haber significado aplicar la misma fase en cada píxel, así que agregué eso a mi respuesta. Sin embargo, no es magia: no se puede separar el componente en fase de la señal y el ruido, por lo que siempre queda un 90 ° residual desplazado del artefacto estimado.
dpwe