2

The function to perform an N-dimensional convolution of arrays A and B in matlab is shown below:

C = convn(A,B) % returns the N-dimensional convolution of arrays A and B.

I am interested in a 3-D convolution with a Gaussian filter. If A is a 3 x 5 x 6 matrix, what do the dimensions of B have to be?

rayryeng
  • 102,964
  • 22
  • 184
  • 193
haxtar
  • 1,962
  • 3
  • 20
  • 44
  • 1
    It depends on the standard deviation of your gaussian and the desired accuracy. `B` will probably be a 3D matrix, but that is not necessary. – m7913d Jun 19 '17 at 20:52

1 Answers1

1

The dimensions of B can be anything you want. There is no set restriction in terms of size. For the Gaussian filter, it can be 1D, 2D or 3D. In 1D, what will happen is that each row gets filtered independently. In 2D, what will happen is that each slice gets filtered independently. Finally, in 3D you will be doing what is expected in 3D convolution. I am assuming you would like a full 3D convolution, not just 1D or 2D.

You may be interested in the output size of convn. If you refer to the documentation, given the two N dimensional matrices, for each dimension k of the output and if nak is the size of dimension k for the matrix A and nbk is the size of dimension k for matrix B, the size of dimension of the output matrix C or nck is such that:

nck = max([nak + nbk - 1, nak, nbk])

nak + nbk - 1 is straight from convolution theory. The final output size of a dimension is simply the sum of the two sizes in dimension k subtracted by 1. However should this value be smaller than either of nak or nbk, we need to make sure that the output size is compatible so that any of the input matrices can fit in the final output. This is why you have the final output size and bounded by both A and B.

To make this easier, you can set the size of the filter guided by the standard deviation of the distribution. I would like to refer you to my previous Stack Overflow post: By which measures should I set the size of my Gaussian filter in MATLAB?

This determines what the output size of a Gaussian filter should be given a standard deviation.

In 2D, the dimensions of the filter are N x N, such that N = ceil(6*sigma + 1) with sigma being the desired standard deviation. Therefore, you would allocate a 3D matrix of size N x N x N with N = ceil(6*sigma + 1);.

Therefore, the code you would want to use to create a 3D Gaussian filter would be something like this:

% Example input
A = rand(3, 5, 6);
sigma = 0.5; % Example

% Find size of Gaussian filter
N = ceil(6*sigma + 1);

% Define grid of centered coordinates of size N x N x N
[X, Y, Z] = meshgrid(-N/2 : N/2);

% Compute Gaussian filter - note normalization step
B = exp(-(X.^2 + Y.^2 + Z.^2) / (2.0*sigma^2));
B = B / sum(B(:));

% Convolve
C = convn(A, B);

One final note is that if the filter you provide has any of its dimensions that are beyond the size of the input matrix A, you will get a matrix using the constraints of each nck value, but then the border elements will be zeroed due to zero-padding.

rayryeng
  • 102,964
  • 22
  • 184
  • 193
  • @rayreng, thank you for your response! I have another question that came up: in creating my filter, I am passing a 3*3 Sigma matrix, to identify the sigma values for each dimension. In this case, how do I determine the overall filter size, since the sigma values are different in each case? Also, the matrix, `A` I am convolving `B` with is pretty large -- one the ranks of 900 * 500 * 1000 -- what are your thoughts? I would appreciate your feedback. – haxtar Jun 22 '17 at 19:02
  • 1
    OK, the equation is now different. You apply `6*sigma + 1` for each direction to find the size of the 3D matrix. Call these `M,N,O` then you change `meshgrid` so that it takes in two additional vectors, one for each dimension. `[X,Y,Z] = meshgrid(-M/2:M/2,-N/2:N/2,-O/2:O/2)` where `M,N,O` is found from above. After, because the sigma is not equal in each dimension, use this form of the Gaussian: https://wikimedia.org/api/rest_v1/media/math/render/svg/32a8040ca1ef2b9e0ffa5664326ac4a150f91cea, but now add an additional sum for the third dimension. After that you can convolve normally. – rayryeng Jun 22 '17 at 19:09
  • Thanks so much @rayryeng! I am running the convolution right now. It is taking a significantly long time (do you have an estimation of how long it would take to convolve a 900*500*1000-sized `A` with a 5*6*7 (this is the size my filter ended up being). I would also appreciate your feedback on a different method of creating a filter. I created a filter `F` with `mvnpdf`, a function for Multivariate Normal Distribution in matlab. In this case, F is what I would end up convolving with A. What's the significant difference between `mvnpdf` and your method of creating a gaussian. – haxtar Jun 22 '17 at 20:28
  • 1
    `mvnpdf` is good because it doesn't assume that the dimensions are independent from one another, or if the off-diagonals of the matrix are non-zero. What I specified assumes the off-diagonals are zero. However, your data is quite large and I would expect it to take a significant amount of time. You have 450M entries in that matrix. I can't give you an exact timing of how long it'll take, but expect it to take a very long time. **Maybe** you could approach this with filtering in the frequency domain but the time required to transform your input data is going to be a killer. – rayryeng Jun 22 '17 at 20:48