Dealing with the Cyclic Boundary Conditions of Frequency Domain Convolution in Order to Apply Linear Convolution

Question

def fft_convolution(image, kernel):
    imageC = image.copy()
    kSize = (kernel.size // 2)+1
    kernelShape = tuple(ti//2 for ti in kernel.shape)
    centery ,centerx = kernelShape
    imagePad = np.pad(image, kSize, mode = "edge")
    imageC = np.fft.fft2(imagePad)
    kernel= np.fft.fft2(kernel, imagePad.shape)
    output = np.real(np.fft.ifft2(np.multiply(imageC, kernel)))
    return output[kSize+centery:-kSize+centery, kSize+centerx:-kSize+centerx]

So my code currently takes an image M x N and kernel P x K Pad on all sides the kernelSize // 2 + 1 and then make the kernel return from fft2 as the correct shape

I then crop using the shape of the kernel from my pad size, I dont really know why, something to do with circular convolutions maybe.

But this is really confusing and makes no sense to me especially when we get to cropping where I had to deal with artifacts of cylic convolution.

I am wondering what is the correct method to apply cropping and padding perfectly, taking into account i want my image to stay the same size and the edges must be interpolated rather than using constant zero values.

Answer 1

I created a MATLAB function which is basically conv2() applied in Frequency Domain:

function [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% Applies Image Convolution in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double' (Single Channel).
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           Sets the convolution shape.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix (Single Channel).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv2()' - https://www.mathworks.com/help/matlab/ref/conv2.html.
% Remarks:
%   1.  A
% TODO:
%   1.  
%   Release Notes:
%   -   1.0.000     29/04/2021  Royi Avital     [email protected]
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

switch(convShape)
    case(CONV_SHAPE_FULL)
        numRowsFft  = numRows + numRowsKernel - 1;
        numColsFft  = numCols + numColsKernel - 1;
        firstRowIdx = 1;
        firstColIdx = 1;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
    case(CONV_SHAPE_SAME)
        numRowsFft  = numRows + numRowsKernel;
        numColsFft  = numCols + numColsKernel;
        firstRowIdx = ceil((numRowsKernel + 1) / 2);
        firstColIdx = ceil((numColsKernel + 1) / 2);
        lastRowIdx  = firstRowIdx + numRows - 1;
        lastColdIdx = firstColIdx + numCols - 1;
    case(CONV_SHAPE_VALID)
        numRowsFft = numRows;
        numColsFft = numCols;
        firstRowIdx = numRowsKernel;
        firstColIdx = numColsKernel;
        % The Kernel when transformed is shifted (Namely its (0, 0) is top
        % left not middle).
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
end

mO = ifft2(fft2(mI, numRowsFft, numColsFft) .* fft2(mH, numRowsFft, numColsFft), 'symmetric');
mO = mO(firstRowIdx:lastRowIdx, firstColIdx:lastColdIdx);


end

It is fully compatible and validated.
The full code is available on my StackExchange Signal Processing Q74803 GitHub Repository (Look at the SignalProcessing\Q74803 folder).

If you want to apply image filter with Constant / Mirror or Replicate boundary mode just pad the image before the function and use CONV_SHAPE_VALID as convShape. Without padding the boundary condition applied is Periodic.