0

I am trying to write a Matlab code to produce motion and Gaussian blur. Here is my code:

f0=imread('cameraman.png');
[Nx,Ny,Nz]=size(f0);
 if Nz>1;f0=double(rgb2gray(f0));
 else
     f0=double(f0);
 end

blurfilter1 = fspecial('gaussian', [7 7], 1);
blurfilter2 = fspecial('motion',5,30);


Au = @(u) imfilter(u,blurfilter1,'symmetric');


Bu = @(u) imfilter(u,blurfilter2,'symmetric');

In above code, operators Au and Bu will produce Gaussian blur and motion blur respectively. However, I need these two operators to be self-adjoint. I have found the following code for producing self-adjoint opearots but I am not sure whether it works for both Gaussian and blur operators. Here is the code:

Au = @(u) imfilter(imfilter(u,blurfilter1,'symmetric'),blurfilter1','symmetric');
Bu = @(u) imfilter(imfilter(u,blurfilter2,'symmetric'),blurfilter2','symmetric');

I really appriciate it if anyone could please help me to build self-adjoint Au and Bu operators for Gaussian and motion blur.

eli
  • 5
  • 2

2 Answers2

1

First there are functions convmtx and convmtx2 which returns the matrix operator form of any convolution kernel.

Second, if you consider the convolution operator, the complex sinusoidal are their eigenfunctions and the amplitude of transfer function is their eigenvalues. So, if the operator is selfadjoint, it must have real eigenvalues and in order to have a selfadjoint convolution operator the Fourier transform of kernel must have real values. If a function is real valued in Fourier space then it must have even symmetry in spatial domain.

Mohammad M
  • 1,327
  • 7
  • 13
0

Gaussian

To start, fspecial('gaussian' is deprecated; the fspecial documentation explains how to do the same thing, but better.

Why do you use an n-dimensional convolution (imfilter) if you're actually after 2 sequences of 1-dimensional convolutions?

Something like (excuse me, this is untested, I uninstalled my Matlab recently)

x = [-3:3]
sigma = 1
gaussian = exp(- x.^ 2 / (2 * sigma ^ 2));
self_adjoint_filter = @(img, impresp1d) conv2(impresp1d, impresp1d, img)
Au = @(img) self_adjoint_filter(img, gaussian)

should do the trick.

Motion blur

That's not a self-adjunct operation. It has a definite direction, so eigenvalues can't be the same magnitude.

Marcus Müller
  • 30,525
  • 4
  • 34
  • 58
  • Thank you so much for your time and valuable answer. Yes, now I understand that motion blur is not a self-adjoint operator and I am working on the structure of its adjoint. Many thanks. – eli Feb 15 '23 at 06:50
  • @Marcus, Hi Marcus, maybe i'm wrong, but what if the motion blur kernel were symmetric around the origin ? – Mohammad M Feb 15 '23 at 08:32
  • then it wouldn't be a motion blur. – Marcus Müller Feb 15 '23 at 09:01
  • as an example, what is wrong with a line centered around origin to be a motion blur kernel? – Mohammad M Feb 15 '23 at 09:50
  • How's convolution with that a self-adjoint operator? The moment you transpose that, you get a complete different result. – Marcus Müller Feb 15 '23 at 11:10
  • the "convolution kernel" is not the "convolution operator", and it's the convolution operator should be self adjoint not the kernel itself. I explained in my answer below – Mohammad M Feb 15 '23 at 14:09
  • @MohammadM exactly! The operation "convolve with a specific kernel" is a linear mapping (==linear operator) and it can be self-adjoint, but that only works, exactly as you say, if its Eigenvalues are real (otherwise, the conjugate of the Eigenvalues can't be the same as the Eigenvalues). Now, take a step back: self-adjoint means: the adjoint of the operator is the same as the operator. That means, for a real-entry operator, that it's diagonally symmetrical. A convolution operator is a band-structure tensor replicating the kernel along its main diagonal,and thus,kernel symmetry is necessary. – Marcus Müller Feb 15 '23 at 14:32
  • sorry marcus, i lost you at "thus,kernel symmetry is necessary". – Mohammad M Feb 15 '23 at 15:36
  • I wrote a code to check it in practice. if the motion blur kernel (or any other kernel) were centered around origin, it will produce a self-adjoint operator. the matlab convmtx2 return the paddded output (the operator is not square), so i cant check the symmetry of the operator itself, but i used some random input to check it. – Mohammad M Feb 15 '23 at 15:38
  • kern = fspecial('motion',6,45); kern_size = size(kern); conv_op = convmtx2(kern,im_size);

    pad_size = (kern_size-1)/2; im_size = [100 100]; im_padded_size = im_size+kern_size-1;

    im = rand(im_size); im_padded = padarray(im,pad_size);

    Op_im_padded = conv_opim(:); adjointOp_im = conv_op.'im_padded(:); adjointOp_im_padded = padarray(reshape(adjointOp_im,im_size),pad_size);

    sum(Op_im_padded(:)~=Op_im_padded(:))

     – Mohammad M Feb 15 '23 at 15:38
  • @MohammadM I don't have matlab to verify your code, so I'll have to understand it by reading it, which really won't work, I'm afraid, since I don't even understand what is what in there; But: the conv_op is what's supposed to be a self-adjoint operator, right. So, as you know how to construct that from the kernel, i.e., what convmtx2 actually does, you see how the kernel structure effects the self-adjointness of the resultant operator, right? – Marcus Müller Feb 15 '23 at 15:47
  • @MohammadM also, what does being a self-adjoint operator mean, for an operator that operates on images? Right, it means that the effect of the operator is the same, whether you apply it directly, or you apply its adjoint. Now, if take the adjoint of, for example, a horizontal motion blur (the operator, not the kernel), I'd expect that to be a vertical motion blur, unless I'm really missing something here. And since these two blurs don't do the same to every image, the horizontal blur operator can't be self-adjoint. – Marcus Müller Feb 15 '23 at 15:49
  • the whole code should be 2 line, but the problem is, output of convmtx is not a square operator, it takes non-padded input but return the output as if the input were zero padded. so i have to handle these paddings myself. – Mohammad M Feb 15 '23 at 15:50
  • all i'm saying is that, if the kernel has a real valued fourier transform, then the resulting operator is self-adjoint. and to have a real valued fourier transform the kernel must have even symmetry. – Mohammad M Feb 15 '23 at 15:52
  • also about the code, it is comparing the output of operator with the adjoint operator and the outputs are the same using a random input. – Mohammad M Feb 15 '23 at 15:58
  • but then you're somehow not doing the adjoing of a non-symmetrical convolution. – Marcus Müller Feb 15 '23 at 16:05
  • transpose of a horizontal motion kernel is vertical one, but the adjoint of a horizontal motion blur operator is a motion blur where the kernel is inverted around the origin. – Mohammad M Feb 15 '23 at 16:06
  • 1
    if you consider it in time domain. the adjoint operator is like the anti-causal version of causal filter. – Mohammad M Feb 15 '23 at 16:07
  • @MohammadM I'll have to write this down as convolution tensor, but that's a good point. – Marcus Müller Feb 15 '23 at 16:11
  • @MohammadM, You may read about the adjoint in https://dsp.stackexchange.com/questions/59089. If you're after a code for the convolution matrix, see https://dsp.stackexchange.com/a/54884. – Royi Feb 17 '23 at 15:42