Gaussian blur sparse matrix with perioded boundary

Hey everyone, I am trying to create an explicit matrix for Gaussian blur without using imfilter or convolution directly. The reason is, that a want to implement a mathematical primal dual algorithm that forces me to have an explicit matrix to blur the image. The theoretical background would be to complicated here.

The idea is to reshape the given image (image -> vector) u and multiply it by a "blurring-matrix" K, such that the reshaped solution Ku (vector -> image) is a blurred image.

Thus, I have computed the Gaussian kernel with the help of fspecial. The computation should depend on the size of the kernel which is variable - and that is my problem. The explizit matrix should have the form in the attached file. I choosed A to be a 10x9 matrix and a variable kernel (3x3, 5x5 and 7x7). The boundary should be treated as periodic which yields this special form of the matrix K.

As the original image will be huge, the blurring matrix K (not the kernel!) has to be sparse.

Is there any other way (hopefully a more effective way) than using spdiags? This is my code so far:

%size of matrix to be blurred
m = 10;
n = 9;
N = m*n;
blur_para = [3,1.0]; %blur_para(1) = kernel size , blur_para(2) = standard deviation

m_blur = blur_para(1) + (1-mod(blur_para(1),2)); %size of filter should be odd for sure
supp = (m_blur-1)/2;

% create kernel for Gaussian filter with size m_blur and
% standard deviation blur_std
H = fspecial('gaussian',[m_blur m_blur],blur_para(2));
    % ! quoted from fspecial.m for 'gaussian' (with own comments):
    %     rows = sze(1); cols = sze(2);       %here: rows = cols
    %     r2 = (rows-1)/2; c2 = (cols-1)/2;   %find center points
    %     [x,y] = meshgrid(-c2:c2, -r2:r2);   %discrete points
    %     radsqrd = x.^2 + y.^2;              %evaluate Gaussian kernel...
    %     f = exp(-radsqrd/(2*sigma^2));      %... at discrete points
    %     f = f/sum(f(:));                    %normalize

%%compute operator matrix K (gaussian filter)
d = ones(m,n,m_blur*m_blur+2);

%create diagonals of K
d(:,:,1) = H(supp+1,supp+1)*d(:,:,1);

d(:,:,2) = H(supp+2,supp+1)*d(:,:,2);
d(m,:,2) = H(supp+2,supp+2);

d(:,:,3) = H(supp,supp+2)*d(:,:,3);
d(1,:,3) = H(supp,supp+1);

d(:,:,4) = H(supp+1,supp+2)*d(:,:,4);

d(:,:,5) = H(supp+2,supp+2)*d(:,:,5);
d(m,:,5) = 0;

d(:,:,6) = zeros(m,n);
d(1,:,6) = H(supp,supp+2);

d(:,:,7) = zeros(m,n);
d(m,:,7) = H(supp+2,supp);

d(:,:,8) = H(supp,supp)*d(:,:,8);
d(1,:,8) = 0;

d(:,:,9) = H(supp+1,supp)*d(:,:,9);

d(:,:,10) = H(supp+2,supp)*d(:,:,10);

d(:,:,11) = H(supp,supp)*d(:,:,11);

d = reshape(d,N,1,m_blur*m_blur+2);

K = spdiags([[d(1,:,11);zeros(N-1,1)],[d(1:m-1,:,10);zeros(N-m+1,1)],[d(1:m,:,9);zeros(N-m,1)],...
    [d(1:m+1,:,8);zeros(N-m-1,1)],[d(1:m+n,:,7);zeros(N-m-n,1)],[d(1:N-m-n,:,6);zeros(m+n,1)],...
    [d(1:N-n-2,:,5);zeros(n+2,1)],[d(1:N-n-1,:,4);zeros(n+1,1)],...
    [d(1:N-n,:,3);zeros(n,1)],[d(1:N-1,:,2);0],d(:,:,1),[0;d(1:N-1,:,2)],...
    [zeros(n,1);d(1:N-n,:,3)],[zeros(n+1,1);d(1:N-n-1,:,4)],...
    [zeros(n+2,1);d(1:N-n-2,:,5)],[zeros(m+n,1);d(1:N-m-n,:,6)],...
    [zeros(N-m-n,1);d(1:m+n,:,7)],[zeros(N-m-1,1);d(1:m+1,:,8)],...
    [zeros(N-m,1);d(1:m,:,9)],[zeros(N-m+1,1);d(1:m-1,:,10)],[zeros(N-1,1);d(1,:,11)]],...
    [-(N-1) -(N-m+1) -(N-m) -(N-m-1) -(N-2*m-1) -(2*m-1) -(m+1) -m -(m-1) -1 0 ...
    1 m-1 m m+1 2*m-1 N-2*m+1 N-m-1 N-m N-m+1 N-1],N,N);

which seems to be way to complicated since this is only for kernel size of 3x3 and not variable at all.

Anyone any idea?

Thanks.

Answer 1

[Image Analyst]

I don't know why you're doing all this complicated stuff. A circular Gaussian matrix is separable. So you can use the method of separable kernels to do two 1D convolutions/filters. Just have 3 nested for loops to pass a 1 by N vector over the image. Then do the same in the perpendicular direction, and you're done.

You say "As the original image will be huge, the blurring matrix has to be sparse." I don't see any reasoning/rationale for this. Just how big is your matrix? Why does the kernel have to be sparse????

You say "The idea is to reshape the given image (image -> vector)" Again I see no reason for this, not even having the boundary wrap around. Why do you want to do this?

Finally you say you don't want to use the highly efficient and optimized imfilter() or conv2(). Why not? Why not use the highly efficient tools that are built in instead of reinventing them yourself?

[Jasmin]

First I can't use the conv2() or imfilter() because I have to implement a special mathematical primal dual algorithm where I need to have the blurring matrix as an explicit operator. It would be way to complicated to explain all the theoretical background here. So please accept my wish to do so.

Second thing, NOT the kernel has to be sparse but the resulting blurring operator as I displayed in the .xls-file which was created with the help of a Gaussian kernel.

I am writing my diploma thesis at the moment and my prof wants me to use that circular boundary. So this is not my choice here. Huge means that I will have an (256x256)-image for example that yields a 65536x65536 blurring matrix. And this is a short example!

The aim is, as I already mentioned, to multiply this matrix with the original image to achieve the blurred image as a result.

So please, does anyone have an idea how to implement such a symmetric matrix for variable size of the kernel?

[Image Analyst]

Well sorry that I can't help you with your Ph.D. but I don't have enough information, and anyway, often graduate theses/dissertations are so specialized that it takes someone with knowledge in that very specialized field (such as your professor) to understand it, not someone like me who's not in that field and just spending a few minutes looking at a one paragraph description. For example, I don't know what a "blurring matrix" is, unless that's the filter (the point spread function), but you said that was 3 by 3 up to 7 by 7, not 65563 pixels by 65536 pixels. Even if "blurring matrix" meant the output "blurred image" I don't know how you'd get that size, you'd only get a (256+N) by (256+N) matrix - less than 100 thousand pixels, not a 4.2 gigapixel image. I have no idea why blurring an image would make an output image 62 thousand times bigger .

And I don't know what a "special mathematical primal dual algorithm" is, and I still don't know why the separable kernels won't work, like http://blogs.mathworks.com/steve/2006/10/04/separable-convolution/ or http://www.mathworks.com/matlabcentral/fileexchange/28218-separate-kernel-in-1d-kernels.

And I don't know what you mean by a variable size kernel. Is the kernel the same size for the whole image, but you can very it from run to run or image to image, or does the kernel change size as you scan one single image?

Anyway, I hope your professor can explain how to do it to you.

[Jasmin]

I am sorry but that was not my intend to force any special knowledge of you. But I don't know how to explain why I can't use convolution. I wanted to ask this question because I thought one doesn't need it as I gave examples in the attached file plus the idea to multiply such a matrix with an image. The kernel size is only set once in the beginning. But from run to run, I want this size variable to be the choice of the user as with bigger size the image will be more blurry. But I just wanted to know how I can implement matrices of this form most efficient. Maybe I should close that question, as it is to difficult to understand what I mean. Nevertheless thanks for your try.