Custom Convolution через FFT стеков изображений: проблема со смещением изображения в реальном пространстве - PullRequest
Новая
       77

Custom Convolution через FFT стеков изображений: проблема со смещением изображения в реальном пространстве

0 голосов
/ 11 апреля 2020

Мне нужен способ сделать более одной свертки более одного раза (20 или 30 раз), итеративно. Для этого я написал собственный способ свертки с помощью преобразования Фурье двух стопок изображений.

Хотя по некоторым причинам я заметил, что в результате свертки в реальном пространстве происходит сдвиг, и я не могу определить, почему и когда это происходит. Я также сопоставил свои результаты с scipy.sgnal.fftconvolve, и там также произошел сдвиг всего объекта изображения.

Кто-нибудь знает почему? И как это исправить?

код выглядит следующим образом: 

stack = np.zeros((203, 203, 203))
stack[50:75, 50:65, 60:70] = 40000
psf0 = np.ones((3, 3, 3))

После этого я использую две функции, чтобы сначала заполнить psf 

def psfpad_fromstack(self, psf):
   n = np.shape(self.stck)
   p = np.shape(psf)
   if len(n) == 2
       px = n[0] - p[0] 
       py = n[1] - p[1]
       px0 = int(px / 2)
       py0 = int(py / 2)
       pxl = int(px0 + 1)
       pyl = int(py0 + 1)
       if px % 2 == 0 and py % 2 == 0: 
           psf_out = np.pad(psf, ((px0, px0), (py0, py0)), 'constant')  
       elif px % 2 != 0 and py % 2 == 0:  
           psf_out = np.pad(psf, ((pxl, px0), (py0, py0)), 'constant') 
       elif px % 2 == 0 and py % 2 != 0: 
           psf_out = np.pad(psf, ((px0, px0), (pyl, py0)), 'constant')
       elif px % 2 != 0 and py % 2 != 0:
           psf_out = np.pad(psf, ((pxl, px0), (pyl, py0)), 'constant')
   if len(n) == 3:  
        pz = n[0] - p[0]
        px = n[1] - p[1]
        py = n[2] - p[2]
        px0 = int(px / 2)
        py0 = int(py / 2)
        pz0 = int(pz / 2)
        pxl = int(px0 + 1)
        pyl = int(py0 + 1)
        pzl = int(pz0 + 1)
        if pz % 2 == 0 and px % 2 == 0 and py % 2 == 0:  
            psf_out = np.pad(psf, ((pz0, pz0), (px0, px0), (py0, py0)), 'constant')  
        elif pz % 2 != 0 and px % 2 == 0 and py % 2 == 0:
            psf_out = np.pad(psf, ((pzl, pz0), (px0, px0), (py0, py0)), 'constant')
        elif pz % 2 == 0 and px % 2 != 0 and py % 2 == 0:
            psf_out = np.pad(psf, ((pz0, pz0), (pxl, px0), (py0, py0)), 'constant')
        elif pz % 2 != 0 and px % 2 != 0 and py % 2 == 0:
            psf_out = np.pad(psf, ((pzl, pz0), (pxl, px0), (py0, py0)), 'constant')
        elif pz % 2 == 0 and px % 2 == 0 and py % 2 != 0:
            psf_out = np.pad(psf, ((pz0, pz0), (px0, px0), (pyl, py0)), 'constant')
        elif pz % 2 != 0 and px % 2 == 0 and py % 2 != 0:
            psf_out = np.pad(psf, ((pzl, pz0), (px0, px0), (pyl, py0)), 'constant')
        elif pz % 2 == 0 and px % 2 != 0 and py % 2 != 0:
            psf_out = np.pad(psf, ((pz0, pz0), (pxl, px0), (pyl, py0)), 'constant')
        elif pz % 2 != 0 and px % 2 != 0 and py % 2 != 0:
            psf_out = np.pad(psf, ((pzl, pz0), (pxl, px0), (pyl, py0)), 'constant')
    return psf_out

и после чтобы нормализовать, из того же класса, который я использую на PSF: 

def psf_normalization(self):            
    psf = self.stck                 
    psfout = psf/np.sum(psf)          # the psf by its calculated norm
    return psfout

Наконец, я делаю эти два типа свертки: как правило, оба они дают небольшой сдвиг на полученном извилистом изображении в реальном пространстве , но в некоторых особых случаях - когда я не знаю, почему это происходит - сдвига нет (как в коде, показанном, когда размеры изображения и psf нечетны). 

stack_ft = fftn(stack)
psf_ft = fftn(psf)
conv = stack_ft * psf_ft
stack_conv = ifftshift(ifftn((conv)))
a = stack
b = psf
cnv = signal.fftconvolve(a, b, mode = 'same')
...