; Copyright (c) 1998-2016 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		AX_INTERP_SQUARE
;
;LAST CHANGED: ----------------------------------- 	26 Aug 2016 (bw)
;
; PURPOSE:
;	This function takes 2D data recorded with non-square pixels
; and interpolates the sampling of the larger pixel axis so that the pixels are square
;
; CATEGORY:
;	Image processing
;
; CALLING SEQUENCE:
;
;	Result = AX_INTERP_SQUARE(TMP)
;
; INPUTS:
; TMP	- axis 2D data structure
;
; KEYWORDS: none
;
; OUTPUTS:
;	returns modified axis 2D data structure with square pixels
;
; COMMON BLOCKS:
;	AXIS_COM	standard set of common blocks
;   SDF.COM     common for sdf
;	NEXUS.COM	common for nexus
;
; MODIFICATION HISTORY:
; (11-Apr-16 aph) first written; taken from read_sdf
; (21-Apr-16 aph) correct error in com file names
; (26-Aug-16 bw ) properly handle images that were scanned in negative x and/or y directions
;-

function ax_interp_square, tmp
@axis_com
@sdf_com
@nexus_com

; --- check if 2d structure
if tmp.t NE '2d' then begin
	print, 'ax_interp_square: not a 2D data structure'
	return, tmp
endif
; load variables
if n_elements(pix_diff) LE 0 then pix_diff = 0.005
x = tmp.x
y = tmp.y
d = tmp.d
; help, tmp, /struct

; --- check if rectangular pixels - if so interpolate to square
if NOT keyword_set(no_interpolate) then begin
	dx = abs(x(1) - x(0))
	dy = abs(y(1) - y(0))
	if abs(dx - dy) GT pix_diff then begin
		if dx GT dy then begin
			axis_log, ' Rectangular pixels: x interpolated'
			nx = fix((x(n_elements(x)-1) - x(0))/dy) + 1
			x = interpol(x,nx)
			tmp.xl = 'interp-square ' + tmp.xl
		endif else begin
			axis_log, ' Rectangular pixels: y interpolated'
			ny = fix((y(n_elements(y)-1) - y(0))/dx) + 1
			y = interpol(y,ny)
			tmp.yl = 'interp-square ' + tmp.yl
		endelse
		nx = n_elements(x) & ny = n_elements(y)
		td = congrid(d,nx,ny, /interp, cubic=-0.5)
		d = td
	endif
endif
; --- check if scan direction is not positive - if so, then flip data array
if NOT keyword_set(no_interpolate) then begin
  dx = x(1) - x(0)
  dy = y(1) - y(0)
  if dx LT 0 then begin
    d = reverse(d, 1)
  endif
  if dy LT 0 then begin
    d = reverse(d, 2)
    if not finite(d(0)) then d(0) = 0    ;There is a bug where display gets strange if the first pixel is NaN (I also saw IDL crash hard)
  endif
endif

tmp = {t:'2d', x:x, y:y, d:d, e:tmp.e, xl:tmp.xl, yl:tmp.yl, dl:tmp.dl}


return, tmp
end