Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.stsci.edu/~mperrin/software/sources/cube_collapse.pro Дата изменения: Tue Jun 22 05:24:02 2010 Дата индексирования: Sat Mar 1 16:57:13 2014 Кодировка:

;+
; NAME: cube_collapse
; PURPOSE:
; Collapse a cube along the wavelength axis.
;
; The default is to average the cube, but there are options to median or
; total combine instead.
;
; There is some support for using quality flags to indicate bad pixels to skip
; when collapsing.
;
; INPUTS:
; cubestruct a cubestruct
; KEYWORDS:
; zrange= which Z axis range to include? Default is all.
; cropmult= crop the collapsed image to multiples of this many pixels
; (for consistency with what cube_rebin does for odd axis
; sizes)
; /median median collapse it
; /total total it, instead of averaging or medianing.
; /silent suppress printed output
; OUTPUTS:
; a modified cubestruct (really an image structure) which contains the
; collpased image, header, and astrometry
;
; HISTORY:
; Began 2007-07-16 17:40:19 by Marshall Perrin
; 2008-02-06 Some support for zrange and quality flags added
;-


FUNCTION cube_collapse, cubestruct, zrange=zrange, silent=silent, cropmult=cropmult, $
median=median,total=total, skipastr=skipastr

waveaxis = cube_getwaveaxis(cubestruct)

if ~(keyword_set(silent)) then print, "Collapsing along wavelength axis "+strc(waveaxis)
if keyword_set(median) then print, "Collapsing using MEDIAN"
if ~tag_exist(cubestruct,'ASTR') then skipastr=1
sz = size(cubestruct.cube)

; TODO add in check for, and use of, optional quality flag element

if keyword_set(total) then method="total"
if keyword_set(median) then method="median"
if ~(keyword_set(method)) then method="mean"


if tag_exist(cubestruct, "flags") then begin
; Use quality flags to ignore bad pixels
if keyword_set(zrange) then begin
;message, "Can't do zrange with quality flags yet"
;if waveaxis ne 2 then message, "Can't do zrange with quality flags, unless wave axis is 3"

; This is semi-convoluted:
; - Transpose the data to have the wavelength axis last.
; - Apply the Z range
; - then transform it back to the original orientation, and
; proceed.

nax = indgen(3)
trans_ax = [nax[where(nax ne waveaxis)], waveaxis]
tmpcube = transpose(cubestruct.cube, trans_ax)
tmpmask = transpose(cubestruct.flags, trans_ax)

tmp = tmpcube[*,*,zrange[0]:zrange[1]]
mask = tmpmask[*,*,zrange[0]:zrange[1]]
back_trans_ax= [where(trans_ax eq 0), where(trans_ax eq 1), where(trans_ax eq 2)]
tmp = transpose(tmp, back_trans_ax)
mask= transpose(mask, back_trans_ax)

endif else begin
mask = (cubestruct.flags ne 0)
tmp = cubestruct.cube
endelse
case method of
"median": begin
tmp[where(mask eq 0)] = !values.f_nan
collapsed = median(tmp, dimen=waveaxis+1)
end
"mean": collapsed = total(tmp*mask, waveaxis+1,/nan) / total(mask, waveaxis+1)
"total": collapsed = total(tmp*mask, waveaxis+1,/nan)
else:
endcase
valid = total (mask, waveaxis+1) ne 0


endif else begin
; Without quality flags, assume there are no bad pixels to ignore.
; Hence include all pixels within zrange.

if keyword_set(zrange) then begin
if ~(keyword_set(silent)) then message,/info, "Using Z range "+printcoo(zrange) +" and method= "+method
subscripts = fltarr(3,2)
subscripts[*,1] = sz[1:3]-1
subscripts[waveaxis,*] = zrange
to_collapse = cubestruct.cube[subscripts[0,0]:subscripts[0,1],subscripts[1,0]:subscripts[1,1],$
subscripts[2,0]:subscripts[2,1]]
case method of
"median": collapsed = median(to_collapse, dimen=waveaxis+1)
"mean": collapsed = total(to_collapse, waveaxis+1) / (zrange[1]-zrange[0]+1)
"total": collapsed = total(to_collapse, waveaxis+1)
else:
endcase
endif else begin
case method of
"median": collapsed = median(cubestruct.cube, dimen=waveaxis+1)
"mean": collapsed = total(cubestruct.cube, waveaxis+1) / sz[waveaxis+1]
"total": collapsed = total(cubestruct.cube, waveaxis+1)
else:
endcase
endelse
endelse

if keyword_set(cropmult) then begin
if ~(keyword_set(silent)) then message,/info, "cropping image size to multiples of "+strc(cropmult)
multsz = mround(sz,cropmult)
wa = where(indgen(n_elements(cubestruct.astr.naxis)) ne waveaxis)
;multsz[waveaxis+1] = sz[waveaxis+1]
collapsed = collapsed[0:multsz[wa[0]+1]-1, 0:multsz[wa[1]+1]-1]

endif

if ~(keyword_set(skipastr)) then begin
; update astrometry info!
newhead = cubestruct.header
astr = cubestruct.astr
wa = where(indgen(n_elements(astr.naxis)) ne waveaxis)
newnaxes = n_elements(astr.naxis)-1

; crop out the right part of the CD matrix. I can't figure out any
; clever way to do this sans for loop. ??
newcd = dblarr(newnaxes, newnaxes)
for ia=0L,newnaxes-1 do begin
newcd[ia,*] = astr.cd[wa[ia], wa]
endfor

NEWASTR = {NAXIS:astr.naxis[wa], CD: astr.cd, CDELT: astr.cdelt[wa], $
CRPIX: astr.crpix[wa], CRVAL:astr.crval[wa], $
CUNIT: astr.cunit[wa], $
CTYPE: string(astr.ctype[wa]), LONGPOLE: double( astr.longpole), $
LATPOLE: double(astr.latpole), PV2: astr.pv2 }
newheader = cubestruct.header
putast3,newheader,newastr
endif else begin
newheader = cubestruct.header
if tag_exist(cubestruct, 'ASTR') then newastr = cubestruct.astr
endelse


sxaddpar, newheader, 'NAXIS', 2,/save
sz = size(collapsed)
sxaddpar, newheader, 'NAXIS1', sz[1],/save
sxaddpar, newheader, 'NAXIS2', sz[2],/save
sxdelpar, newheader, 'NAXIS3'


retval = {image: collapsed, header: newheader}
if keyword_set(newastr) then retval = create_struct(retval, "ASTR", newastr)
if keyword_set(valid) then retval = create_struct(retval, "Flags", valid)

return, retval


end