python - NumPy IFFT introducing black bars in OaA Convolution Algorithm -

i'm having trouble diagnosing , fixing error. i'm trying write oaa algorithm, described in paper.

#!/usr/bin/env python # -*- coding: utf-8 -*-  """ quick implementation of several convolution algorithms compare times  """  import numpy np import _kernel tqdm import trange, tqdm pil import image scipy.misc import imsave time import time, sleep   class convolve(object):     """ contains methods convolve 2 images """     def __init__(self, image_array, kernel, back_same_size=true):         self.array = image_array         self.kernel = kernel          # store these values accessed _lot_         self.__rangex_ = self.array.shape[0]         self.__rangey_ = self.array.shape[1]         self.__rangekx_ = self.kernel.shape[0]         self.__rangeky_ = self.kernel.shape[1]          # ensure kernel suitable convolve image         if (self.__rangekx_ >= self.__rangex_ or \             self.__rangeky_ >= self.__rangey_):             raise valueerror('must submit suitably-sized arrays')          if (back_same_size):             # pad array convolution             self.__offsetx_ = self.__rangekx_ // 2             self.__offsety_ = self.__rangeky_ // 2              self.array = np.lib.pad(self.array,      \                 [(self.__offsety_, self.__offsety_), \                  (self.__offsetx_, self.__offsetx_)],\                  mode='constant', constant_values=0)              # update these             self.__rangex_ = self.array.shape[0]             self.__rangey_ = self.array.shape[1]         else:             self.__offsetx_ = 0             self.__offsety_ = 0          # returned instead of originals         self.__arr_ = np.zeros([self.__rangex_, self.__rangey_])\      def spaceconv(self):         """ normal convolution, o(n^2*n^2). slow """          # o(n^2) part of algorithm         in trange(self.__rangex_):             j in xrange(self.__rangey_):                 # o(n^2) portion                 total = 0.0                 k in xrange(self.__rangekx_):                     t in xrange(self.__rangeky_):                         total += \                   self.kernel[k][t] * self.array[i+k][j+t]                  # update entry in self.__arr_, returned                 # http://stackoverflow.com/a/38320467/3928184                 self.__arr_[i][j] = total          return self.__arr_[self.__offsetx_\                            :self.__rangex_ - self.__offsetx_,\                            self.__offsety_\                            :self.__rangey_ - self.__offsety_]      def spaceconvdot(self):         """ same former method """          def dot(ind, jnd):             """ perform simple 'dot product' between 2              dimensional image subsets. """             total = 0.0              # o(n^2) part of algorithm             k in xrange(self.__rangekx_):                 t in xrange(self.__rangeky_):                     total += \               self.kernel[k][t] * self.array[k + ind, t + jnd]             return total          # o(n^2) part of algorithm         in trange(self.__rangex_):             j in xrange(self.__rangey_):                 self.__arr_[i][j] = dot(i, j)          return self.__arr_[self.__offsetx_\                            :self.__rangex_ - self.__offsetx_,\                            self.__offsety_\                            :self.__rangey_ - self.__offsety_]      def oaconv(self):         """ faster convolution algorithm, o(n^2*log(n)). """         numpy.fft import fft2 fft, ifft2 ifft          # solve 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_          # padding on each side, i.e. left, right, top , bottom;          # centered possible         right = diffx // 2         left = diffx - right         bottom = diffy // 2         top = diffy - bottom          # pad array         self.array = np.lib.pad(self.array,                 \                             ((left, right), (top, bottom)), \                            mode='constant', constant_values=0)          divx = self.array.shape[0] / float(self.__rangekx_)         divy = self.array.shape[1] / float(self.__rangeky_)          # let's make sure...         if not (divx % 1.0 == 0.0 or divy % 1.0 == 0.0):             raise valueerror('image not partitionable (?)')         else:             divx = int(divx)             divy = int(divy)          # list of tuples partition array         subsets = [(i*self.__rangekx_, (i + 1)*self.__rangekx_,\                     j*self.__rangeky_, (j + 1)*self.__rangeky_)\                         in xrange(divx)                  \                         j in xrange(divy)]          # padding individual blocks in subsets list         padx = self.__rangekx_ // 2         pady = self.__rangeky_ // 2          self.__arr_ = np.lib.pad(self.__arr_,              \                             ((left + padx, right + padx),  \                              (top + pady, bottom + pady)), \                            mode='constant', constant_values=0)          kernel = np.pad(self.kernel,                  \                         [(pady, pady), (padx, padx)], \                       mode='constant', constant_values=0)          # need once         trans_kernel = fft(kernel)          # transform each partition , oa on conv_image         tup in tqdm(subsets):             # slice , pad array subset             subset = self.array[tup[0]:tup[1], tup[2]:tup[3]]              subset = np.lib.pad(subset,      \                 [(pady, pady), (padx, padx)],\                 mode='constant', constant_values=0)              trans_subset = fft(subset)              # multiply 2 arrays entrywise             subset = trans_kernel * trans_subset             space = ifft(subset).real              # overlap indices , add them             self.__arr_[tup[0]:tup[1] + 2 * padx, \                         tup[2]:tup[3] + 2 * pady] += space          # crop image , back, convolved         return self.__arr_[self.__offsetx_ + padx + left   \                            :padx + left + self.__rangex_   \                               - self.__offsetx_,           \                            self.__offsety_ + pady + bottom \                            :pady + bottom + self.__rangey_ \                               - self.__offsety_]      def osconv(self):         """ convolve image using os """         numpy.fft import fft2 fft, ifft2 ifft         pass      def builtin(self):         """ convolves using scipy's convolution function - extremely         fast """         scipy.ndimage.filters import convolve         return convolve(self.array, self.kernel)   if __name__ == '__main__':     try:         import pyplot plt     except importerror:         import matplotlib.pyplot plt      image = np.array(image.open('spider.jpg'))      image = np.rot90(np.rot90(np.rot90(image.t[0])))      times = []      #for in range(3, 21, 2):     kern = _kernel.kernel()     kern = kern.kg2(11, 11, sigma=2.5, mux=0.0, muy=0.0)     kern /= np.sum(kern)        # normalize volume      conv = convolve(image, kern)     #     #    # time result of increasing kernel size     #    _start = time()     convolved = conv.oaconv()     #convolved = conv.builtin()     #    _end = time()     #    times.append(_end - _start)      #x = np.array(range(3, 21, 2))     #plt.plot(range(3, 21, 2), times)     #plt.title('kernel size vs. spaceconv time', fontsize=12)     #plt.xlabel('kernel size (px)', fontsize=12)     #plt.ylabel('time (s)', fontsize=12)     #plt.xticks(x, x)     #plt.show()      #conv = convolve(image[:2*kern.shape[0],:5*kern.shape[1]], kern)      plt.imshow(convolved, interpolation='none', cmap='gray')     plt.show()     #imsave('spider2', convolved, format='png') 

but now, when call it, black bars follows:

here's example gaussian kernel 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 problem down multiplication , ifft

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

i think has ifft of discrete gaussian kernel (for reason). it's odd because if plot

space = np.real(ifft(transformed_subset)) 

i following (no black bars, , pieces fine):

and if plot opposite, i.e.

space = np.real(ifft(transformed_kernel)) 

i again no black bars, , appears place them in right places.

what missing? i've been staring @ days, editing indices , whatnot, can't rid of tessellation !!

your problem is, kernel seems in middle:

but not - there shift (11, 11) involved, becomes obvious when @ results of convolution. center of kernel should in (0,0) , (because of modulo) kernel should this:

so changed code little (sorry never used np much, hope can gist):

 ... code ...  kernel = np.pad(self.kernel,                  \                         [(pady, pady), (padx, padx)], \                       mode='constant', constant_values=0)  #move kernel center origin:                     new_kernel=np.full_like(kernel, 0)  x,y=kernel.shape  x_2=x//2  y_2=y//2  x in xrange(x):      y in xrange(y):           n_x=(x+x_2)%x           n_y=(y+y_2)%y                   new_kernel[n_x,n_y]=kernel[x,y]     # need once             trans_kernel = fft(new_kernel)# take transform of shifted kernel   .... code ...... 

and voilà:

there no black grid!