Thread Subject: Image Convolution Equations

Hi Everyone, I am trying to figure out how to do L = H*B, where L is an image I have, and B will be some created filter. I was told I had to write the convolution equation for each pixel, and then solve Ax = b for A (in this case A corresponds to all H's, X to all B's and b to all L's). However, I only have access to L and B, so I need to get H. I don't know how to loop through the matrix L and create a convolution equation using H and B, when I do not have H.

If anyone has any suggestions, please message me.

Thank you,

-Bc

"Brenda " <bcortez23@gmail.com> wrote in message <hkojdf$d8k$1@fred.mathworks.com>...
> Hi Everyone, I am trying to figure out how to do L = H*B, where L is an image I have, and B will be some created filter. I was told I had to write the convolution equation for each pixel, and then solve Ax = b for A (in this case A corresponds to all H's, X to all B's and b to all L's). However, I only have access to L and B, so I need to get H. I don't know how to loop through the matrix L and create a convolution equation using H and B, when I do not have H.
================

The equation L=H*B can be solved by doing H=L/B

However, if B is a known filter, it is computationally cheaper to solve the equation in the Fourier domain. For example, if you have no noise in L, you would do this

 fft2(H)=fft2(L)./fft2(B)

Since you probably do have noise in L, you more likely want to do this

 fft2(H)=fft2(L)./(fft2(B) +beta*fft(C))

where C is a high-pass filter and beta is a regularizing weight.

Thank you for your comment. However, I would like to elaborate on this problem. I am supposed to write an equation for each pixel as follows:

L(m,n) = sum(i=1-inf) sum(j=1-imf) H(i,j)*B(m-i,n-j) ------ (* means multiplication now)

and for L(1,1), where I would get B to be a 3x3 filter:

L(1,1) = H[0,0]*B[1,1] + H[1,0]*B[0,1] + H[2,0]B[-1,1]...

So here, we take into account having to sum all the elements around the pixel of interest L(1,). I am supposed to construct this for all pixels of this image - 96x96.

Then having L = H*B equations, I think of these as Ax = b and then use what you suggested to get the values of A, which in this case is my H. I know what the L values are and can opt for a guess of B, so how do I create these equation in matlab? I cannot come up with 96x96 equations by hand.

Thank you all...

This is what the convolution function, conv2(), does. Why don't you
just use it? Or are you supposed to basically rewrite it on your
own? If you have to do that, there are intelligent ways to optimize
it that are faster than the brute force, but intuitive, dumb way of
doing it. (I'm still talking about the spatial domain here.)

"Brenda " <bcortez23@gmail.com> wrote in message <hkpia9$o0t$1@fred.mathworks.com>...
> Thank you for your comment. However, I would like to elaborate on this problem. I am supposed to write an equation for each pixel as follows:
>
> L(m,n) = sum(i=1-inf) sum(j=1-imf) H(i,j)*B(m-i,n-j) ------ (* means multiplication now)
>
> and for L(1,1), where I would get B to be a 3x3 filter:
>
> L(1,1) = H[0,0]*B[1,1] + H[1,0]*B[0,1] + H[2,0]B[-1,1]...
>
> So here, we take into account having to sum all the elements around the pixel of interest L(1,). I am supposed to construct this for all pixels of this image - 96x96.
>
> Then having L = H*B equations, I think of these as Ax = b and then use what you suggested to get the values of A, which in this case is my H. I know what the L values are and can opt for a guess of B, so how do I create these equation in matlab? I cannot come up with 96x96 equations by hand.
=====================

I still don't understand what's forcing you to construct these equations when you have a closed form formula for the solution:

zp=96*2; %zero-padding

H=ifft2(fft2(L,zp,zp)./fft2(B,zp,zp));

H=H(1:96,1:96); %unzero-pad

I thin it is really cool I got to know about some of the things you are mentioning, but since this is an instructional problem to solve, my teacher said something like:

"It would be instructional for you to see what the H matrix would look like"

Ok, so assuming I do what you suggest, and the fact that I have L and I can create a B using fspecial. If the equation I am dealing with is L = H*B, then you are saying I can find H from L and B by using the Fourier functions and/or doing H = L/B?

Ok, so I tried using the Fourier functions, and I get an error. I first created filter using fspecial of type Gaussian of the same size as L (96x96). So I want L = H*B, and not that I have B, and since I have L, then I want to find H. Therefore, I did:

>> fft2(H)=fft2(im)./fft2(B);
Warning: Divide by zero.
??? Subscript indices must either be real positive integers or logicals.
 
As you can see below, when I look for the minimum in the columns, I see that there are many zeros, and this is giving me the error, apparently, so what do I do now?
>> min(im(:,:))

ans =

  Columns 49 through 64

   21 6 0 2 1 2 1 4 0 0 2 7 8 10 5 1

On Feb 8, 1:56 pm, "B C" <bcorte...@gmail.com> wrote:
> I thin it is really cool I got to know about some of the things you are mentioning, but since this is an instructional problem to solve, my teacher said something like:
>
> "It would be instructional for you to see what the H matrix would look like"
>
> Ok, so assuming I do what you suggest, and the fact that I have L and I can create a B using fspecial. If the equation I am dealing with is L = H*B, then you are saying I can find H from L and B by using the Fourier functions and/or doing H = L/B?

------------------------------------------------------------------------------
No, that is NOT what he is saying. Look -- he has fft in there but
you don't. Basically you're trying to do what is called an "inverse
filter." Where
fft(H) = fft(L) ./ fft(B)
then to get H you have to inverse fft that. Just what Matt J showed.
This is not used in practice because it's so lousy. A better way is
to use a Wiener (or Wiener-Helstrom) filter which is more stable
because it takes into account the spectrum of the noise (which is
always present in real life). But it may be instructive for you if
you're looking at noise-free synthetic signals, or just to see how it
goes haywire if you add a little bit of noise.

Hi there, well, conv2 give one the convolution of 2 matrices, right? But in my case for the convolution L = H*B, I already have the convolution of those matrices - H and B, which is my L. If I create a filter matrix B using fspecial, I want to be able to get H from the L and B I already know, but I don't think conv2 does this. There is the deconv function, but this function results in the deconvolution of some A from B, and in my case, I get this error:

>> deconv(L,H)
??? Undefined function or method 'filter' for input arguments of type 'uint8'.

Any ideas?

"B C" <bcortez23@gmail.com> wrote in message <hkpnco$5nc$1@fred.mathworks.com>...
> Ok, so I tried using the Fourier functions, and I get an error. I first created filter using fspecial of type Gaussian of the same size as L (96x96). So I want L = H*B, and not that I have B, and since I have L, then I want to find H. Therefore, I did:
>
> >> fft2(H)=fft2(im)./fft2(B);
> Warning: Divide by zero.
> ??? Subscript indices must either be real positive integers or logicals.
> As you can see below, when I look for the minimum in the columns, I see that there are many zeros, and this is giving me the error, apparently, so what do I do now?
> >> min(im(:,:))

First of all, you skipped the zero padding, which was important. Secondly
min(im) is not relevant because you are not dividing by im. You are dividing by
fft2(B). Clearly, fft2(B) has zeros in it somewhere. If your filter spectrum is not bounded well away from zero, then, L=H*B has no unique solution, and
as ImageAnalyst was saying, needs to be modified/regularized. Deriving a Wiener filter from B is one way to do so, but your prof. may have intended something else for you.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hkpoeh$b5v$1@fred.mathworks.com>...
> "B C" <bcortez23@gmail.com> wrote in message <hkpnco$5nc$1@fred.mathworks.com>...
> > Ok, so I tried using the Fourier functions, and I get an error. I first created filter using fspecial of type Gaussian of the same size as L (96x96). So I want L = H*B, and not that I have B, and since I have L, then I want to find H. Therefore, I did:
> >
> > >> fft2(H)=fft2(im)./fft2(B);
> > Warning: Divide by zero.
> > ??? Subscript indices must either be real positive integers or logicals.
> > As you can see below, when I look for the minimum in the columns, I see that there are many zeros, and this is giving me the error, apparently, so what do I do now?
> > >> min(im(:,:))
>
> First of all, you skipped the zero padding, which was important. Secondly
> min(im) is not relevant because you are not dividing by im. You are dividing by
> fft2(B). Clearly, fft2(B) has zeros in it somewhere. If your filter spectrum is not bounded well away from zero, then, L=H*B has no unique solution, and
> as ImageAnalyst was saying, needs to be modified/regularized. Deriving a Wiener filter from B is one way to do so, but your prof. may have intended something else for you.
===================================================
Hi there, I did do it, let me show you:

L = imread('testIm1.jpg');
[r,c] = size(im); % image is 96x96
B = fspecial('gaussian', [r,c], 0.9); % gaussian filter

zp=96*2; %zero-padding

H=ifft2(fft2(L,zp,zp)./fft2(B,zp,zp));

H=H(1:96,1:96); %unzero-pad

%Subsampled because I want to use other matrices to construct one larger one
n = 2;
output = H(1:n:end, 1:n:end);
subSam_size = size(output);
figure(1)
imshow(output)
title('Sub-sample image')
xlabel(int2str(subSam_size))

PROBLEM IS: The result from ifft2 gives me a black screen, and this can't be right.

"B C" <bcortez23@gmail.com> wrote in message <hkpp8m$247$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hkpoeh$b5v$1@fred.mathworks.com>...
> > "B C" <bcortez23@gmail.com> wrote in message <hkpnco$5nc$1@fred.mathworks.com>...
> > > Ok, so I tried using the Fourier functions, and I get an error. I first created filter using fspecial of type Gaussian of the same size as L (96x96). So I want L = H*B, and not that I have B, and since I have L, then I want to find H. Therefore, I did:
> > >
> > > >> fft2(H)=fft2(im)./fft2(B);
> > > Warning: Divide by zero.
> > > ??? Subscript indices must either be real positive integers or logicals.
> > > As you can see below, when I look for the minimum in the columns, I see that there are many zeros, and this is giving me the error, apparently, so what do I do now?
> > > >> min(im(:,:))
> >
> > First of all, you skipped the zero padding, which was important. Secondly
> > min(im) is not relevant because you are not dividing by im. You are dividing by
> > fft2(B). Clearly, fft2(B) has zeros in it somewhere. If your filter spectrum is not bounded well away from zero, then, L=H*B has no unique solution, and
> > as ImageAnalyst was saying, needs to be modified/regularized. Deriving a Wiener filter from B is one way to do so, but your prof. may have intended something else for you.
> ===================================================
> Hi there, I did do it, let me show you:
>
> L = imread('testIm1.jpg');
> [r,c] = size(im); % image is 96x96
> B = fspecial('gaussian', [r,c], 0.9); % gaussian filter
>
> zp=96*2; %zero-padding
>
> H=ifft2(fft2(L,zp,zp)./fft2(B,zp,zp));

Very well. But again, the solution assumes that
fft2(B,zp,zp))
is bounded away from zero. Because of the warning message, we know that this is not the case, and that your problem is ill-posed. You need to reformulate it so that a solution actually exists.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hkpqgl$ipc$1@fred.mathworks.com>...
> "B C" <bcortez23@gmail.com> wrote in message <hkpp8m$247$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hkpoeh$b5v$1@fred.mathworks.com>...
> > > "B C" <bcortez23@gmail.com> wrote in message <hkpnco$5nc$1@fred.mathworks.com>...
> > > > Ok, so I tried using the Fourier functions, and I get an error. I first created filter using fspecial of type Gaussian of the same size as L (96x96). So I want L = H*B, and not that I have B, and since I have L, then I want to find H. Therefore, I did:
> > > >
> > > > >> fft2(H)=fft2(im)./fft2(B);
> > > > Warning: Divide by zero.
> > > > ??? Subscript indices must either be real positive integers or logicals.
> > > > As you can see below, when I look for the minimum in the columns, I see that there are many zeros, and this is giving me the error, apparently, so what do I do now?
> > > > >> min(im(:,:))
> > >
> > > First of all, you skipped the zero padding, which was important. Secondly
> > > min(im) is not relevant because you are not dividing by im. You are dividing by
> > > fft2(B). Clearly, fft2(B) has zeros in it somewhere. If your filter spectrum is not bounded well away from zero, then, L=H*B has no unique solution, and
> > > as ImageAnalyst was saying, needs to be modified/regularized. Deriving a Wiener filter from B is one way to do so, but your prof. may have intended something else for you.
> > ===================================================
> > Hi there, I did do it, let me show you:
> >
> > L = imread('testIm1.jpg');
> > [r,c] = size(im); % image is 96x96
> > B = fspecial('gaussian', [r,c], 0.9); % gaussian filter
> >
> > zp=96*2; %zero-padding
> >
> > H=ifft2(fft2(L,zp,zp)./fft2(B,zp,zp));
>
> Very well. But again, the solution assumes that
> fft2(B,zp,zp))
> is bounded away from zero. Because of the warning message, we know that this is not the case, and that your problem is ill-posed. You need to reformulate it so that a solution actually exists.

==============================
Thank for your time and comments Matt, I didn't think this would be so hard. I am learning LSE and DSV, and so I thought it would be easy to use matlab to create all equations to form Ax=b, but this seems to not be the best way. Now that I have tried some of the things you and the other person have suggested, I would like to look further into this. I was reading about convolution and the whole thing about rotating the filter mask seems complicated to do from scratch, so I think I would like to use these Fourier functions. One thing that I am not understanding is the size of the filter - in my case B. I was reading that the filter can be any size, that is a could have some 3x3 filter and use that in the convolution, and in my code above you can see that B is actually rxc. So it is not 3x3, instead it is 96x96. This is when I get the black image, but if I change it to 3x3, then I see
something, BUT NOTHING that looks like L. If I need to pose my problem correctly, what am I missing? If B is 3x3, then fft2 pads the 3x3 with zeros up to r and c dimensions. If B is 96x96, fft2 pads it to 96*2, so how can I bound away from zero like you suggest?

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <900d5bfc-5f97-416e-be98-686ecf04c367@n33g2000yqb.googlegroups.com>...
> On Feb 8, 1:56 pm, "B C" <bcorte...@gmail.com> wrote:
> > I thin it is really cool I got to know about some of the things you are mentioning, but since this is an instructional problem to solve, my teacher said something like:
> >
> > "It would be instructional for you to see what the H matrix would look like"
> >
> > Ok, so assuming I do what you suggest, and the fact that I have L and I can create a B using fspecial. If the equation I am dealing with is L = H*B, then you are saying I can find H from L and B by using the Fourier functions and/or doing H = L/B?
>
> ------------------------------------------------------------------------------
> No, that is NOT what he is saying. Look -- he has fft in there but
> you don't. Basically you're trying to do what is called an "inverse
> filter." Where
> fft(H) = fft(L) ./ fft(B)
> then to get H you have to inverse fft that. Just what Matt J showed.
> This is not used in practice because it's so lousy. A better way is
> to use a Wiener (or Wiener-Helstrom) filter which is more stable
> because it takes into account the spectrum of the noise (which is
> always present in real life). But it may be instructive for you if
> you're looking at noise-free synthetic signals, or just to see how it
> goes haywire if you add a little bit of noise.
===============================================
When I was reading the descriptions of some of the convolution functions in matlab, I remember reading something about regularized, and so I went and looked for it and found it. So this is my code for that:

L = imread('test1.jpg');
B = fspecial('gaussian', [3,3], 0.9);
H = deconvreg(L,B);
imshow(H)

This give me an ok H, but I think it works for now. However, I wanted to know if you could explain what the regularization means for me in terms of an image matrix or in this case the deconvolution of an image matrix with a blur Gaussian matrix. Thank you, Bc

"B C" <bcortez23@gmail.com> wrote in message <hkpse0$gd2$1@fred.mathworks.com>...

 I was reading about convolution and the whole thing about rotating the filter mask seems complicated to do from scratch, so I think I would like to use these Fourier functions.
=================
It's also a very inefficient solution.


 One thing that I am not understanding is the size of the filter - in my case B. I was reading that the filter can be any size, that is a could have some 3x3 filter and use that in the convolution, and in my code above you can see that B is actually rxc. So it is not 3x3, instead it is 96x96. This is when I get the black image, but if I change it to 3x3, then I see
> something, BUT NOTHING that looks like L. If I need to pose my problem correctly, what am I missing? If B is 3x3, then fft2 pads the 3x3 with zeros up to r and c dimensions.
=======================
It's more of an issue of what B's spectrum looks like. For example, if B=0 then it doesn't matter whether is is 3x3 or rxc. There is no solution to the problem.


>so how can I bound away from zero like you suggest?
===================

As ImageAnalyst suggested, but posing the problem as a Wiener filtering of L.

 

Ok, great, so I think I am seeing results. I used both deconvreg and deconvvwnr and compared the results, and they seem slightly different, but at least I am getting somewhere. Now, as you saw from my code:

n = 2;
output = H(1:n:end, 1:n:end);
subSam_size = size(output);

I subsampled because the purpose is that if my original matrix is 96x96, I can subsample, let say 8 L matrices as above to get 8 48x48 matrices. Now, what I would like to be able to do is combine these matrices, so my question is: If I have the 8, 48x48 matrices, and I want to combine them, should I subsample at different n values? like decimals? so that I can merger the matrices and there is no overlap in them?

If you know of some reference to this, I would really appreciate some guidance.

Thanks,

Bc

Bc: When you display real arrays, I usually use
imshow(imageArray, []);
If your numbers are small, like in the range of 0-1 which MATLAB likes
to use a lot, then it may appear black.

Matt: Any chance ifft() would return a complex number array instead
of a real valued array? Or is it always going to be guaranteed to be
real, not complex?

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <83161cf8-b358-4627-8640-35fd0bfdaecd@19g2000yql.googlegroups.com>...
 
> Matt: Any chance ifft() would return a complex number array instead
> of a real valued array? Or is it always going to be guaranteed to be
> real, not complex?

It would definitely return a result with imaginary-valued round-off noise. You would need to strip that away.

"Bc " <bcortez23@gmail.com> wrote in message <hkpup0$c57$1@fred.mathworks.com>...

> I subsampled because the purpose is that if my original matrix is 96x96, I can subsample, let say 8 L matrices as above to get 8 48x48 matrices. Now, what I would like to be able to do is combine these matrices, so my question is: If I have the 8, 48x48 matrices, and I want to combine them, should I subsample at different n values? like decimals? so that I can merger the matrices and there is no overlap in them?
================

I didn't really understand that. Which matrices are 48x48 and which are 96x96? Is it the original matrix H is 96x96 but you are trying to restore it from several observations of L which is 48x48? So this is a super-resolution problem?

"Bc " <bcortez23@gmail.com> wrote in message <hkptjh$rhk$1@fred.mathworks.com>...

However, I wanted to know if you could explain what the regularization means for me in terms of an image matrix or in this case the deconvolution of an image matrix with a blur Gaussian matrix.
========================

They are basically doing what I said earlier in Message #2,

iift2(H)=fft2(L)./(fft2(B)+beta*fft2(C)); %assumes everything is zero-padded

where C is a high-pass filter. The Wiener filter fits into this family where
beta*fft2(C) is the known noise-spectrum (if you know it).

However, deconvreg gives you other options for C, which you can read about by doing

>>doc deconvreg

It explains that the default C is a Laplacian operator, which is a filter such that
conv2(C,H) approximates the 2nd derivative of H.

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hkq27l$k1c$1@fred.mathworks.com>...
> "Bc " <bcortez23@gmail.com> wrote in message <hkpup0$c57$1@fred.mathworks.com>...
>
> > I subsampled because the purpose is that if my original matrix is 96x96, I can subsample, let say 8 L matrices as above to get 8 48x48 matrices. Now, what I would like to be able to do is combine these matrices, so my question is: If I have the 8, 48x48 matrices, and I want to combine them, should I subsample at different n values? like decimals? so that I can merger the matrices and there is no overlap in them?
> ================
>
> I didn't really understand that. Which matrices are 48x48 and which are 96x96? Is it the original matrix H is 96x96 but you are trying to restore it from several observations of L which is 48x48? So this is a super-resolution problem?

=======================================================
Yes Matt, I am doing SR for an image. I have slightly, very slightly different images, and so I wanted to create the HR image from several LR images by solving L=H*B subsampled. My original images are 96x96 and if I subsample as shown in my code before, then I will have images that are 48x48.

Now, if I had actually created these equation for the pixels, then I was expecting to do this for several images (lets say 8), so that I would have 8 subsample images which I could then "put together" to created a "larger" version of the same image. What do you think?

I'm going to to try the fft2 function again while I get some feedback. Thanks.

"Bc " <bcortez23@gmail.com> wrote in message <hkqjdj$b54$1@fred.mathworks.com>...

> Yes Matt, I am doing SR for an image. I have slightly, very slightly different images, and so I wanted to create the HR image from several LR images by solving L=H*B subsampled. My original images are 96x96 and if I subsample as shown in my code before, then I will have images that are 48x48.
>
> Now, if I had actually created these equation for the pixels, then I was expecting to do this for several images (lets say 8), so that I would have 8 subsample images which I could then "put together" to created a "larger" version of the same image. What do you think?
=================

It depends on what the differences are among the 8 observations. If they all differ by a known translation than taking the Fourier transform will have an easily processible effect