| Re: faster convol on local subsets? [message #78621 is a reply to message #78613] |
Mon, 05 December 2011 15:17   |
Andre
Messages: 4
Registered: June 2011
|
Junior Member |
|
|
On Dec 5, 11:54 am, Yngvar Larsen <larsen.yng...@gmail.com> wrote:
> On Dec 5, 1:37 am, Andre <note....@gmail.com> wrote:
>
>> Hello experts,
>
>> Maybe somebody has an easy solution for this?
>> I have a 2D array (img) and want the filter response from kernels that vary according to the image position. In a second array (loc, same dimensions as img) I have the information which kernel should be used at each pixel. My current approach is to first convolve the full image with the j-th kernel and take the response only at the positions with the current j indexed in the loc array:
>
>> for j=0, n do begin
>> kernel=kernel_store[*,*,j]
>> response_temp = convol(img, kernel, /edge_zero, /NAN)
>> index=where(loc eq j)
>> if (index[0] gt -1)then response[index]=response_temp[index]
>> endfor
>
>> I works fine, but it is relatively slow and I wonder if there is a smarter (faster) to apply only the convolutions that are really needed?
>
>> Thanks in advance for any help!
>
> Yes, it seems like IDL does not implement 2D convolution very
> efficiently. I found out that a straight forward implementation by
> zeropadding to a power-of-2 length followed by multiplication in the
> FFT domain is much faster unless the convolution kernel is very small.
> Something like this (when /EDGE_ZERO and /NORMALIZE is set, some more
> work for other EDGE_* keywords):
>
> sizeA = size(array, /dimensions)
> sizeB = size(kernel, /dimensions)
>
> dim1 = sizeA[0] + sizeB[0] - 1
> dim2 = sizeA[1] + sizeB[1] - 1
>
> s1 = 2L^ceil(alog(dim1)/alog(2))
> s2 = 2L^ceil(alog(dim2)/alog(2))
>
> A = dcomplexarr(s1, s2)
> B = dcomplexarr(s1, s2)
>
> A[0,0] = array
> B[0,0] = kernel
>
> convol = fft(fft(A)*fft(B), /inverse)*s1*s2
> convol = convol[sizeB[0]/2:sizeB[0]/2+sizeA[0]-1, $
> sizeB[1]/2:sizeB[1]/2+sizeA[1]-1]
>
> convol = double(convol)/total(abs(kernel))
>
> --
> Yngvar
Thanks for all the suggestions so far.
I tried it with the changes that Jeremy suggested but for some reason
it runs even a little bit slower than the original version.
On a 2300x2900 array the original loop runs for 322.46600s while with
REVERSEINDICES it needs 394.51800s (even when precomputing the kernel
outside the loop). My guess is that it takes more time because calling
the routine in each loop is expensive (http://ross.iasfbo.inaf.it/IDL/
Robishaw/idlfast.html).
I did not yet find time to check the implementation that Yngvar
suggested but tried http://idlastro.gsfc.nasa.gov/ftp/pro/image/convolve.pro
which also implements convolution in the Fourier domain. Still its
slower than the native IDL convolution. According to a comment in
their code IDL 8.1 has a CONVOL_FFT() which seems worth a further try
after I got the update.
Last I also tried to convolve at each position only with desired
kernel. The code looks more or less like this
m=half_kernel_size
nc= number of columns of the input
nr = number of rows of the input
for i=m, nc - m-1 do begin
for j=m, nr - m-1 do begin
patch=img[i-m:i+m, j-m:j+m]
kernel=kernel_store[*,*, (max_loc[i,j])]
temp = convol(patch, kernel])
response[i,j] = temp[m, m]
endfor
endfor
As expected the convolution runs much quicker than on the full image
but the large number of loops eats up all that speed gains and in the
end its even 409.56500s for the same array as before.
...to be continued...
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
