Re: Q: Help processing large arrays?

*To*: mathgroup at smc.vnet.net*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*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