;;;; -*- Syntax: IDL; Mode: Lisp -*-

pro utilities
;+
;.NAME:
;	utilities.pro
;
;.PURPOSE:
;	Define a bunch of handy IDL utilities.
;
;.COMMENT:
;	Compile with ".run utilities" from IDL V2.
;
;.AUTHOR:
;	Hans-Martin Adorf
;	Space Telescope - European Coordinating Facility
;	European Southern Observatory
;	D-8046 Garching bei Muenchen
;	Karl-Schwarzschild-Str. 2
;	F.R. Germany
;	Internet: adorf@eso.org		SPAN:  ESO::ADORF or adorf@eso.span
;
;.MODIFICATION HISTORY:
;	Jul 1990	HMA	Initial implementation
;	Aug 1990	HMA	Further improvements
;	Sep 1991	HMA	Some comments provided
;       Oct 1991	HMA	Added new features
;       Feb 1992        HMA     Added new utilities
;       jan 1994        M.B     Truncated version for IDL v 3.1.1 <==> LASCO
;
;.TO DO
;	Introduce global variables to control, e.g. the verbosity of
;	graphics and image display of intermediate results.
;-
  tools					;require tools.pro
  return
end

;;;----------------------------------------------------------------------------
;;; Math functions
;;; It's a shame that these have to be implemented by the user.
;;;----------------------------------------------------------------------------
function  floor1, n, base
  if (n_params(0) lt 2) then base = 1
  return, (fix(base)*fix(n/base))
end

function my_mod, n, base
;; IDL refused to accept MOD here.
   return, (n - floor1(n, base))
end

function sgn, x
  ;; Returns the sign of an array.
  result = 0*x
  pos_supp = where(x gt 0, pos_cnt)
  if (pos_cnt gt 0) then result(pos_supp) = +1
  neg_supp = where(x lt 0, neg_cnt)
  if (neg_cnt gt 0) then result(neg_supp) = -1
  return, result
end

function my_atan, x, y
  ;; We can't use the built-in single-argumentatan, since it looses the 
  ;; information about quadrants etc.
  ;; see Steele, p. 306
  result = 0.*y
  supp_y = where(y gt 0, cnt)
  if (cnt gt 0) then begin
    supp = where(x(supp_y) gt 0, cnt)	;Quadrant I
    if (cnt gt 0) then result(supp) = atan(y(supp)/x(supp))
    supp = where(x(supp_y) eq 0, cnt)	;Positive y-axis
    if (cnt gt 0) then result(supp) = !pi/2.0
    supp = where(x(supp_y) lt 0, cnt)	;Quadrant II
    if (cnt gt 0) then result(supp) = atan(y(supp)/x(supp)) + !pi
  endif

  supp_y = where(y eq 0, cnt)
  if (cnt gt 0) then begin
    supp = where(x(supp_y) gt 0, cnt)	;Positive x-axis
    if (cnt gt 0) then result(supp) = 0.
    supp = where(x(supp_y) eq 0, cnt)	;Origin
    if (cnt gt 0) then result(supp) = 0.	;should be error
    supp = where(x(supp_y) lt 0, cnt)	;Negative x-axis
    if (cnt gt 0) then result(supp) = !pi
  endif

  supp_y = where(y lt 0, cnt)
  if (cnt gt 0) then begin
    supp = where(x(supp_y) gt 0, cnt)	;Quadrant IV
    if (cnt gt 0) then result(supp) = atan(y(supp)/x(supp))
    supp = where(x(supp_y) eq 0, cnt)	;Negative y-axis
    if (cnt gt 0) then result(supp) = -!pi/2.0
    supp = where(x(supp_y) lt 0, cnt)	;Quadrant III
    if (cnt gt 0) then result(supp) = atan(y(supp)/x(supp)) - !pi
  endif
  return, result
end

function arg, x
  ;; Returns the phase of an array.
  z = complex(x)
  return, my_atan(re(z), im(z))
end

function imaginary_unit
  return, sqrt(complex(-1))
end

function cis, phi
  ;; Returns exp(i phi) = cos(phi) + i*sin(phi)
  i = imaginary_unit()
  return, cos(phi) + i*sin(phi)
end

function symmetrize, array
  ;; works only for odd frames
  return, (array + rotate(array, 2))/2.0
end

function antisymmetrize, array
  return, (array - rotate(array, 2))/2.0
end
  
;;;----------------------------------------------------------------------------
;;; Predicates
;;; should be vectorized, too.
;;;----------------------------------------------------------------------------
function evenp, n
  return, (my_mod(n, 2) eq 0)
end

function oddp, n
  return, (my_mod(n, 2) eq 1)
end

function even_or_odd, n
;; The (-1)^n function
  if (my_mod(n, 2) eq 0) then result = 1 else result = -1
  return, result
end

;;;----------------------------------------------------------------------------
;;; Auxiliary functions
;;;----------------------------------------------------------------------------
;; This function should be replaced by the system function N_ELEMENTS
function no_of_pixels, array
;; Return the no of pixels of a 1D or 2D data array.
  return, n_elements(array)       ; modified by M.B 24/01/94 :LAS
end

function origin, frame
  ;; Defines the origin of even- and odd-sized frames.
  n = size(frame)
  result = intarr(2)
  result(0) = n(1)/2
  if (n(0) eq 2) then result(1) = n(2)/2 else result(1) = 0
  return, result
end

function centre, image
  ;; Shift the origin of the image from the bottom left corner to the centre.
  o = origin(image)
  s = size(image)
  if (s(0) eq 1) then $
    return, shift(image, o(0)) $
  else $
    return, shift(image, o(0), o(1))
end

function decentre, image
  ;; Shift the origin of the image from the centre to the bottom left corner.
  o = origin(image)
  s = size(image)
  if (s(0) eq 1) then $
    return, shift(image, -o(0)) $
  else $
    return, shift(image, -o(0), -o(1))
end

function unit_pulse, nx, ny
  ;; Puts a 1 at the origin of an artificial frame. (centered Dirac)
  if (n_params(0) lt 2) then ny = 1
  x = fltarr(nx, ny)
  o = origin(x)
  x(o(0), o(1)) = 1.0
  return, x
end

;;;----------------------------------------------------------------------------
;;; Frame insertion and extraction
;;;----------------------------------------------------------------------------
function insert, child, parent
  ;; Insert child array into parent array (centred)
  ;; This f. is not tested for odd frames and does presumably not work as is.
  nc = size(child)
  if (nc(0) eq 1) then begin
      oc = origin(child) & ocx = oc(0)
      op = origin(parent) & opx = op(0)
      parent(opx-ocx:opx-ocx+nc(1)-1) = child 
    endif else begin
      oc = origin(child) & ocx = oc(0) & ocy = oc(1)
      op = origin(parent) & opx = op(0) & opy = op(1)
      parent(opx-ocx:opx-ocx+nc(1)-1,opy-ocy:opy-ocy+nc(2)-1) = child
  endelse
  return, parent
end

function extract, image, nx, ny
  ;; Extract a subframe of dimensions (nx ny) centred on the origin
  if (n_params(0) lt 3) then ny = 1
  n = size(image)
  o = origin(image)
  nxm = nx/2
  nym = ny/2
  if oddp(nx) then nxp = nxm else nxp = nxm-1
  if oddp(ny) then nyp = nym else nyp = nym-1
  if (n(0) eq 1) then begin
    result = image(o(0)-nxm:o(0)+nxp)
  endif else begin
    result = image(o(0)-nxm:o(0)+nxp,o(1)-nym:o(0)+nyp)
  endelse
  return, result
end

function normalize, data
  ;; Put sum over all values to one.
  return, data/total(data)
end

function upsample, array, step
  ;; Insert zeros in between data values.
  ;; Attention: the origin is not properly mapped.
  if (n_params(0) lt 2) then step = 2
  n = size(array)
  case n(0) of
  1: begin
       result = complexarr(step*n(1))
       for i = 0, n(1)-1 do result(step*i) = array(i)
     end
  endcase
  return, result
end

function downsample, array, step
  if (n_params(0) lt 2) then step = 2
  n = size(array)
  case n(0) of
  1: begin
       result = complexarr(n(1)/step)
       for i = 0, n(1)/step-1 do result(i) = array(step*i)
     end
  endcase
  return, result
end

function clip, array
  ; result = array
  ; supp = where(array lt 0, cnt)
  ; if (cnt gt 0) then result(supp) = 0.
  ; return, result
  return, array > 0.
end

;;;----------------------------------------------------------------------------
;;; Fourier transformation, convolution etc.
;;;----------------------------------------------------------------------------

function dft, data
  ;; Discrete Fourier transform
  ;; The IDL FFT normalizes during the forward transformation.
  ;; return, fft(data, -1)
  return, centre(fft(decentre(data), -1))
end

function dfti, ft
  ;; Inverse discrete Fourier transform
  ;; return, fft(ft, +1)
  return, centre(fft(decentre(ft), +1))
end

function convolve, x, y
  ;; Convolve two (real) arrays.
  ;; We have shift the result, since the IDL origin is in the bottom left 
  ;; corner.
  ;; For efficiency reasons, an odd dimensioned frame is embedded in a 
  ;; larger even frame
  n = size(x)  
  nx = n(1)/2
  ny = n(2)/2
  if (oddp(n(1)) and oddp(n(2))) then begin
    odd = 1
    x = insert(x, fltarr(n(1)+1, n(2)+1))
    y = insert(y, fltarr(n(1)+1, n(2)+1))
  endif else begin
    odd = 0
  endelse
  z = no_of_pixels(x)*float(dfti(dft(x)*dft(y)))
  ;; x = 0 & y = 0				;save space
  if odd then z = z(0:n(1)-1,0:n(2)-1)	;remove border
  return, z
end

function correlate, x, y
  ;; Correlate two (real) arrays.
  ;; We have shift the result, because the origin is in the 
  ;; bottom left corner.
  ;; Apart from taking the conjugate of y, this function is identical to 
  ;; the convolution
  return, no_of_pixels(x)*float(dfti(dft(x)*conj(dft(y))))
end
;; x = [[1,2,3],[4,5,6],[7,8,9]]
;; x = [[1,2,3,4], [5,6,7,8], [9,10,11,12] ,[13,14,15,16]]

function autocorrelate, x
  ft = dft(x)
  return, no_of_pixels(x)*float(dfti(ft*conj(ft)))
end

function power_spectrum, x
  ;; The squared one
  ft = dft(x)
  return, float(ft*conj(ft))
end

function log_power_spectrum, x
  return, alog(power_spectrum(x))
end

;;;----------------------------------------------------------------------------
;;; Image manipulation
;;;----------------------------------------------------------------------------
function zero_nyquist, ft
  ;; Zero the Nyquist-frequencies of an even image
  ;; Assume the origin is centred and the Nyquist-frequencies are 
  ;; at row and column 0
  n = size(ft)
  if evenp(n(1)) then ft(0, *) = 0
  if ((n(0) eq 2) and evenp(n(2))) then ft(*, 0) = 0
  return, ft
end


function complement_nyquist, fx
  ;; Complements the missing Nyquist frequency components
  ;; required for upsampling/interpolating an even frame.
  ; fx = dft(x)				;Fourier transform of a frame
  n = size(fx) & nx1 = n(1) & nx2 = n(2)
  ny1 = nx1 + 1 & ny2 = nx2 + 1
  fy = insert(fx, complexarr(ny1, ny2))
  one = intarr(nx1, nx2)+1
  mask = insert(one, intarr(ny1, ny2))
  supp = where(mask eq 0)
  dummy = rotate(fy, 2)
  fy(supp) = conj(dummy(supp))		;complement values
  corner_value = fy(0, 0)/4		;factor valid only for 2D arrays
  fy(0, 0) = corner_value		;  according to Antoine Llebaria, LAS
  fy(0, ny2-1) = corner_value
  fy(ny1-1, 0) = corner_value
  fy(ny1-1, ny2-1) = corner_value
  ; y = dfti(fy)
  return, fy
end

function resample, image, nx, ny
  ;; Resample the given image to (nx, ny) dimensions using 
  ;; aliased sinc-interpolation
  case n_params() of
    2: ny = 1
    3:
  endcase
  factor = float(no_of_pixels(image))/nx/ny ;normalizing factor
  print, "... FFT given (source) image ..."
  old_f = dft(image)
  ; old_f = zero_nyquist(old_f)		; suppress Nyquist frequency comp.
  old_f = complement_nyquist(old_f)	; complement Nyquist frequency comp.
  print, "... insert source into destination ..."
  ;; IDL normalizes during the forward transform.
  new_f = factor*insert(old_f, complexarr(nx,ny))
  ; old_f = 0				          ; save memory space
  print, "... inverse fft ..."
  new = dfti(new_f)
  ; new_f = 0				          ; save memory space
  return, float(new)
end

function replicate_image, image
  ;; Replicate an image into a frame with twice the linear size
  ;; by replicating each pixel.
  s = size(image)
  print, s
  image_repl = fltarr(2*s(1), 2*s(2))
  for i = 0, s(1)-1 do begin
    ir = 2*i
    for j = 0, s(2)-1 do begin
      jr = 2*j
      v = image(i,j)
      image_repl(ir, jr) = v
      image_repl(ir+1, jr) = v
      image_repl(ir, jr+1) = v
      image_repl(ir+1, jr+1) = v
    endfor
  endfor
  return, image_repl
end

;;;----------------------------------------------------------------------------
;;; Utilities for creating masks
;;;----------------------------------------------------------------------------
function band_mask, npix, width
  ;; Create a band mask of size npix*npix with a horizontal band of width.
  mask = fltarr(npix, npix)
  o = origin(mask) & ox = o(0)
  w2 = width/2
  mask(ox-w2:ox+w2, *) = 1.0
  return, mask
end
    
function circle_mask, npix, radius
  ;; Create a circular mask, 1 inside, 0 outside a circle of radius.
  nh = fix(radius+0.5)
  rs_max = (float(nh))^2
  n_max = 2*nh
  np = n_max + 1
  ns = np*np
  mask = fltarr(np, np)
  o = origin(mask) & ox = o(0) & oy = o(1)
  for i = 0, n_max do begin
;    print, "i = ", i
    inh = i - ox
    is = (float(inh))^2
    for j = 0, n_max do begin
      jnh = j - oy
      rs = is + (float(jnh))^2
      if (rs le rs_max) then mask(i,j) = 1.
    endfor
  endfor
  return, insert(mask, fltarr(npix, npix))
end

function ellipse_mask, npix, a, b, theta
  ;; modified by M.B on 25/01/94 : LAS
  ;; Create an ellipse mask, 1 inside, 0 outside the ellipe
  ;; a		= major half axis
  ;; b		= minor half axis
  ;; theta	= angle of between x-axis and major axis (degrees)
  as = (float(a))^2
  bs = (float(b))^2
  ct = cos(theta*!dpi/180.)
  st = sin(theta*!dpi/180.)
  nh = npix/2
  n_max = 2*nh
  mask = fltarr(npix, npix)
  for i = 0, n_max do begin
;    print, "i = ", i
    inh = i - nh
    x = float(inh)
    for j = 0, n_max do begin
      jnh = j - nh
      y = float(jnh)
      xp =  ct*x + st*y
      yp = -st*x + ct*y
      rs = sqrt((xp/a)*(xp/a) + (yp/b)*(yp/b))
      if (rs le 1.) then mask(i,j) = 1.
    endfor
  endfor
  return, mask
end

function gauss_mask, npix, sigma
  ;; Create a circular mask with gaussian values
  a = 1./sqrt(2.*3.14159)/sigma
  b = 0.5/sigma^2
  n_max = npix - 1
  nh = float(npix/2)-0.5
  ns = npix*npix
  mask = fltarr(npix, npix)
  for i = 0, n_max do begin
;    print, "i = ", i
    inh = i - nh
    is = (float(inh))^2
    for j = 0, n_max do begin
      jnh = j - nh
      rs = is + (float(jnh))^2
      mask(i, j) = a*exp(-b*rs)
    endfor
  endfor
  return, mask
end

function cone_mask, npix, r
  ;; Create a circular mask with cone values
  n = fix(2*r)
  mask = fltarr(n, n)
  o = origin(mask) & ox = o(0) & oy = o(1)
  i_max = n - 1
  ns = n*n
  for i = 0, i_max do begin
;    print, "i = ", i
    for j = 0, i_max do begin
      x = float(i - ox)
      y = float(j - oy)
      mask(i, j) = sqrt(x*x + y*y)
    endfor
  endfor
  return, insert(clip(-mask+r), fltarr(npix, npix))
end

function rolloff_mask, npix, radius, sigma
  ;; Create a circular mask, 1 inside, gaussian values outside a cirle 
  ;; npix	= linear number of pixels
  ;; radius	= radius of circular disk
  ;; sigma	= std. dev. of gaussian roll off
; a = 1./sqrt(2.*3.14159)/sigma
  b = 0.5/sigma^2
  nh = npix/2
  rs_max = (float(radius))^2
  n_max = npix - 1
  mask = fltarr(npix, npix)
  for i = 0, n_max do begin
;    print, "i = ", i
    inh = i - nh
    is = (float(inh))^2
    for j = 0, n_max do begin
      jnh = j - nh
      rs = is + (float(jnh))^2
      if (rs le rs_max) then $
	mask(i,j) = 1. $
      else $
        mask(i,j) = exp(-b*(sqrt(rs)-radius)^2)
    endfor
  endfor
  return, mask
end

;;;----------------------------------------------------------------------------
;;; Some functions
;;;----------------------------------------------------------------------------
function euclid, x, y
  return, sqrt(x*x+y*y)
end

function gauss, x, mu, sigma
  c = 2.50663	; sqrt(2*pi)
  return, exp(-((x-mu)/sigma)^2/2)/c/sigma
end

;;;----------------------------------------------------------------------------
;;; Statistics, e.g. for background estimation.
;;;----------------------------------------------------------------------------
function mean, array
  ;; Return the average of a 1D or 2D data array.
  npix = size(array)
  return, total(array) / no_of_pixels(array)
end

function maxpos, array
  ;; Return the positions of the absolute maxima
  n = size(array) & nx = n(1) & ny = n(2)
  m = max(array)
  w = where(array eq m)
  ix = my_mod(w, nx)
  iy = floor1(w/nx)
  print, "          max1        max2 ..."
  return, [[ix], [iy]]
end

function statistics, array
  ;; modified by M.B on 25/01/94 : LAS
  ;; Compute simple statistics of the array.
  n = no_of_pixels(array)
  fluxtot = total(array)
  min = min(array)
  max = max(array)
  mean = mean(array)
  sigma = stdev(array)
  sn = mean / sigma
  med = median(array)
  print, 'no of pixels = ', n
  print, 'total flux = ',fluxtot
  print, 'min          = ', min
  print, 'max          = ', max
  print, 'mean         = ', mean
  print, 'stdev        = ', sigma
  print, 'S/N          = ', sn
  print, 'median       = ', med
  return, [n, fluxtot, min, max, mean, sigma, sn, med]
end

;;;----------------------------------------------------------------------------
;;; Plot utilities  Reviewed & reconceptualized by M.B on 25/01/94 : LAS
;;; for VAX/VMS v 5
;;;----------------------------------------------------------------------------
pro copy_disp,image,minval,maxval
  set_plot,'ps'
  device,/portrait,bits_per_pixel=4
  device,filename='idlpost.ps'
  tvscl,bytscl(image,min=minval,max=maxval)
  device, /close
  set_plot,'x'
  cmd = 'lasera post idlpost.ps'
  spawn, cmd
end

pro send_plot,image,col,row,minval,maxval
  ln03,'idlplot',0
  if (col lt 0.) then begin $
    plot,image(*,row),linestyle=0,yrange=[minval,maxval]
  endif
  if (row lt 0.) then begin $
    plot,image(col,*),linestyle=0,yrange=[minval,maxval]
  endif
  if ((col ge 0.) and (row ge 0.)) then begin $
    plot,image(*,row),linestyle=0,yrange=[minval,maxval]
    oplot,image(col,*),linestyle=1
  endif
  ln03
  cmd = 'lasera tek4014 idlplot.dat'
  spawn, cmd
end

function mtf, ima
  ;; Modulation transfer function
  result = abs(dft(ima))
  shade_surf, result
  return, result
end