Как правильно рассчитать график PSD (график плотности спектра мощности) для изображений, чтобы удалить периодический шум?

Я пытаюсь удалить периодический шум из изображения с помощью PSDP, мне это удалось, но я не уверен, что то, что я делаю, на 100% правильно.
По сути, это своего рода продолжение этой видеолекции, в которой обсуждается именно эта тема на 1д сигналах.

Что я сделал до сих пор:

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

Это первая попытка:

# img link https://github.com/VladKarpushin/Periodic-noise-removing-filter/blob/master/www/images/period_input.jpg?raw=true
img = cv2.imread('./img/periodic_noisy_image2.jpg',0)
img_flattened = img.flatten()
n = img_flattened.shape[0] # 447561
fft = np.fft.fft(img_flattened, img_flattened.shape[0])

# the values range is just absurdly large, so 
# we have to use log at some point to get the
# values range to become sensible!
psd = fft*np.conj(fft)/n
freq = 1/n * np.arange(n)
L = np.arange(1,np.floor(n/2),dtype='int')

# use log so we have a sensible range!
psd_log = np.log(psd)
print(f'{psd_log.min()=} {psd_log.max()=}')
# cut off range to remove noise!
indexes = psd_log<15
# use exp to get the original vlaues for plotting comparison
psd_cleaned = np.exp(psd_log * indexes)
# get the denoised fft
fft_cleaned = fft * indexes

# in case the initial parts were affected, 
# lets restore it from fft so the final image looks well
span = 10
fft_cleaned[:span] = fft[:span]

# get back the image
denoised_img = np.fft.ifftn(fft_cleaned).real.clip(0,255).astype(np.uint8).reshape(img.shape)

plt.subplot(2,2,1), plt.imshow(img,cmap='gray'), plt.title('original image')
plt.subplot(2,2,2), plt.imshow(denoised_img, cmap='gray'), plt.title('denoise image')
plt.subplot(2,2,3), plt.plot(freq[L],psd[L]), plt.title('PSD')
plt.subplot(2,2,4), plt.plot(freq[L],psd_cleaned[L]), plt.title('PSD clean')
plt.show()

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

2. во второй попытке я просто рассчитал спектр мощности обычным способом и получил гораздо лучший результат, имхо!:

# Read the image in grayscale
img = cv2.imread('./img/periodic_noisy_image2.jpg', 0)

# Perform 2D Fourier transform
fft = np.fft.fftn(img)
fft_shift = np.fft.fftshift(fft)

# Calculate Power Spectrum Density, its the same as doing fft_shift*np.conj(fft_shift)
# note that abs(fft_shitf) calculates square root of powerspectrum, so we **2 it to get the actual power spectrum!
# but we still need to divide it by the frequency to get the plot (for visualization only)!
# this is what we do next!
# I use log to make large numbers smaller and small numbers larger so they show up properly in visualization
psd = np.log(np.abs(fft_shift)**2)

# now I can filter out the bright spots which signal noise
# take the indexes belonging to these large values
# and then use tha to set them in the actual fft to 0 to suppress them
# 20-22 image gets too smoothed out, and >24, its still visibly noisy
ind = psd<23
psd2 = psd*ind
fft_shift2 = ind * fft_shift
# since this is not accurate, we may very well endup destroying 
# the center of the fft which contains low freq important image information
# (it has large values there as well) so we grab that area from fft and copy
# it back to restore the lost values this way!
cx,cy = img.shape[0]//2, img.shape[1]//2
area = 20
# restore the center in case it was overwritten!
fft_shift2[cx-area:cx+area,cy-area:cy+area] = fft_shift[cx-area:cx+area,cy-area:cy+area]

ifft_shift2 = np.fft.ifftshift(fft_shift2)
denoised_img = np.fft.ifftn(ifft_shift2).real.clip(0,255).astype(np.uint8)

# Get frequencies for each dimension
freq_x = np.fft.fftfreq(img.shape[0])
freq_y = np.fft.fftfreq(img.shape[1])

# Create a meshgrid of frequencies
freq_x, freq_y = np.meshgrid(freq_x, freq_y)

# Plot the PSD
plt.figure(figsize=(10, 7))
plt.subplot(2,2,1), plt.imshow(img, cmap='gray'), plt.title('img')
plt.subplot(2,2,2), plt.imshow(denoised_img, cmap='gray'), plt.title('denoised image')
#plt.subplot(2,2,3), plt.imshow(((1-ind)*255)), plt.title('mask-inv')
plt.subplot(2,2,3), plt.imshow(psd2, extent=(np.min(freq_x), np.max(freq_x), np.min(freq_y), np.max(freq_y))), plt.title('Power Spectrum Density[cleaned]')
plt.subplot(2,2,4), plt.imshow(psd, extent=(np.min(freq_x), np.max(freq_x), np.min(freq_y), np.max(freq_y))),plt.title('Power Spectrum Density[default]')
plt.xlabel('Frequency (X)')
plt.ylabel('Frequency (Y)')
plt.colorbar()
plt.show()

Кажется, это работает, но я не получаю хорошего результата, я не уверен, делаю ли я здесь что-то не так, или это лучшее, чего можно достичь!

3. Далее я попытался полностью обвести прямоугольником все яркие пятна и установить для них все нули. Таким образом я гарантирую, что окружающие значения также будут учтены в максимально возможной степени, и это то, что я получаю выход:

img = cv2.imread('./img/periodic_noisy_image2.jpg')
while (True):
    # calculate the dft
    ffts = np.fft.fftn(img)
    # now shift to center for better interpretation
    ffts_shifted = np.fft.fftshift(ffts) 
    # power spectrum
    ffts_shifted_mag = (20*np.log(np.abs(ffts_shifted))).astype(np.uint8)
    # use selectROI to select the spots we want to set to 0!
    noise_rois = cv2.selectROIs('select periodic noise spots(press Spc to take selection, press esc to end selection)', ffts_shifted_mag,False, False,False)
    print(f'{noise_rois=}')
    # now set the area in fft_shifted to zero 
    for y,x,h,w in noise_rois:
        # we need to provide a complex number!
        ffts_shifted[x:x+w,y:y+h] = 0+0j

    # shift back
    iffts_shifted = np.fft.ifftshift(ffts_shifted)
    iffts = np.fft.ifftn(iffts_shifted)

    # getback the image
    img_denoised = iffts.real.clip(0,255).astype(np.uint8)

    # lets calculate the new image magnitude
    denoise_ffts = np.fft.fftn(img_denoised)
    denoise_ffts_shifted = np.fft.fftshift(denoise_ffts)
    denoise_mag = (20*np.log(np.abs(denoise_ffts_shifted))).astype(np.uint8)

    cv2.imshow('img-with-periodic-noise', img)
    cv2.imshow('ffts_shifted_mag', ffts_shifted_mag)
    cv2.imshow('denoise_mag',denoise_mag)
    cv2.imshow('img_denoised', img_denoised)
    # note we are using 0 so it only goes next when we press it, otherwise we cant see the result!
    key = cv2.waitKey(0)&0xFF
    cv2.destroyAllWindows()

    if key == ord('q'):
        break

Опять же, я предположил, что, удалив эти периодические шумы, изображение будет выглядеть намного лучше, но я все равно вижу узоры, а это значит, что они не удалены полностью! но в то же время я убрал яркие пятна!

Еще сложнее (пока невозможно) обесшумить это изображение с помощью этого метода:

Это явно периодический шум, так что же я здесь упускаю или делаю неправильно?

Для справки, это другое изображение с периодическим шумом, с которым я экспериментировал:

Обновлять :

Прочитав комментарии и предложения, я пришел к следующей версии, которая в целом работает прилично, но я столкнулся со следующими проблемами:

1. Я не получаю крошечных мнимых значений, даже если результат выглядит довольно хорошо! Кажется, я не могу полагаться на эту проверку, чтобы увидеть, что пошло не так, она существует, когда шума очень мало/едва заметно, и когда шум повсюду!
2. На некоторых изображениях по-прежнему наблюдается значительное количество шума (приведен пример). Было бы здорово узнать, ожидается ли это, и мне следует двигаться дальше, или есть что-то не так, и нужно устранить это.

def onchange(x):pass
cv2.namedWindow('options')
cv2.createTrackbar('threshold', 'options', 130, 255, onchange)
cv2.createTrackbar('area', 'options', 40, max(*img.shape[:2]), onchange)
cv2.createTrackbar('pad', 'options', 0, max(*img.shape[:2]), onchange)
cv2.createTrackbar('normalize_output', 'options', 0, 1, onchange)

while(True):

    threshold = cv2.getTrackbarPos('threshold', 'options')
    area = cv2.getTrackbarPos('area', 'options')
    pad = cv2.getTrackbarPos('pad', 'options')
    normalize_output = cv2.getTrackbarPos('normalize_output', 'options')
    
    input_img = cv2.copyMakeBorder(img, pad, pad, pad, pad, cv2.BORDER_REFLECT) if pad>0 else img
    
    fft = np.fft.fftn(input_img)
    fft_shift = np.fft.fftshift(fft)
    # note since we plan on normalizing the magnitude spectrum,
    # we dont clip and we dont cast here!
    # +1 so for the images that have 0s we dont get -inf down the road and dont face issues when we want to normalize and create a mask out of it!
    fft_shift_mag = 20*np.log(np.abs(fft_shift)+1)
    # now lets normalize and get a mask out of it, 
    # the idea is to identify bright spot and set them to 0
    # while retaining the center of the fft as it has a lot
    # of image information 
    fft_shift_mag_norm = cv2.normalize(fft_shift_mag, None, 0,255, cv2.NORM_MINMAX)
    # now lets threshold and get our mask
    if img.ndim>2:
        mask = np.array([cv2.threshold(fft_shift_mag_norm[...,i], threshold, 255, cv2.THRESH_BINARY)[1] for i in range(3)])
        # the mask/img needs to be contiguous, (a simple .copy() would work as well!)
        mask = np.ascontiguousarray(mask.transpose((1,2,0)))
    else:
        ret, mask = cv2.threshold(fft_shift_mag_norm, threshold, 255, cv2.THRESH_BINARY)
        
    w,h = input_img.shape[:2]
    cx,cy = w//2, h//2
    mask = cv2.circle(mask, (cy,cx), radius=area, color=0, thickness=cv2.FILLED)

    # now that we have our mask prepared, we can simply use it with the actual fft to 
    # set all these bright places to 0
    fft_shift[mask!=0] = 0+0j
    ifft_shift = np.fft.ifftshift(fft_shift)
    img_denoised = np.fft.ifftn(ifft_shift).real.clip(0,255).astype(np.uint8)
    img_denoised = img_denoised[pad:w-pad,pad:h-pad]
    
    # check the ifft imaginary parts are close to zero otherwise sth is wrong!
    almost_zero = np.all(np.isclose(ifft_shift.imag,0,atol=1e-8))
    if not almost_zero:
        print('imaginary components not close to 0, something is wrong!')
    else:
        print(f'all is good!')
        
    # do a final contrast stretching: 
    if normalize_output:
        p2, p98 = np.percentile(img_denoised, (2, 98))
        img_denoised = img_denoised.clip(p2, p98)
        img_denoised = cv2.normalize(img_denoised, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)
        
    cv2.imshow('input_img', input_img)
    cv2.imshow('fft-shift-mag-norm', fft_shift_mag_norm)
    cv2.imshow('fft_shift_mag', ((fft_shift_mag.real/fft_shift_mag.real.max())*255).clip(0,255).astype(np.uint8))
    cv2.imshow('mask', mask)
    cv2.imshow('denoised', img_denoised)

    key = cv2.waitKey(30)&0xFF
    if key == ord('q') or key == 27:
        cv2.destroyAllWindows()
        break

относительно хороший результат:

Не так много! Это единственный пример, когда я до сих пор слышу много шума. Я не уверен, что это лучшее, на что я могу рассчитывать, или еще есть куда улучшаться:

Есть и другие образцы, где мне тоже не удалось убрать весь шум, например этот (я мог бы его немного подправить, но артефакты всё равно остались бы):

Я объяснил это низким качеством самого изображения и принял его. Однако я ожидал, что во втором примере есть место для улучшений, и думал, что в конечном итоге смогу получить что-то вроде этого или близко к этому:

• Мои предположения неверны?
• Являются ли эти артефакты/шумы, которые мы видим на выходе, периодическим шумом или каким-то другим типом шума?
• Условно говоря, Это лучшее, чего можно достичь/на что надеяться при использовании этого метода? Я имею в виду просто удаление периодического шума и не прибегая к чему-то продвинутому?