;Downsample an image.
; Definition of: downsample (1) To make a digital audio signal smaller by lowering
;  its sampling rate or sample size (bits per sample). Downsampling is done to
;  decrease the bit rate when transmitting over a limited bandwidth or to convert
;  to a more limited audio format. (2) To decrease the color depth of a digital image;
;  for example, from 24 bits to 16 bits per pixel.
;  
; Inputs
;   image -  two-dimensional array of data to be downsampled. This has a size [x,y]
;   amount - integer factor by which to downsample. For instance,
;            amount=10 means one row in 10 and one column in 10 is kept,
;            leading to a total of one pixel in 100 being kept.
;            IMPORTANT! amount must be a factor of both the number of rows
;            and number of columns in the input image. If not,
;            downsample will fail with an error in rebin. Special case: if amount is 1,
;            the input image is immediately returned unchanged.
; Returns
;   A downsampled version of the image. This will have a size [x/amount+1,y/amount+1]
function downsample,image,amount,longitude=longitude
  if amount eq 1 then return,image
  if keyword_set(longitude) then begin
    ;compute the sines and cosines of the input, then interpolate each of them
    ;and then take the atan of the results. This avoids all singlarities and branch
    ;cuts related to longitude.
    slon=sin(image*!const.dtor)
    clon=cos(image*!const.dtor)
    slon_downsample=downsample(slon,amount) ;We specifically do NOT set /longitude here
    clon_downsample=downsample(clon,amount) ;otherwise this is an infinite loop.
    image_downsample=atan(slon_downsample,clon_downsample)*!radeg
    return,image_downsample
  end else begin
    x=(size(image,/dimensions))[0]
    y=(size(image,/dimensions))[1]
    xd=x/amount
    yd=y/amount

    ;interpolate to [x+1,y+1]
    image_resize=interpolate(image,findgen(x+1)*(x-1)/x,findgen(y+1)*(y-1)/y,/grid)

    ; then rebin to [xd+1,yd+1]
    image_downsample=make_array(xd+1,yd+1,type=size(image,/type))
    image_downsample[0:xd-1,0:yd-1]=rebin(image_resize[0:x-1,0:y-1],xd,yd,/sample)
    image_downsample[xd,0:yd-1]=rebin(image_resize[x,0:y-1],1,yd,/sample)
    image_downsample[0:xd-1,yd]=rebin(image_resize[0:x-1,y],xd,1,/sample)
    image_downsample[xd,yd]=image_resize[x,y]

    return,image_downsample
  end
end