| Convolution [message #26609] |
Tue, 11 September 2001 05:22  |
Kay Bente
Messages: 9
Registered: September 2001
|
Junior Member |
|
|
Hi
I have to convolute a 256x256x128 Floating Point array with a 3D Gaussian
Kernel of ~ 30x30x30, this lasts round about 45Minutes. So my question is,
if there is any way how i can speed this up. I tried to separate this in
each dimension with a 1D Kernel, but I don�t know if I have done this
correct (cause the procedure hangs up after a few loops)
I know that the Convolution of two functions is a Multiplication in Fourier
Space, but how can I do this with discrete arrays, do I have to enlarge my
kernel to the size of the array i want to smooth? If so, the creation of the
kernel with the dimensions of my array nearly lasts as long as the normal
convolution :-(
So i would be very happy about any hints!!
This is how I tried it
FUNCTION Gauss3D, mat, n, s
;n = Size Of Array
;s = Sigma
;Gauss Kernel is completely circularly symmetric operator
;=> Convol in each dimension with a 1D Kernel
;Create Kernel
a = DindGen(n)-(n/2)
print,a
kernel = Exp(-(Temporary(a))^2/(2.*s^2))/(Sqrt(2.*!PI)*s)
tstart=SysTime(1)
s=Size(mat,/Dimensions)
mat2=FltArr(s(0),s(1),s(2))
;Convol in X,Y
FOR i=0, s(2) DO BEGIN
buffer=mat(*,*,i)
buffer1=Convol(buffer,kernel,Total(kernel))
buffer1=Convol(buffer1,Rotate(kernel,1),Total(kernel))
mat2(*,*,i)=Temporary(buffer1)
ENDFOR
;Convol in Z
FOR i=0, s(1) DO BEGIN
buffer=mat2(*,i,*)
buffer1=Convol(buffer,kernel,Total(kernel))
mat2(*,i,*)=Temporary(buffer)
ENDFOR
Return, mat2
END
---
Posted via news://freenews.netfront.net
Complaints to news@netfront.net
|
|
|
|
| Re: convolution [message #34607 is a reply to message #26609] |
Fri, 28 March 2003 08:17  |
Chris[1]
Messages: 23
Registered: January 2003
|
Junior Member |
|
|
Hi Larry;
Someday I'm going to get IDL installed on the same machine I write mail
from, but till then....
W/r to your problem, formally one of the forward fft's should be conjugated;
and when you plot the results, plot the absolute values, or explicitly the
real and complex parts, and see if that helps.
Chris
"Larry Morgan" <lkm8@ukc.ac.uk> wrote in message
news:b61e00$khl$1@athena.ukc.ac.uk...
> Hi,
> I am at a loss to explain the output from the program below and although
> it's not strictly an idl problem I was wondering if anyone could help me.
> I want to convolve the two functions in the left half of the plot window
> together. When I multiply their fourier transforms together and inverse
> transform the result back I get what appears in the right hand window.
This
> is not at all what I expected although from everything I've read there is
> nothing wrong with the method I have used.
> Can anyone help me?
> cheers
> Larry
>
> Pro convolve
>
> xxx=((DINDGEN(20000))*0.01)
>
>
beamlong=(0.960944*exp(-((xxx-100.)^2)/(2*(10.6337^2))))+(0. 0390565*exp(-((x
xx-100.)^2)/(2*(33.2939^2))))
>
> loz_850=0.00060018403/(4.0*(((xxx-100.)/27.269890)^2.0)+1.0)
>
> loz_850=loz_850/max(loz_850)
>
> !p.multi=[0,2,1,0,0]
> plot,xxx,loz_850,linestyle=1
> oplot,xxx,beamlong
> imconv_850=fft(fft(loz_850)*fft(beamlong),/inverse)
>
> plot,xxx,imconv_850
> !p.multi=0
> end
>
|
|
|
|
|
|
|
|
| Re: convolution [message #60323 is a reply to message #26609] |
Tue, 20 May 2008 04:31 |
Chris[5]
Messages: 16
Registered: May 2008
|
Junior Member |
|
|
On May 20, 12:49 am, Chris <cnb4s...@gmail.com> wrote:
> On May 19, 12:32 pm, sarah <sarahwiddeco...@yahoo.com> wrote:
>
>
>
>> On May 12, 9:31 pm, David Fanning <n...@dfanning.com> wrote:
>
>>> sarah writes:
>>>> x=make_array(1024)
>>>> sigma=15
>>>> mu =15
>>>> const=1/(sigma*sqrt(2*!pi))
>>>> for i = 0,1024 do x[i]= array[0,*]
>>>> f= const*( EXP(-1.0*(x - mu)^2/(2*sigma^2)))
>
>>>> z = convol(array,array2,/center)
>>>> z = z*2
>>>> print,f
>>>> end
>
>>>> here is the message I get:% Out of range subscript encountered: X.
>>>> % Execution halted at: CONV1 29
>>>> /Users/Dave/Desktop/conv1.pro
>>>> % $MAIN$
>
>>>> I don't see why this doesn't work? I am very frustrated
>
>>> The problem is on this line:
>
>>> for I = 0,1024 do x[I]= array[0,*]
>
>>> 0 to 1024 is 1025 numbers. (Count them if you
>>> don't believe me.) But X is only big enough to
>>> hold 1024 numbers, so you are, uh, going out
>>> of its subscript range, as the error message
>>> suggests.
>
>>> But this line of code is completely unnecessary.
>>> Simply typing this is enough:
>
>>> x = Reform(array[0,*])
>
>>> Cheers,
>
>>> David
>
>>> --
>>> David Fanning, Ph.D.
>>> Fanning Software Consulting, Inc.
>>> Coyote's Guide to IDL Programming:http://www.dfanning.com/
>>> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
>
>> Thank you for your help! This did indeed solve my problem.
>
>> I have developed a new problem in my convolution. It seesm i need to
>> convolve with a kernel.
>> I can only convolve two arrays and do not seem to be able to
>> incorporate the gaussian kernel I need into the convolution.
>
>> Is this a three way convolution? I do not know how to do this. I am
>> trying to convolve 2 datasets with a kernel.
>> I have tried the code below:
>
>> pro conv_nokern
>
>> Openr, lun, 'model.dat', /Get_Lun
>
>> Point_Lun, lun, 0
>> ReadF, lun, adim, bdim, num_columns
>> spec = fltarr(2, 1024)
>> readf,lun,spec
>> a = spec(0,*)
>> b = spec(1,*)
>> Free_Lun, lun
>> window,2,xsize=500,ysize=500,retain=2
>> plot,a,b,yrange=[0,1],xrange=[4265,4200]
>
>> openr,lun,'aataunorm.dat',/get_lun
>> Point_Lun, lun, 1
>> ReadF, lun, cdim, ddim, num_columns
>> data = fltarr(2, 1024)
>> readf,lun,data
>> c = data(0,*)
>> d = data(1,*)
>> window,4,xsize=500,ysize=500,retain=2
>> plot,c,d,yrange=[0,1],xrange=[4265,4200]
>> print,data
>
>> fconv=convol(b,d,/edge_truncate)
>> ;define convoution function
>> print,fconv
>
>> openw,1,'data.dat'
>> printf,1,a,fconv
>
>> fconv2=fconv/92.4259
>> window,6,xsize=500,ysize=500,retain=2
>> plot,a,fconv2,yrange=[0,1],xrange=[4265,4200]
>> close,1
>
>> end
>
> Say a and b are your vectors containing data, and you want to convolve
> with a gaussian of width sigma. Add the following code:
>
> sigma=31. ; whatever you want- make it odd though
>
> ker=findgen(ker)-(ker-1)/2.
> ker=exp(-ker^2/(2.*sigma^2)) ; turn it into a gaussian
> ker/=total(ker) ; and normalize it
>
> ;convolve a and b with ker
>
> aconvol=convol(a,ker,/edge_truncate)
> bconvol=convol(b,ker,/edge_truncate)
>
> hope this helps
>
> chris
woops, sorry. I made a typo in defining ker. Replace
sigma=31. ; whatever you want- make it odd though
ker=findgen(ker)-(ker-1)/2.
with this:
sigma=30 ; width of gaussian, even or odd
n=8*sigma+1. ;gaussian goes to +/- 4 sigma
ker=findgen(n)-(n-1)/2.
The rest should work. By the way, if you want to check what happened
to your spectrum afterwards, try this
plot,a
oplot,aconvol,linestyle=2
|
|
|
|
| Re: convolution [message #60325 is a reply to message #26609] |
Tue, 20 May 2008 03:49 |
Chris[5]
Messages: 16
Registered: May 2008
|
Junior Member |
|
|
On May 19, 12:32 pm, sarah <sarahwiddeco...@yahoo.com> wrote:
> On May 12, 9:31 pm, David Fanning <n...@dfanning.com> wrote:
>
>
>
>> sarah writes:
>>> x=make_array(1024)
>>> sigma=15
>>> mu =15
>>> const=1/(sigma*sqrt(2*!pi))
>>> for i = 0,1024 do x[i]= array[0,*]
>>> f= const*( EXP(-1.0*(x - mu)^2/(2*sigma^2)))
>
>>> z = convol(array,array2,/center)
>>> z = z*2
>>> print,f
>>> end
>
>>> here is the message I get:% Out of range subscript encountered: X.
>>> % Execution halted at: CONV1 29
>>> /Users/Dave/Desktop/conv1.pro
>>> % $MAIN$
>
>>> I don't see why this doesn't work? I am very frustrated
>
>> The problem is on this line:
>
>> for I = 0,1024 do x[I]= array[0,*]
>
>> 0 to 1024 is 1025 numbers. (Count them if you
>> don't believe me.) But X is only big enough to
>> hold 1024 numbers, so you are, uh, going out
>> of its subscript range, as the error message
>> suggests.
>
>> But this line of code is completely unnecessary.
>> Simply typing this is enough:
>
>> x = Reform(array[0,*])
>
>> Cheers,
>
>> David
>
>> --
>> David Fanning, Ph.D.
>> Fanning Software Consulting, Inc.
>> Coyote's Guide to IDL Programming:http://www.dfanning.com/
>> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
>
> Thank you for your help! This did indeed solve my problem.
>
> I have developed a new problem in my convolution. It seesm i need to
> convolve with a kernel.
> I can only convolve two arrays and do not seem to be able to
> incorporate the gaussian kernel I need into the convolution.
>
> Is this a three way convolution? I do not know how to do this. I am
> trying to convolve 2 datasets with a kernel.
> I have tried the code below:
>
> pro conv_nokern
>
> Openr, lun, 'model.dat', /Get_Lun
>
> Point_Lun, lun, 0
> ReadF, lun, adim, bdim, num_columns
> spec = fltarr(2, 1024)
> readf,lun,spec
> a = spec(0,*)
> b = spec(1,*)
> Free_Lun, lun
> window,2,xsize=500,ysize=500,retain=2
> plot,a,b,yrange=[0,1],xrange=[4265,4200]
>
> openr,lun,'aataunorm.dat',/get_lun
> Point_Lun, lun, 1
> ReadF, lun, cdim, ddim, num_columns
> data = fltarr(2, 1024)
> readf,lun,data
> c = data(0,*)
> d = data(1,*)
> window,4,xsize=500,ysize=500,retain=2
> plot,c,d,yrange=[0,1],xrange=[4265,4200]
> print,data
>
> fconv=convol(b,d,/edge_truncate)
> ;define convoution function
> print,fconv
>
> openw,1,'data.dat'
> printf,1,a,fconv
>
> fconv2=fconv/92.4259
> window,6,xsize=500,ysize=500,retain=2
> plot,a,fconv2,yrange=[0,1],xrange=[4265,4200]
> close,1
>
> end
Say a and b are your vectors containing data, and you want to convolve
with a gaussian of width sigma. Add the following code:
sigma=31. ; whatever you want- make it odd though
ker=findgen(ker)-(ker-1)/2.
ker=exp(-ker^2/(2.*sigma^2)) ; turn it into a gaussian
ker/=total(ker) ; and normalize it
;convolve a and b with ker
aconvol=convol(a,ker,/edge_truncate)
bconvol=convol(b,ker,/edge_truncate)
hope this helps
chris
|
|
|
|
| Re: convolution [message #60336 is a reply to message #26609] |
Mon, 19 May 2008 15:32 |
sarah[1]
Messages: 3
Registered: April 2008
|
Junior Member |
|
|
On May 12, 9:31 pm, David Fanning <n...@dfanning.com> wrote:
> sarah writes:
>> x=make_array(1024)
>> sigma=15
>> mu =15
>> const=1/(sigma*sqrt(2*!pi))
>> for i = 0,1024 do x[i]= array[0,*]
>> f= const*( EXP(-1.0*(x - mu)^2/(2*sigma^2)))
>
>> z = convol(array,array2,/center)
>> z = z*2
>> print,f
>> end
>
>> here is the message I get:% Out of range subscript encountered: X.
>> % Execution halted at: CONV1 29
>> /Users/Dave/Desktop/conv1.pro
>> % $MAIN$
>
>> I don't see why this doesn't work? I am very frustrated
>
> The problem is on this line:
>
> for I = 0,1024 do x[I]= array[0,*]
>
> 0 to 1024 is 1025 numbers. (Count them if you
> don't believe me.) But X is only big enough to
> hold 1024 numbers, so you are, uh, going out
> of its subscript range, as the error message
> suggests.
>
> But this line of code is completely unnecessary.
> Simply typing this is enough:
>
> x = Reform(array[0,*])
>
> Cheers,
>
> David
>
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming:http://www.dfanning.com/
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
Thank you for your help! This did indeed solve my problem.
I have developed a new problem in my convolution. It seesm i need to
convolve with a kernel.
I can only convolve two arrays and do not seem to be able to
incorporate the gaussian kernel I need into the convolution.
Is this a three way convolution? I do not know how to do this. I am
trying to convolve 2 datasets with a kernel.
I have tried the code below:
pro conv_nokern
Openr, lun, 'model.dat', /Get_Lun
Point_Lun, lun, 0
ReadF, lun, adim, bdim, num_columns
spec = fltarr(2, 1024)
readf,lun,spec
a = spec(0,*)
b = spec(1,*)
Free_Lun, lun
window,2,xsize=500,ysize=500,retain=2
plot,a,b,yrange=[0,1],xrange=[4265,4200]
openr,lun,'aataunorm.dat',/get_lun
Point_Lun, lun, 1
ReadF, lun, cdim, ddim, num_columns
data = fltarr(2, 1024)
readf,lun,data
c = data(0,*)
d = data(1,*)
window,4,xsize=500,ysize=500,retain=2
plot,c,d,yrange=[0,1],xrange=[4265,4200]
print,data
fconv=convol(b,d,/edge_truncate)
;define convoution function
print,fconv
openw,1,'data.dat'
printf,1,a,fconv
fconv2=fconv/92.4259
window,6,xsize=500,ysize=500,retain=2
plot,a,fconv2,yrange=[0,1],xrange=[4265,4200]
close,1
end
|
|
|
|
| Re: convolution [message #60415 is a reply to message #26609] |
Mon, 12 May 2008 19:31 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
sarah writes:
> x=make_array(1024)
> sigma=15
> mu =15
> const=1/(sigma*sqrt(2*!pi))
> for i = 0,1024 do x[i]= array[0,*]
> f= const*( EXP(-1.0*(x - mu)^2/(2*sigma^2)))
>
> z = convol(array,array2,/center)
> z = z*2
> print,f
> end
>
> here is the message I get:% Out of range subscript encountered: X.
> % Execution halted at: CONV1 29
> /Users/Dave/Desktop/conv1.pro
> % $MAIN$
>
> I don't see why this doesn't work? I am very frustrated
The problem is on this line:
for I = 0,1024 do x[I]= array[0,*]
0 to 1024 is 1025 numbers. (Count them if you
don't believe me.) But X is only big enough to
hold 1024 numbers, so you are, uh, going out
of its subscript range, as the error message
suggests.
But this line of code is completely unnecessary.
Simply typing this is enough:
x = Reform(array[0,*])
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: convolution [message #60416 is a reply to message #60053] |
Mon, 12 May 2008 18:56 |
sarah[1]
Messages: 3
Registered: April 2008
|
Junior Member |
|
|
On Apr 25, 8:29 am, Vince Hradil <hrad...@yahoo.com> wrote:
> On Apr 24, 2:17 pm, sarah <sarahwiddeco...@yahoo.com> wrote:
>
>> Hi, i am IDL beginner and i do not have much support. I am trying to
>> convolve 2 spectral datasets using the convol function and I cannot
>> get it to work. I have tried everything! Does anyone know the best way
>> to do this? maybe convol is the wrong thing to use I am very lost.
>
>> Thanks
>
> What _exactly_ have you tried? What was the result? i.e. explain "I
> cannot get it to work.". Code snippets always help.
Sorry it has taken me so long to reply. (I have been sick). I have
tried writing my own convolution function and using CONVOL. It is a
very short code, here it is:
pro conv1
Openr, lun, 'model.dat', /Get_Lun
Point_Lun, lun, 0
ReadF, lun, adim, bdim, num_columns
array = fltarr(2, 1024)
readf,lun,array
a = array(0,*)
b = array(1,*)
Free_Lun, lun
Openr, lun, 'data.dat', /Get_Lun
Point_Lun, lun, 1
ReadF, lun, cdim, ddim, num_columns
array2 = fltarr(2, 1024)
readf,lun,array
c = array(0,*)
d = array(1,*)
x=make_array(1024)
sigma=15
mu =15
const=1/(sigma*sqrt(2*!pi))
for i = 0,1024 do x[i]= array[0,*]
f= const*( EXP(-1.0*(x - mu)^2/(2*sigma^2)))
z = convol(array,array2,/center)
z = z*2
print,f
end
here is the message I get:% Out of range subscript encountered: X.
% Execution halted at: CONV1 29
/Users/Dave/Desktop/conv1.pro
% $MAIN$
I don't see why this doesn't work? I am very frustrated
Sarah
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
