# 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 toREFORMputs them adjacent to each other, but in the wrong dimension, thenTRANSPOSEmakes 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/OVERWRITEkeywords for bothREFORMstatements, 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

Copyright © 2006 David W. Fanning

Last Updated 11 January 2006