Maximum Value Array Resampling

QUESTION: I'm trying to figure out an IDL-efficient way to resample a series of images. I know how to do this is a ruinously laborious fashion using loops but I know there's any easier way.

Consider the following 4x4 array as an example:

   x = [$
   [0,3,4,5],$
   [1,2,7,0],$
   [3,2,9,0],$
   [7,0,5,6]]

I would like to resample this to a 2x2 array named Y. Each element would contain the maximum value of the corresponding four pixels. Y would then look like this:

 [3,7]
 [7,9]

So element [0,0] in Y is Max(x[0:1,0:1]) etc. As I understand it, Rebin/Congrid won't do this. Each image is about 5000x2000 and there are several hundred to process.

ANSWER: You can find this entire discussion in the IDL newsgroup archives. Search for the topic "Maximum value array resampling" on 5 August 2005. This is a highly edited version of that discussion, with answers provided by Dick French and JD Smith.

First, Dick French.

Try this. The basic idea is to create a vector of indices for each cell. Once you know the indices of each of the four entries in each cell, you can do a direct compare. It took 7-8 seconds on my Powerbook G4 for a 4000 x 5000 image, compared to 40-45 seconds using loops.

   PRO Example

   nx=4000L
   ny=5000L

   ; put in random values (this is the slow step!)

   X=rebin(fix(1000*(randomu(seed,nx*ny))),nx,ny)
   nx2=nx/2L
   ny2=ny/2L

   ; get indices of upper left element of each 2x2 cell

   l=rebin(nx*2#lindgen(ny2),nx2,ny2)+rebin(lindgen(nx2)#2,nx2,ny2)


   ; Index method: compare with indices of ul,ur, ll, lr of each 2x2 cell

   print,'Start....'
   T10=systime(1)
   xfinal=x[l]>x[l+1]>x[l+nx]>x[l+nx+1]
   print,'Time=',systime(1)-t10

   ; Loop method:

   t20=systime(1)
   y=intarr(nx2,ny2)
     for i=0,nx2-1 do x[i,*]=x[2*i,*] > x[2*i+1,*]
     for i=0,ny2-1 do y[*,i]=x[0:nx2-1,2*i] > x[0:nx2-1,2*i+1]
   print,'Time=',systime(1)-t20

   print,max(abs(xfinal-y)) ; confirm that we get same result

   END

JD Smith then suggested an alternative (and typically compact) method to find Y, and suggests a small change to Dick's method.

For arbitrary images with both dimensions even:
d = size(x,/DIMENSIONS) & nx = d[0]/2 & ny = d[1]/2
y = transpose(max(reform(transpose(reform(x,2,nx,2*ny),[0,2,1]), 4,ny,nx),DIMENSION=1))
How does it work? It juggles dimensions so that the indices of all the 2x2 sub-arrays are next to each other in memory, and then uses Max(DIMENSION=1) to collapse over them. The inner call to REFORM puts them adjacent to each other, but in the wrong dimension, then TRANSPOSE makes them adjacent in the fast-changing dimension. The rest is straightforward.

There may be a quicker way with only one call to TRANSPOSE, but I couldn't find it (anyone?). Also, if you don't care to keep X, throw a couple of /OVERWRITE keywords for both REFORM statements, to save some memory and time.

By the way, since for processing many images of the same size you can pre-compute the indices, you might consider a small change to the modified loop method:
   ;; Pre-compute the indices:
   d=size(x,/DIMENSIONS) & nx=d[0] & ny=d[1]
   nx2=nx/2 & ny2=ny/2

   inds1=rebin(lindgen(nx2)*2L,nx2,ny2,/SAMPLE)+ $
         rebin(transpose(lindgen(ny2)*2L*nx),nx2,ny2,/SAMPLE)
   inds2=inds1+1L
   inds3=inds1+nx
   inds4=inds1+nx+1L

   ;; Form the sub-sampled image (for each image)
   xmax=x[inds1]>x[inds2]>x[inds3]>x[inds4]
This brings the total processing time per 5000x4000 image to under 0.5s on my not-so-fast Linux box (and doesn't have any loops ;).

And then finally, not being able to let it go, JD comes up with a general purpose routine that can be used to find the maximum or minimum resampled value and leaves the rest of us with a further challenge.

Here's a version which can do MAX or MIN, and any box size (not just 2x2). Keep in mind that for arrays which are not multiples of the box size, the edges will be lost. Another interesting question is how to do a sliding box MAX/MIN efficiently, ala MEDIAN.

   ;; BOX_MAX: Compute local maximum (or minimum. with /MIN) box
   ;; downsampling, with box size boxx, boxy (default 2,2).  Pre-computed
   ;; INDS may be passed.
   ;;   JD Smith (c) 2005.
   function box_max,array,boxx,boxy,INDS=inds, MIN=min
     if n_elements(boxx) eq 0 then boxx=2
     if n_elements(boxy) eq 0 then boxy=2
     min=keyword_set(min) 
     d=size(array,/DIMENSIONS) 
     nx=d[0] & ny=d[1]
  
     if n_elements(inds) eq 0 then begin 
        nx_out=nx/boxx & ny_out=ny/boxy
        inds=rebin(lindgen(nx_out)*boxx,nx_out,ny_out,/SAMPLE)+ $
             rebin(transpose(lindgen(ny_out)*boxy*nx),nx_out,ny_out,/SAMPLE)
     endif 
  
     ret=array[inds]
     for i=0L,boxx-1L do begin 
        for j=0L,boxy-1L do begin 
           if i eq 0 && j eq 0 then continue
           if min then ret <= array[inds+i+j*nx] else ret >= array[inds+i+j*nx]
        endfor
     endfor
  
     return,ret
   end