The problem is to simulate an impulse function convolving with a 2D Gaussian
function in matlab. So discretization is the first step. Suppose there is a
0.4m x 0.4m block, the impulse function is located at (1/6*0.4m,
1/2*0.4m)=(0.0666m, 0.2m) on this 0.4m x 0.4m block. The Guassian function
is assumed to centered at origin and its sigma is 1m. After convolution, we
are going to find the result value at points, e.g. (0.1m, 0.1m), (0.3m,
0.3m), etc.
I want to ask if the sampling step size for this 0.4m x 0.4m block, where
the impulse function is locaed, should be the same as the sampling step size
for the Gaussian function?
I assumed the step size of the impulse function's 0.4m x 0.4m block should
be the same as the step size of the Gaussian function. Am I correct?
The problem I met is: in order to capture fully the energy of the Gaussian
function, suppose I sample the Gaussian function into [-N .. N, -N .. N]
grids with certain stepsize,
then the following should be satisfied. (maybe 99.5% energy will be
maintained by 4*sigma):
-----> N*stepsize > 4*sigma
In order to make the Gaussian function small, I should try to minimize N
under the above condition, the Matlab can handle N = several hundreds
reasonably fast.
In order to reduce N, the stepsize should be large.
Assuming the plane where impulse function resides should be sampled at the
same stepsize as this Gaussian function,
The 0.4m x 0.4m block will be sample very rough since the stepsize is large.
So I begin to doubt if my sampling approach is correct.
%-------------- Source code ------------------------------
tt=zeros(11, 11);
%using 10 samples for the 0.4m x 0.4m block where the impulse function
resides
tt(floor(10/2)+1, floor(10/6)+1)=1;
%place an impulse function at (1/6*0.4m, 1/2*0.4m) of the0.4m x0.4m block
uu=mygaussian([201, 201], 0.04, 1);
%the sampling stepsize is 0.4m/10=0.04m, the sigma is 1m. take 2N+1=201
%N=100, 100*0.04 > = 4*sigma
tic
tt1=conv2(tt, uu, 'same');
%convolution in Matlab
toc
tt1(floor(10/4)+1, floor(10/4)+1)
%result value at (0.1m, 0.1m) in the 0.4m x 0.4m block
tt1(floor(30/4)+1, floor(30/4)+1)
%result value at (0.3m, 0.3m) in the 0.4m x 0.4m block
function h=mygaussian(p2, p3, p4);
% p3 is stepsize, p4 is sigma, p2 is [2N+1, 2N+1]
siz = (p2-1)/2;
std = p4;
[x,y] = meshgrid(-p3*siz(2):p3:p3*siz(2),-p3*siz(1):p3:p3*siz(1));
arg = -(x.*x + y.*y)/(2*std*std);
h = exp(arg)/(2*pi*std*std);
%----------Execution results---------------------
Elapsed time is 0.015000 seconds.
ans =
0.1579
ans =
0.1541
Correct answer should be:
>> fz1(0.1, 0.1)
ans =
0.1583
>> fz1(0.3, 0.3)
ans =
0.1541
where fz1= 1/(2*pi)*exp(-0.5*((x-1/6*0.4).^2+(y-0.5*0.4).^2));
So the result value at (0.1m, 0.1m) is 0.25% off. the result at (0.3m, 0.3m)
is correct.
So I should increase the number of samples on the 0.4m x 0.4m block, by
decreasing the stepsize.
But decreasing the stepsize will increase the number of samples on the
Gaussian function, the performance will be severely degraded if we go from
201 to 801, etc.
Please tell me is there any way that I can set a good number of samples on
the 0.4m x 0.4 m block where the impulse function resides, by theoratical
derivation? Instead of my trial and error?
Any thoughts? Please help me! Your advice/suggestions/comments are highly
appreciated!
-Lucy
Difficult problem: how to simulate an impulse response convolve with a 2D Gaussian function in Matlab?
Started by ●August 16, 2004
Reply by ●August 16, 20042004-08-16
lucy wrote: ...> Please tell me is there any way that I can set a good number of samples on > the 0.4m x 0.4 m block where the impulse function resides, by theoratical > derivation? Instead of my trial and error? > > Any thoughts? Please help me! Your advice/suggestions/comments are highly > appreciated! > > -LucyA 2-D Gaussian kernel is separable, which greatly simplifies the computation. (The computation becomes O(N) instead of O(N�)). That's a major reason for choosing it in the first place. Take advantage of that. Jerry -- Where are we going? Why are we in this hand basket? ������������������������������������������������������������������������
Reply by ●August 16, 20042004-08-16
Jerry Avins <jya@ieee.org> writes:> [...] > A 2-D Gaussian kernel is separable, which greatly simplifies the > computation. (The computation becomes O(N) instead of O(N�)). That's a > major reason for choosing it in the first place. Take advantage of that.Nice catch, Jerry! -- Randy Yates Sony Ericsson Mobile Communications Research Triangle Park, NC, USA randy.yates@sonyericsson.com, 919-472-1124
Reply by ●August 16, 20042004-08-16
"Jerry Avins" <jya@ieee.org> wrote in message news:412107f9$0$5922$61fed72c@news.rcn.com...> lucy wrote: > > ... > > > Please tell me is there any way that I can set a good number of sampleson> > the 0.4m x 0.4 m block where the impulse function resides, bytheoratical> > derivation? Instead of my trial and error? > > > > Any thoughts? Please help me! Your advice/suggestions/comments arehighly> > appreciated! > > > > -Lucy > > A 2-D Gaussian kernel is separable, which greatly simplifies the > computation. (The computation becomes O(N) instead of O(N�)). That's a > major reason for choosing it in the first place. Take advantage of that. > > Jerry > -- > Where are we going? Why are we in this hand basket? > ������������������������������������������������������������������������ >Hi Jerry, This is a nice hint! Thank you very much! For one impulse function convolving with 2D Gaussian, I can do 1D then multiplying out the 2D result from the 1D result... that's ok. But what if there is an array of impulse functions, scattered in a 2D comb fashion, maybe even diamond shaped 2D comb fashion, can this impulse function array convolving with 2D gaussian function still be splitted into 1D impulse array convolving with 1D Gaussian, then multiply out the 2D result? I am not sure. What do you think? Any more thoughts? Thank you a lot!
Reply by ●August 16, 20042004-08-16
"Jerry Avins" <jya@ieee.org> wrote in message news:412107f9$0$5922$61fed72c@news.rcn.com...> lucy wrote: > > ... > > > Please tell me is there any way that I can set a good number of sampleson> > the 0.4m x 0.4 m block where the impulse function resides, bytheoratical> > derivation? Instead of my trial and error? > > > > Any thoughts? Please help me! Your advice/suggestions/comments arehighly> > appreciated! > > > > -Lucy > > A 2-D Gaussian kernel is separable, which greatly simplifies the > computation. (The computation becomes O(N) instead of O(N�)). That's a > major reason for choosing it in the first place. Take advantage of that. > > Jerry > -- > Where are we going? Why are we in this hand basket? > ������������������������������������������������������������������������ >Hi Jerry, thank you for your hint. I've found the following: So the procedure is first to do two 1D Guassian convolution saperately, then do bilinear interpolation... Don't know how approximate and how inaccuracy it will be yet. Will try and see... However, I think this works for 1 impluse function convolving with the Gaussian function, not sure if it is applicable for an array of impulse functions convolving with the Gaussian function, esp. if the array of impulse functions does not align with x and y axes, instead, it's in diamond shaped comb fashion... ---------------------------------------------------------------------------- --------- If bilinear is "yes" and the major axis of the Gaussian kernel is aligned along either the x or y axis, GAUSS uses the fact that the Gaussian function is mathematically separable (bilinear) in x and y to speed up the convolution process. A bilinear 2D convolution kernel in x and y is one which can be separated into two equivalent 1D convolution kernels in x and y respectively. Although the bilinear approximation and the full 2D convolution are mathematically equivalent, the user will actually see SMALL differences between an image convolved with the full 2D kernel and the same image convolved with the equivalent bilinear kernel. These differences are the result of the finite size of the convolution kernel (the integration does not extend to infinity in either direction), and the fact that off-axis kernel elements outside the nsigma limit cannot be set to 0 in the bilinear case as they are in the full 2D case. Therefore the bilinear kernel is less radially symmetric than the full 2D kernel. In most cases the differences are small and more than made up for by the greatly decreased execution time.
Reply by ●August 16, 20042004-08-16
lucy wrote:> > The problem is to simulate an impulse function convolving with a 2D Gaussian > function in matlab.I didn't wade thru everything below but what you describe above is the process of multiplying every coefficient in your 2d array by the number 1 (or maybe some other number if the magnitude of the impulse is not 1). The only issue is where and how your array is stored before and after. -jim>So discretization is the first step. Suppose there is a > 0.4m x 0.4m block, the impulse function is located at (1/6*0.4m, > 1/2*0.4m)=(0.0666m, 0.2m) on this 0.4m x 0.4m block. The Guassian function > is assumed to centered at origin and its sigma is 1m. After convolution, we > are going to find the result value at points, e.g. (0.1m, 0.1m), (0.3m, > 0.3m), etc. > > I want to ask if the sampling step size for this 0.4m x 0.4m block, where > the impulse function is locaed, should be the same as the sampling step size > for the Gaussian function? > > I assumed the step size of the impulse function's 0.4m x 0.4m block should > be the same as the step size of the Gaussian function. Am I correct? > > The problem I met is: in order to capture fully the energy of the Gaussian > function, suppose I sample the Gaussian function into [-N .. N, -N .. N] > grids with certain stepsize, > > then the following should be satisfied. (maybe 99.5% energy will be > maintained by 4*sigma): > > -----> N*stepsize > 4*sigma > > In order to make the Gaussian function small, I should try to minimize N > under the above condition, the Matlab can handle N = several hundreds > reasonably fast. > > In order to reduce N, the stepsize should be large. > > Assuming the plane where impulse function resides should be sampled at the > same stepsize as this Gaussian function, > > The 0.4m x 0.4m block will be sample very rough since the stepsize is large. > > So I begin to doubt if my sampling approach is correct. > > %-------------- Source code ------------------------------ > tt=zeros(11, 11); > %using 10 samples for the 0.4m x 0.4m block where the impulse function > resides > tt(floor(10/2)+1, floor(10/6)+1)=1; > %place an impulse function at (1/6*0.4m, 1/2*0.4m) of the0.4m x0.4m block > uu=mygaussian([201, 201], 0.04, 1); > %the sampling stepsize is 0.4m/10=0.04m, the sigma is 1m. take 2N+1=201 > %N=100, 100*0.04 > = 4*sigma > tic > tt1=conv2(tt, uu, 'same'); > %convolution in Matlab > toc > tt1(floor(10/4)+1, floor(10/4)+1) > %result value at (0.1m, 0.1m) in the 0.4m x 0.4m block > tt1(floor(30/4)+1, floor(30/4)+1) > %result value at (0.3m, 0.3m) in the 0.4m x 0.4m block > > function h=mygaussian(p2, p3, p4); > % p3 is stepsize, p4 is sigma, p2 is [2N+1, 2N+1] > > siz = (p2-1)/2; > std = p4; > > [x,y] = meshgrid(-p3*siz(2):p3:p3*siz(2),-p3*siz(1):p3:p3*siz(1)); > arg = -(x.*x + y.*y)/(2*std*std); > > h = exp(arg)/(2*pi*std*std); > > %----------Execution results--------------------- > > Elapsed time is 0.015000 seconds. > > ans = > > 0.1579 > > ans = > > 0.1541 > > Correct answer should be: > > >> fz1(0.1, 0.1) > > ans = > > 0.1583 > > >> fz1(0.3, 0.3) > > ans = > > 0.1541 > > where fz1= 1/(2*pi)*exp(-0.5*((x-1/6*0.4).^2+(y-0.5*0.4).^2)); > > So the result value at (0.1m, 0.1m) is 0.25% off. the result at (0.3m, 0.3m) > is correct. > > So I should increase the number of samples on the 0.4m x 0.4m block, by > decreasing the stepsize. > > But decreasing the stepsize will increase the number of samples on the > Gaussian function, the performance will be severely degraded if we go from > 201 to 801, etc. > > Please tell me is there any way that I can set a good number of samples on > the 0.4m x 0.4 m block where the impulse function resides, by theoratical > derivation? Instead of my trial and error? > > Any thoughts? Please help me! Your advice/suggestions/comments are highly > appreciated! > > -Lucy-----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
Reply by ●August 16, 20042004-08-16
"jim" <"N0sp"@m.sjedging@mwt.net> wrote in message news:4121150a_3@corp.newsgroups.com...> > > lucy wrote: > > > > The problem is to simulate an impulse function convolving with a 2DGaussian> > function in matlab. > > I didn't wade thru everything below but what you describe above is the > process of multiplying every coefficient in your 2d array by the number > 1 (or maybe some other number if the magnitude of the impulse is not 1). > The only issue is where and how your array is stored before and after. > > -jim > >Hi Jim, yeah, you've seen through everything. For one impulse function convolving with the Gaussian function, it is nothing but a storage problem(copying the function values from the filter to the image plane where the impulse function resides...) But as you can see from my followup posts: there will be an array of impulse function convolving with Gaussian function and overlapping will occur and there is a serious problem of tradeoff the sampling grid size and the accuracy... that's why I was asking the questions... any more thoughts?
Reply by ●August 16, 20042004-08-16
lucy wrote:> > "jim" <"N0sp"@m.sjedging@mwt.net> wrote in message > news:4121150a_3@corp.newsgroups.com... > > > > > > lucy wrote: > > > > > > The problem is to simulate an impulse function convolving with a 2D > Gaussian > > > function in matlab. > > > > I didn't wade thru everything below but what you describe above is the > > process of multiplying every coefficient in your 2d array by the number > > 1 (or maybe some other number if the magnitude of the impulse is not 1). > > The only issue is where and how your array is stored before and after. > > > > -jim > > > > > > Hi Jim, yeah, you've seen through everything. For one impulse function > convolving with the Gaussian function, it is nothing but a storage > problem(copying the function values from the filter to the image plane where > the impulse function resides...) > > But as you can see from my followup posts: there will be an array of impulse > function convolving with Gaussian function and overlapping will occur and > there is a serious problem of tradeoff the sampling grid size and the > accuracy... that's why I was asking the questions... any more thoughts?I take it you are implementing a gaussian blur with adjustable radius. Separating into 1d kernels may help if the size is large, but modern computers are bound to there cache behavior so that can be a problem with separate 1d filters. The other issue you will face is - what are you going to do when your kernel hits the boundaries of the image? -jim -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
Reply by ●August 17, 20042004-08-17
"jim" <"N0sp"@m.sjedging@mwt.net> wrote in message news:41211acb_3@corp.newsgroups.com...> > > lucy wrote: > > > > "jim" <"N0sp"@m.sjedging@mwt.net> wrote in message > > news:4121150a_3@corp.newsgroups.com... > > > > > > > > > lucy wrote: > > > > > > > > The problem is to simulate an impulse function convolving with a 2D > > Gaussian > > > > function in matlab. > > > > > > I didn't wade thru everything below but what you describe above is the > > > process of multiplying every coefficient in your 2d array by thenumber> > > 1 (or maybe some other number if the magnitude of the impulse is not1).> > > The only issue is where and how your array is stored before and after. > > > > > > -jim > > > > > > > > > > Hi Jim, yeah, you've seen through everything. For one impulse function > > convolving with the Gaussian function, it is nothing but a storage > > problem(copying the function values from the filter to the image planewhere> > the impulse function resides...) > > > > But as you can see from my followup posts: there will be an array ofimpulse> > function convolving with Gaussian function and overlapping will occurand> > there is a serious problem of tradeoff the sampling grid size and the > > accuracy... that's why I was asking the questions... any more thoughts? > > I take it you are implementing a gaussian blur with adjustable radius. > Separating into 1d kernels may help if the size is large, but modern > computers are bound to there cache behavior so that can be a problem > with separate 1d filters. The other issue you will face is - what are > you going to do when your kernel hits the boundaries of the image? > > -jim >Hi separable 1d filter is really faster than 2D, the performance is 1000 times faster... but I really don't know what I am going to do when my kernel hits the boundary of the image. I used "conv2(A, B, 'same') in Matlab, it did not give any error message so I have even not noticed the problem... what do you think is the solution? thanks a lot,
Reply by ●August 17, 20042004-08-17
lucy wrote: ...> but I really don't know what I am going to do when my kernel hits the > boundary of the image.... This is a general problem. It doesn't depend on 1-D, 2-D, or the method chosen. There are several solutions. (When many ways are listed, you know immediately that none are very good.) Here are some: 1) Assume a value for off-the-edge samples (black, white, gray, etc.). 2) Allow the data set to shrink with each pass. 3) Use modified kernels near the edge that are more appropriate to their positions. 4) Compute probable values for off-the-edge samples. You can show that this can be equivalent to (3). To repeat, the problem is inherent in the process of convolution. It goes unnoticed in audio work because most material starts and ends with silence, and so extending the data with null samples works fine. When the data don't fit that pattern, the result becomes what we call a "transient". Jerry -- ... the worst possible design that just meets the specification - almost a definition of practical engineering. .. Chris Bore ������������������������������������������������������������������������
