Re: Q: Help processing large arrays?
- Subject: [mg2571] Re: Q: Help processing large arrays?
- From: wagner at bullwinkle.cs.Colorado.EDU (Dave Wagner)
- Date: Tue, 21 Nov 1995 09:25:20 -0500
- Approved: usenet@wri.com
- Distribution: local
- Newsgroups: wri.mathgroup
- Organization: University of Colorado, Boulder
In article <DIBxnD.D4G at wri.com>,
robert fuentes <robert at mps.lfwc.lockheed.com> wrote:
> In this case, I have a 381 by 289 image that I need to run a
>Sobel operator over (i.e., convolve with a 3 by 3 Sobel operator). This
>involves taking a 3 by 3 array and moving it over the image pixel by
>pixel, applying the multipliers of the Sobel operator to the
>corresponding 3 by 3 area of the image and replacing the value of the
>image with the sum of the operation. Here is my code that I use, it seems
>to be very slow (680+ seconds on a P90 with 32Mb ram):
>
>RPConvolve[image_,rows_,cols_,kernel_,ki_,kj_] :=
> Module[
> {iadjust,jadjust,mistart,mistop,mjstart,mjstop,
> divider,nuimage,i,j,maskarea,istart,jstart,m,n},
>iadjust = (ki-1)/2;
>jadjust = (kj-1)/2;
>mistart = iadjust+1;
>mistop = rows - iadjust;
>mjstart=jadjust+1;
>mjstop = cols - jadjust;
>nuimage = image;
>For[i=mistart,i<=mistop,i++,
> For[j=mjstart,j<=mjstop,j++,
> maskarea = Table[image[[m,n]],
> {m,i-iadjust,i+iadjust},
> {n,j-jadjust,j+jadjust}];
> nuimage[[i,j]] = Apply[Plus,Flatten[maskarea*kernel]];
> ];
> ];
>nuimage
>]
>
>thanks for any help,
>Robert Fuentes
>robert at mps.lfwc.lockheed.com
To start with, a simple change you can make that will speed things up
by a small constant factor is to replace the pair of For loops with
a single Do loop:
Do[ (* body of code *), {i, mistart, mistop}, {j, mjstart, mjstop}]
This isn't changing your programming paradigm but Do's are more efficient
than For's.
Now if you want to make a radical paradigm shift, you can try this
technique, which is due to Richard Gaylord. It's sort of like
"bringing the mountain to Mohammed" rather than vice-versa:
Instead of moving the convolution kernel along each pixel of the image,
move the entire image over the convolution kernel! (You will need
gobs of memory to do this.)
The basic idea is, suppose your kernel is
a b c
d e f
g h i
What you want to do first is multiply the entire image by e, which can be
done in Mathematica using the simple statement e*image.
Next, you want to rotate the entire image down one row and multiply
that by b. Since the image is a list of rows, you rotate the
entire image down using a simple RotateRight[image, 1].
Then, rotate up one row (RotateLeft[image, 1]) and multiply by h.
Rotate left one column (Map[RotateLeft[#,1]&, image]) and multiply
by d. Similarly for the other elements of the kernel.
When you're done, you'll have 9 partially convolved images (hence the
need for memory). Add them up, and you have your entire convolved image.
Of course, you might want to do something a bit more sophisticated
than simple rotation in order to account for the edge effects.
By the way, note that there are no explicit loops, and nowhere does the
size of the image appear in the code. Functional programming at its
finest! (The size of the convolution kernel, however, is "hard-wired" in.)
Richard uses this technique to simulate physical systems. See his
book on the subject.
Dave Wagner
Principia Consulting
(303) 786-8371
dbwagner at princon.com
http://www.princon.com/princon