lundi 13 juin 2016

NumPy IFFT introducing black bars in OaA Convolution Algorithm

I'm having trouble diagnosing and fixing this error. I'm trying to write the OaA algorithm, described in this paper.

import numpy as np
from numpy.fft import fft2 as FFT, ifft2 as IFFT

class convolve(object):

    """ contains methods to convolve two images with different convolution 
    algorithms """

    def __init__(self, image_array, kernel):
        self.array = image_array
        self.kernel = kernel

        self.__rangeX_ = self.array.shape[0]
        self.__rangeY_ = self.array.shape[1]

        self.__rangeKX_ = self.kernel.shape[0]
        self.__rangeKY_ = self.kernel.shape[1]

   def OaAconvolve(self):
        """ faster convolution algorithm, O(N^2*log(n)). """

        # solve for the total padding along each axis
        diffX = (self.__rangeKX_ -                      
                (self.__rangeX_ - self.__rangeKX_*( 
                    self.__rangeX_ // self.__rangeKX_)))
                    % self.__rangeKX_

        diffY = (self.__rangeKY_ -                      
                    (self.__rangeY_ - self.__rangeKY_*( 
                    self.__rangeY_ // self.__rangeKY_)))
                    % self.__rangeKY_

        # each side, i.e. left, right, top and bottom
        right = diffX // 2
        left = diffX - right
        bottom = diffY // 2
        top = diffY - bottom

        # pad the array [(top, bottom), (left, right)]
        self.array = np.pad(self.array, [(left, right), (top, bottom)], 
                mode='constant', constant_values=0)

        # return a list of tuples to partition the array
        subsets = [(i*self.__rangeKX_, (i + 1)*self.__rangeKX_,           
                 j*self.__rangeKY_, (j + 1)*self.__rangeKY_)              
                 for i in xrange(self.array.shape[0] // self.__rangeKX_)  
                 for j in xrange(self.array.shape[1] // self.__rangeKY_)]

        # set up OaA method
        padX = self.__rangeKX_ // 2
        padY = self.__rangeKY_ // 2

        self.kernel = np.pad(self.kernel, [(padY, padY), (padX, padX)], 
                    mode='constant', constant_values=0)

        transformed_kernel = FFT(self.kernel)

        # create a blank array of the same size
        convolved_image = np.zeros([self.array.shape[0] + 2*padX, 
                                    self.array.shape[1] + 2*padY])

        print self.array.shape

        # transform each partition and OaA on the convolved_image
        for tup in tqdm(subsets):
            # slice and pad the array subset
            subset = np.pad(self.array[tup[0]:tup[1], tup[2]:tup[3]], 
                            [(padY, padY), (padX, padX)], 
                            mode='constant', constant_values=0)

            # transform
            transformed_subset = FFT(subset)

            # multiply the two arrays entrywise and take the IFFT. np.real()
            # is used because some residual/negligible imaginary terms are 
            # left over after the IFFT.
            space = np.real(IFFT(transformed_kernel*transformed_subset))

            # overlap with indices and add them together to build the image
            convolved_image[tup[0]:tup[1] + 2*padX,tup[2]:tup[3] + 2*padY] += space

        # crop image and get it back, convolved
        return convolved_image[padX + left:padX + left + self.__rangeX_, 
                               padY + bottom:padY + bottom + self.__rangeY_]

But now, when I call it, I get black bars as follows:

enter image description here

Here's an example Gaussian kernel I am using.

[[ 0.          0.02390753  0.03476507  0.02390753  0.        ]
 [ 0.02390753  0.06241541  0.07990366  0.06241541  0.02390753]
 [ 0.03476507  0.07990366  0.10040324  0.07990366  0.03476507]
 [ 0.02390753  0.06241541  0.07990366  0.06241541  0.02390753]
 [ 0.          0.02390753  0.03476507  0.02390753  0.        ]]

I believe I've narrowed the problem down to the multiplication and IFFT

space = np.real(IFFT(transformed_kernel*transformed_subset))

I think it has something to do with the IFFT of discrete Gaussian kernel (for some reason). It's odd because if I just plot

space = np.real(IFFT(transformed_subset))

I get the following (no black bars, and it pieces it back together fine):

okay

And if I plot the opposite, i.e.

space = np.real(IFFT(transformed_kernel))

I again get no black bars, and it appears to place them in the right places.

enter image description here

What am I missing? I've been staring at this for days, editing indices and whatnot, but I can't get rid of this tessellation !!

Aucun commentaire:

Enregistrer un commentaire