# How to exclude diffuse moving features in a stack of images

I have some stack of grayscale reconstructed images represented as a 3-dimensional matrix DPC with dimensions [Nz,Nx,Nt]. A single page/frame/image would be an NxMx1 slice of DPC. Each frame has some diffuse artifacts spread out over the image swamping out the signal of interest. If I were to animate the data where I am displaying a new frame every 0.1 seconds, the signal of interest is faintly visible as a static, non-moving feature of varying intensity while the diffuse artifact spreads outward and upwards

I am looking to filter out all the "moving" diffuse artifacts so I'm just left with the signal. If I take an average over all the frames, I get something that looks like this:

While the "static" components have a much higher contrast when I take the average, I suspect the "moving" components of the image are still degrading the image quality.

I processed this data in matlab and c++ but this is more of a general image processing question.

Clarification

@ChristophRackwitz without going into overwhelming detail, each frame is a phase difference image. The moving noise (or artifacts) in each frame corresponds to an increasingly defocused copy of the fully in focus image. The algorithm that reconstructs any given frame of the image produces a dataset where the defocused artifact is mixed in to the actual signal, so fixing it before it gets corrupted is non-trivial.

@AnderBiguri in the comments is likely correct that this might be more of a research problem than a forum question. However, I was hoping to double check if I had missed any common image processing algorithms that might be applicable in this situation. (aside: since I was thinking of image processing algorithms, I thought SO was the better site, but signal processing works too)

Essentially, I wanted to view this problem as some sort of image processing problem where I'm trying to filter moving objects and extract the static objects. However, a better description might be a moving, hazy background with static objects.

For example, if you had a timelapse camera shot of some busy public square and you wanted to remove all the people from the image, you could apply median filtering. Median filtering is used to filter out moving objects in a series of images to extract a static background. Is there a comparatively simple algorithm to extract static objects from a moving background?

• the animation and temporal sum/mean images greatly help understand this! good visualization is rare on this site. -- tough problem. garbage data, noise drowns out signal. can you fix this signal before it gets corrupted this badly? what can you tell us about how this signal was formed? Jun 6 at 19:46
• This looks like one of those problems where you need to search the scientific literature for answers, not a forum on how to solve programming bugs, unforutnately. Jun 7 at 10:39

So I tried Gillespie's suggestion of TV.

First, I got the individual frames from your GIF out using ffmpeg:

ffmpeg -i Q83349.gif -vsync 0 temp%d.png

Then I tried Gaussian smoothing all the individual frames and finding the average of the result.

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import scipy.ndimage as ndimage

mean_picture = np.zeros([338,428,4])

for idx in np.linspace(1,51):
file_name = 'temp%d.png' % idx
img = ndimage.gaussian_filter(mpimg.imread(file_name), sigma=(1, 1, 0), order=0)
mean_picture = mean_picture + img/51

plt.figure(1)
plt.imshow(noisy)
plt.title('Original')
plt.figure(2)
plt.imshow(img)
plt.title('Last image processed')
plt.figure(3)
plt.imshow(mean_picture)
plt.title('Overall mean of Gaussian smoothed')


Then I tried the Bregman TV approach on each frame before also finding the mean. To do this I replaced the gaussian_filter line above with

from skimage.restoration import denoise_tv_chambolle, denoise_tv_bregman
img = denoise_tv_bregman(noisy, weight=2)


This didn't seem to improve things much:

Then I wondered what a three-image median might do across the image train.

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

median_picture = np.zeros([338,428,4])

for idx in np.linspace(2,50):
file_name_1 = 'temp%d.png' % (idx-1)
file_name_2 = 'temp%d.png' % idx
file_name_3 = 'temp%d.png' % (idx+1)
img = np.median([noisy_1, noisy_2, noisy_3], axis=0)
median_picture = median_picture + img/49

plt.imshow(median_picture)
plt.title('Overall mean of medians')


This seemed to improve on the others slightly, as there is some finer detail visible.

I don't think these are as good$$^1$$ as your image. If you can share how you arrived at that, it would probably help.

$$^1$$ For certain values of good.

• 1) You can't just naively apply TVD. It requires tuning. 2) Try applying TVD between images, rather than on individual ones and then averging them. Jun 7 at 19:59
• @Gillespie I tried different smoothing parameters (weight) but it didn't seem to help. I'll see about doing it across images, though I'm not sure how to do that! :-)
– Peter K.
Jun 7 at 20:28
• Thanks for the suggestions and experimentation! I did take the median of my entire image stack, but it didn't occur to me to take a sliding median filter. That works much better. Regarding the total variation approaches: I suspect something along these lines will be the most useful. Since there is a "flowing" of the background upwards and outwards with time (or in other words, since the background has some dependence on the time axis of the animation) I suspect incorporating this into the regularization norm might yield good results. Jun 7 at 20:30
• I can't share much more information about the problem without getting too deep into the physics/science behind my experiment, and I suspect that would be far outside the scope of any reasonable stack exchange question. Jun 7 at 20:33
• @drakon101 OK! Thanks for the info. I'll play a bit more with TV to see if I can get it to work across images. I'll update if I manage to get it.
– Peter K.
Jun 7 at 20:41

Assuming there is no underlying data-specific mechanism that can be exploited, I think your best bet would be Total Variation Denoising. It reduces noise by minimizing the variance of the data, subject to the solution fitting the data. This has advantages over methods such as averaging in that it preserves edges, which is especially important for images.

Searching for packages that implement it, I came across one in MATLAB, though I'm sure there are others you could find.

It seems that the «artifacts» roll vertically at constant speed? I wonder if that could be exploited somehow. I.e. a signal model consisting of desired signal + uncorrelated noise + «slow innovation, constant velocity noise»?