Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
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
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 >
On Apr 24, 9: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 did you try fft?or fftw? Cheers
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.
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
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.")
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
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
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