Python: полосовой фильтр изображения

У меня есть образ данных с артефактом изображения, который выводится как синусоидальный фон, который я хочу удалить. Поскольку это одночастотная синусоидальная волна, кажется естественным преобразование Фурье и либо полосовой фильтр, либо "фильтр надреза" (где, я думаю, я использовал бы гауссовский фильтр на + -omega).

Мои данные. Красные пятна - это то, что я хочу, фоновая синусоидальная волна в kx + ky нежелательно.

Пытаясь сделать это, я замечаю две вещи:

1), просто выполнив fft и назад, я уменьшил компонент синусоидальной волны, показанный ниже. Кажется, есть какая-то фильтрация высоких частот данных, просто отправляясь туда и обратно?

import numpy as np

f = np.fft.fft2(img)                  #do the fourier transform
fshift1 = np.fft.fftshift(f)          #shift the zero to the center
f_ishift = np.fft.ifftshift(fshift1)  #inverse shift
img_back = np.fft.ifft2(f_ishift)     #inverse fourier transform
img_back = np.abs(img_back)

Это изображение img_back:

Обратное преобразование Фурье, без фильтра.

Возможно, фильтрация здесь достаточно хороша для меня, но я не уверен в этом, так как у меня нет хорошего понимания подавления фона.

2) Чтобы быть более уверенным в подавлении на нежелательных частотах, я сделал булевскую "полосу пропускания" и применил ее к данным, но преобразование Фурье игнорирует маску.

a = shape(fshift1)[0]
b = shape(fshift1)[1]

ro = 8
ri = 5
y,x = np.ogrid[-a/2:a/2, -b/2:b/2] 
m1 = x*x + y*y >= ro*ro 
m2 = x*x + y*y <= ri*ri
m3=np.dstack((m1,m2))       
maskcomb =[]
for r in m3:
    maskcomb.append([any(c) for c in r])  #probably not pythonic, sorry
newma = np.invert(maskcomb)
filtdat = ma.array(fshift1,mask=newma) 
imshow(abs(filtdat))
f_ishift = np.fft.ifftshift(filtdat) 
img_back2 = np.fft.ifft2(f_ishift) 
img_back2 = np.abs(img_back2)

Здесь результат такой же, как и раньше, потому что np.fft игнорирует маски. Исправить это было просто:

filtdat2 = filtdat.filled(filtdat.mean())

К сожалению, (но при отражении также неудивительно) результат показан здесь:

Результат фильтра полосового фильтра Brickwall.

Левый график имеет амплитуду БПФ, при этом применяется полосовой фильтр. Это темное кольцо вокруг центрального (DC) компонента. Эта фаза не отображается.

Ясно, что фильтр "brickwall" не является правильным решением. Явление создания колец из этого фильтра хорошо объяснено здесь: Что происходит, когда вы применяете фильтр кирпичной стены к 1D-набору данных.

Итак, теперь я застрял. Возможно, было бы лучше использовать один из встроенных методов scipy, но они, похоже, для 1d данных, как в этой реализации фильтра butterworth. Возможно, правильная вещь заключается в использовании fftconvolve(), как это делается здесь, чтобы размыть изображение. Мой вопрос о fftconvolve заключается в следующем: требует ли это оба изображения (изображение и фильтр) находятся в реальном пространстве? Я думаю, да, но в примере они используют гауссовский, поэтому он неоднозначный (fft (гауссовый) = гауссовый). Если это так, тогда кажется неправильным попытаться создать настоящий пространственный полосовой фильтр. Возможно, в правильной стратегии используется convolve2d() с пространственным изображением и пространственным фильтром. Если да, знаете ли вы, как сделать хороший 2d-фильтр?

Ответы

Ответ 1

Итак, одна из проблем заключается в том, что ваш синусоид имеет период, который не сильно отличается от компонентов сигнала, которые вы пытаетесь сохранить. то есть расстояние между сигнальными пиками примерно совпадает с периодом фона. Это затруднит фильтрацию.

Мой первый вопрос: действительно ли этот фон является постоянным от эксперимента к эксперименту или зависит от выборки и экспериментальной установки? Если он является постоянным, то вычитание фоновых кадров будет работать лучше, чем фильтрация.

Большинство стандартных функций фильтра scipy.signal(bessel, chebychev и т.д.), как вы говорите, предназначены для одномерных данных. Но вы можете легко расширить их до изотропной фильтрации в 2-D. Каждый фильтр в частотном пространстве является рациональной функцией f. Два представления: [a, b], которые являются коэффициентами числителя и многозначного многочлена, или [z, p, k], который является факторизованным представлением многочлена, т.е.: H(f) = k(f-z0)*(f-z1)/(f-p0)*(f-p1) Вы можете просто взять многочлен от одного алгоритмов проектирования фильтра, оценивайте его как функцию sqrt (x ^ 2 + y ^ 2) и применяйте его к данным вашей частотной области.

Можете ли вы разместить ссылку на исходные данные изображения?