;+
; NAME: ircal_dewarp
; PURPOSE:
; Remove distortion from an IRCAL image
;
; NOTES:
; Right now this is -incredibly- basic and just handles the anisomorphic
; magnification term, and nothing higher order than that.
;
; INPUTS:
; image0 an image or stack of images.
; KEYWORDS:
; /NOFLUXSCALE Don't rescale image counts to preserve total flux.
; (default is to do scaling)
; pixelscale= new pixelscale to use. Default is 0.040 arcsec
; CUBIC= use cubic resampling. See docs for INTERPOLATE.
; /BYTE return a byte image, not float. Useful for dewarping
; masks
; /CROPTOP Drop the top 5 rows of the IRCAL detector. These are the
; rows full of many bad pixels and sometimes it's just
; easiest to get rid of them.
; OUTPUTS:
;
; HISTORY:
; Began 2005-11-07 15:24:51 by Marshall Perrin
; 2006-05-09 Actually got it working finally. Major code updates.
; 2006-06-09 Added /croptop option
;-
function ircal_dewarp,image0,nofluxscale=nofluxscale,pixelscale=pixelscale,$
cubic=cubic,byte=byte,croptop=croptop
@instruments.include
scalex = ircal.pixelscale
scaley = ircal.pixelscale_y
; what should be the new (regridded) pixelscale?
;if ~(keyword_set(pixelscale)) then pixelscale=scalex
if ~(keyword_set(pixelscale)) then pixelscale=0.040 ; Nyquist at J for Lick.
; TODO would this be better at twice Nyquist?
image = image0
sz = size(image)
; because pixel coordinates are the CENTERS of pixels, an image with
; N pixels is actually (N-1)*pixelscale units wide.
; (well, not counting the half-pixel border on either side)
; new image pixel counts
nsx = floor((sz[1]-1)*scalex/pixelscale)+1
nsy = floor((sz[2]-1)*scaley/pixelscale)+1
; should we just drop the top 5 rows (with tons of bad pixels?)
ncrop=5
if keyword_set(croptop) then nsy = floor((sz[2]-ncrop-1)*scaley/pixelscale)+1
; handle 3D arrays if necessary via recursion
if sz[0] eq 3 then begin
out = dblarr(nsx,nsy,sz[3])
for i=0L,sz[3]-1 do begin
statusline, "Dewarping "+strc(i+1)+"/"+strc(sz[3])
out[0,0,i] = ircal_dewarp(image[*,*,i],nofluxscale=nofluxscale,$
pixelscale=pixelscale,cubic=cubic,croptop=croptop)
endfor
return,out
endif
; now actually handle 2D arrays.
;
; indices of new pixels relative to old
ix = findgen(nsx)*pixelscale/scalex
iy = findgen(nsy)*pixelscale/scaley
; can't use interpolate with NANs, so fix if needed.
wnans = where(~finite(image),nanct)
if nanct gt 0 then begin
nanmask = finite(image)
;TODO add an option for more clever nan repair here.
image[wnans] = median(image)
ix2 = ix ; needed because bilinear stomps on its input IX and IY arrays.
iy2 = iy
newnanmask = bilinear(nanmask,ix2,iy2)
endif
; On tests with linear ramp data, the following code works very well on all
; central pixels, but has slight errors on pixels within 2 rows of the edge.
; This is presumably due to how the cubic interpolater deals with edge
; effects. For now I am going to IGNORE this.
;TODO how well will this deal with IRCAL's undersampled data? Is there something
; we can do which will improve this?
; ANS: it definitely rings off of bad pixels. But it does OK on other stuff,
; as far as I can tell. So do the best job of removing hot/cold/bad pixels
; first and this should be OK.
;im_out = interpolate(image,ix,iy,cubic=-0.5,/grid)
im_out = interpolate(image,ix,iy,cubic=cubic,/grid)
if ~(keyword_set(nofluxscale)) then begin
; Rescale the image to preserve the total flux
;
; Pixel area conversion factor, to preserve flux
pixelarearatio = pixelscale^2 / (scalex*scaley)
im_out *= pixelarearatio
endif
; TODO James points out that for just rescaling the axes, with no spatially
; varying distortions, then one can use the Jacobian to compensate for the
; change of variables between the different coordinate systems.
; See eqn 6 on http://mathworld.wolfram.com/Jacobian.html
; But in practice here - we distort all pixels equally, so it's just some
; scalar coefficient, and then we measure our photometry on standards
; reduced using the code, so the fixed scalar coefficients should cancel out
; in the end. - MDP 2006-06-13
if nanct gt 0 then begin
newwnan = where(newnanmask eq 0)
im_out[newwnan] = !values.f_nan
endif
if keyword_set(byte) then im_out = byte(round(im_out))
return,im_out
end