Блог

Я выполняю некоторую обработку изображений на периодических изображениях в python. Мне нужно проанализировать реальное изображение с довольно четко определенной периодичностью. Долгосрочная цель этого проекта состоит в том, чтобы мой код обнаружил незначительные искажения в изображении и исправил эти искажения. Для этого мне нужна информация о фазе для периодичности на изображении.

Мне нужно сделать это, умножив исходное реальное изображение на сложную экспоненту, колеблющуюся вокруг пика в БПФ изображения, соответствующего одной определенной периодичности. Затем примените фильтрацию нижних частот с соответствующим значением для сигмы, где сигма-радиус гауссовой маски, используемой для сглаживания. Абсолютное значение отфильтрованного изображения должно дать мне амплитуду этой конкретной периодичности, а np.угол(фильтрованное изображение) или np.arctan2(мнимое/реальное) должны дать мне фазу, соответствующую этой конкретной периодичности.

Мой код отлично работает для определения амплитуд, но не для фазы. Я застрял на этом уже несколько недель и не могу понять, почему это работает неправильно.

Ниже приведен мой самый обновленный код для конкретного примера. Zraw в данном случае представлял собой всего лишь (256 x 256) 2D-изображение плоских волн. На этом изображении Zraw, который я использую для этого примера, нанесен рядом с его FFT. Я выбрал v1 на основе самого правого пика БПФ.

 from scipy.ndimage import gaussian_filter


v1 = np.array([145,128])   # Pixel position of a peak in Zraw's FFT
x0, y0 = 151, 130          # Pixel position of a maxima in the original image
sigma = 10                 # Radius in pixels of gaussian mask

N = len(Zraw[0])           # Get the number of pixels of the image (NxN)
ctr = np.floor(N/2)   1    # Get the center pixel position


# Reciprocal lattice vectors
q1 = ((v1 - np.array([ctr,ctr])) * 2* np.pi) / N

# Create meshgrid 
x = np.linspace(0,N,N)
[X,Y] = np.meshgrid(x,x)

# Create complex image --  multiply real Zraw image by complex exponential
Z1 = Zraw * np.exp(1j*(q1[0]*(X-x0)   q1[1]*(Y-y0)))

# Low pass gaussian filtering
# Filter real and imaginary parts separately because complex inputs not allowed in scipy's gaussian_filter function
Z1_real = gaussian_filter(np.real(Z1), sigma)
Z1_imag = gaussian_filter(np.imag(Z1), sigma)

# Combine filtered real amp; imag parts to get back a complex image
Z1_filt = Z1_real   1j*Z1_imag


# Amplitudes
Z1_amplitudes = np.abs(Z1_filt)

# Phase
Z1_phase = np.angle(Z1_filt)
 

I really need help with figuring out why my code isn’t working correctly — it’s not extracting the correct phase. Here’s an example of my code’s output using sigma=10. I know that the images I get depend a lot on the value I choose for sigma, the smoothing radius, but I find that regardless of the sigma I use, I don’t get a correct phase image.

I’ve tried many different techniques I’ve found online, I just can’t see where exactly I’m going wrong. Previously, I also tried doing the low-pass filtering directly in Fourier space:

     fftZ1 = npf.fftshift(npf.fft2(Z1))
    mask = np.zeros((N,N))
    mask  = np.exp(-((X-ctr)**2   (Y-ctr)**2)/(2*sigma**2))
    
    fftZ1_filt = mask * fftZ1
    
    Z1_filt = npf.ifft2(npf.ifftshift(fftZ1_filt))

    Z1_amplitudes = np.abs(Z1_filt)
    Z1_phase = np.angle(Z1_filt)
 

which gave me very similar results for the phase, but I found that my other approach worked better for retrieving the correct values for the amplitudes, so I’ve been sticking with that.

Я не могу понять, почему фаза выглядит так, я просмотрел весь Google и stack overflow, пробуя то, что другие предлагают для подобных проблем, но я не могу понять это для своего случая. Извините за очень длинный пост, я хотел все прояснить. Я довольно новичок в обработке изображений и Фурье, поэтому буду признателен за любую помощь, большое вам спасибо!