pyFFTW中的卷积产品与scipy不同

Question

我正在尝试使用 pyfftw 来实现模仿 scipy.fftconvolve 的 FFT 卷积，以获得性能和图片作为输入：

import numpy as np
import pyfftw
a = np.ones((6000, 4000), dtype='float32')
b = np.kaiser(25, 8)
b = np.outer(b, b).astype('float32')

class fftconvolve:
    def __init__(self, A, B, domain, threads=8):
        MK =  B.shape[0]
        NK = B.shape[1]
        M = A.shape[0]
        N = A.shape[1]

        if domain =="same":
            Y = M
            X = N
        elif domain == "valid":
            Y = M - MK + 1
            X = N - NK + 1
        elif domain == "full":
            Y = M + MK - 1
            X = N + NK - 1

        self.fft_A_obj = pyfftw.builders.rfft2(A, s=(M + MK -1, N + NK -1), threads=threads)
        self.fft_B_obj = pyfftw.builders.rfft2(B, s=(M + MK -1, N + NK -1), threads=threads)
        self.ifft_obj = pyfftw.builders.irfft2(self.fft_A_obj.output_array, s=(Y, X), threads=threads)

    def __call__(self, A, B):
        return self.ifft_obj(np.fft.ifftshift(
            np.fft.fftshift(self.fft_A_obj(A)) * np.fft.fftshift(self.fft_B_obj(B))
        ))

调用它：

plan = fftconvolve(a, b, "full", threads=8)
c_1 = plan(a, b)
c_1

输出：

array([[  3.89971137e-06,   3.51986018e-05,   1.24518745e-04, ...,
          1.25271297e-04,   3.56316777e-05,   4.04627326e-06],
       [  4.91737483e-05,   2.60021159e-04,   8.61040782e-04, ...,
          8.63055116e-04,   2.61142646e-04,   4.95371969e-05],
       [  1.26523402e-04,   8.49825097e-04,   2.90915114e-03, ...,
          2.90881563e-03,   8.49568460e-04,   1.26304061e-04],
       ..., 
       [  1.28503540e-04,   8.52331228e-04,   2.91197700e-03, ...,
          2.91016186e-03,   8.51134886e-04,   1.28111642e-04],
       [  2.14206957e-05,   2.32703838e-04,   8.34190170e-04, ...,
          8.34319100e-04,   2.32750244e-04,   2.14206957e-05],
       [ -8.42595455e-06,   2.29651105e-05,   1.12404508e-04, ...,
          1.12760317e-04,   2.31778213e-05,  -8.35505125e-06]], dtype=float32)

调用 scipy 等效项：

c_2 = scipy.signal.fftconvolve(a, b, "full").astype(np.float32)
c_2

输出：

array([[  5.47012860e-06,   3.68362089e-05,   1.26135841e-04, ...,
          1.26135841e-04,   3.68362089e-05,   5.47012769e-06],
       [  3.68362089e-05,   2.48057506e-04,   8.49407224e-04, ...,
          8.49407224e-04,   2.48057506e-04,   3.68362089e-05],
       [  1.26135841e-04,   8.49407224e-04,   2.90856976e-03, ...,
          2.90856976e-03,   8.49407224e-04,   1.26135841e-04],
       ..., 
       [  1.26135841e-04,   8.49407224e-04,   2.90856976e-03, ...,
          2.90856976e-03,   8.49407224e-04,   1.26135841e-04],
       [  3.68362089e-05,   2.48057506e-04,   8.49407224e-04, ...,
          8.49407224e-04,   2.48057506e-04,   3.68362089e-05],
       [  5.47012814e-06,   3.68362089e-05,   1.26135841e-04, ...,
          1.26135841e-04,   3.68362089e-05,   5.47012814e-06]], dtype=float32)

检查输出：

c_1 == c_2

给：

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ..., 
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]], dtype=bool)

还有：

np.allclose(c_1, c_2)

给：

False

所以输出不正确。删除 fftshift 并不会改变任何事情。

在我的project 中，scipy 版本给出了正确的图像，我实现的 pyfftw 给出了模糊的输出。

编辑

我也在 double 类型 (np.float64) 中进行了测试，虽然卷积的原始结果足够接近（实际上，scipy 将卷积进行了两倍），但图片仍然很糟糕：

用 scipy 反卷积：

在这里使用自定义卷积进行反卷积：它不仅模糊，而且边缘会出现条纹：

a = np.ones((6000, 4000), dtype='float64')
b = np.kaiser(25, 8)
b = np.outer(b, b).astype('float64')

现在：

np.allclose(c_1, c_2)

返回：

True

什么可以给出这个结果？

Answer 1

这按预期工作：

class fftconvolve:
    def __init__(self, A, B, domain, threads=8):
        MK =  B.shape[0]
        NK = B.shape[1]
        M = A.shape[0]
        N = A.shape[1]

        if domain =="same":
            self.Y = M
            self.X = N
        elif domain == "valid":
            self.Y = M - MK + 1
            self.X = N - NK + 1
        elif domain == "full":
            self.Y = M + MK - 1
            self.X = N + NK - 1

        self.M = M + MK - 1
        self.N = N + NK - 1

        a = np.pad(A, ((0, MK - 1), (0, NK - 1)), mode='constant')
        b = np.pad(B, ((0, M - 1), (0, N - 1)), mode='constant')

        self.fft_A_obj = pyfftw.builders.rfft2(a, s=(self.M, self.N), threads=threads)
        self.fft_B_obj = pyfftw.builders.rfft2(b, s=(self.M, self.N), threads=threads)
        self.ifft_obj = pyfftw.builders.irfft2(self.fft_A_obj.output_array, s=(self.M, self.N), threads=threads)


        self.offset_Y = int(np.floor((self.M - self.Y)/2))
        self.offset_X = int(np.floor((self.N - self.X)/2))

    def __call__(self, A, B):
        MK =  B.shape[0]
        NK = B.shape[1]
        M = A.shape[0]
        N = A.shape[1]

        a = np.pad(A, ((0, MK - 1), (0, NK - 1)), mode='constant')
        b = np.pad(B, ((0, M - 1), (0, N - 1)), mode='constant')

        return self.ifft_obj(self.fft_A_obj(a) * self.fft_B_obj(b))[self.offset_Y:self.offset_Y + self.Y, self.offset_X:self.offset_X + self.X]