Viewing contents of file '../idllib/astron/contrib/varosi/code/allpro/calc_im_stack.pro'
;+
; NAME:
; calc_im_stack
; PURPOSE:
; Compute an algebraic expression (entered by user)
; involving the intersection regions of a stack of mosaic images.
; CALLING:
; image_result = calc_im_stack( image_List, calc, names, imsel )
; INPUT:
; image_List = array of structures containing image location info,
; and either the image data or pointers to image data.
; KEYWORDS:
; DATA = pooled-array containing image data (e.g. image_data).
; OUTPUT:
; Function returns the result of algebraic expression of images.
; OPTIONAL OUTPUTS:
; calc = string, the algebraic expression that was computed.
; names = string array, the names of images used in computation.
; imsel = integer array, indices of images used in computation.
; EXTERNAL CALLS:
; pro overlap_images
; pro write_images
; pro printw
; pro box_erase
; function box_create
; function Trp3D
; function intersection
; function N_struct
; function round_off
; function conv_ascii_time
; function next_word
; function get_image
; function frac_pix_extrac
; PROCEDURE:
; Parse the text typed in by user, substitute arrays for capital letters,
; extract the images, then execute the statement in current environment.
; HISTORY:
; Written: Frank Varosi NASA/GSFC 1990.
; F.V.1991, compute arbitrary algebraic expression, as entered by user.
; F.V.1992, use new function frac_pix_extrac( image ).
; F.V.1993, added MINIMAL, SPECIFIC, and VARIABLE intersection options.
; F.V.1994, modified to reject data less than thresholds specified by >.
; F.V.1994, added keyword options which bypass menus.
;-
function calc_im_stack, image_List, calc, names, imsel, DATA=image_data, $
INTERSECT=intersecop, REGION=sizreg, CALC_OPTION=calc_op
Nmath = check_struct( image_List, Magf, sizims )
if (Nmath LE 1) then begin
print,' need more than one image for overlap'
return,[-1]
endif
sz = size( sizims )
if (sz(0) EQ 1) then sizims = sizims # replicate( 1, Nmath )
if keyword_set( intersecop ) then begin
intersecop = strupcase( intersecop )
endif else begin
menu = [" ", "MINIMAL overlap (intersection of all images)", $
"SPECIFIC intersection (depending on chosen images)", $
"VARIABLE intersection (with user selected box)", $
"SIMPLE intersection (assume images aligned at (0,0))" ]
intersecop = next_word( menu( wmenu( menu,INIT=3,TIT=0 ) > 0 ) )
verbose = 1
set_cursor = 1
if strlen( intersecop ) LE 1 then return,[-1]
endelse
wset, image_List(0).windo
wshow, image_List(0).windo
LF = string( 10b )
BELL = string( 7b )
if (intersecop EQ "VARIABLE") then begin
printw,["select corners of intersection box " + $
"with LEFT -> MIDDLE buttons, ", $
"or RIGHT button to QUIT" ]
if keyword_set( set_cursor ) then tvcrs,0.5,0.5,/NORM
BOXC: CASE box_create( xL,yL, xH,yH ) OF
(-2): return,[-2]
(-4): BEGIN
if NOT keyword_set( set_cursor ) then return,[-4]
sizebox = get_words( get_text_input( "size of box ?" ) )
if max( strlen( sizebox ) ) GT 0 then begin
sizreg = fix( sizebox )
if N_elements( sizreg ) EQ 1 then $
sizreg = [ sizreg, sizreg ]
if min( sizreg ) LE 1 then return,[-1]
endif else return,[-1]
xL = !D.x_size/2
yL = !D.y_size/2
s = sizreg * Magf
box_erase
box_cursor, xL,yL, s(0),s(1), /ADJACENT,/CONTIN,/INIT
box_draw, POS_XY=[xL,yL], SIZE_XY=s
xL = round( xL/float( Magf ) )
yL = round( yL/float( Magf ) )
posxy = [ [xL,xL+sizreg(0)], [yL,yL+sizreg(0)] ]
END
else: BEGIN
posxy = round( [ [xL,xH], [yL,yH] ]/Magf )
sizreg = reform( posxy(1,*)-posxy(0,*) )
END
ENDCASE
if keyword_set( verbose ) then $
print," size of region (x,y) :",sizreg
if keyword_set( set_cursor ) then printw,[" "," "," "],/ERASE
xmin = image_List.Locx
ymin = image_List.Locy
xmax = xmin + sizims(1,*) -1
ymax = ymin + sizims(2,*) -1
posxyall = Trp3D( [ [ [xmin], [xmax] ], $
[ [ymin], [ymax] ] ], [2,3,1] )
intersecf = intarr( Nmath )
sec_to = intarr( 2, 2, Nmath )
sec_from = fltarr( 2, 2, Nmath )
for i=0,Nmath-1 do begin
intersecf(i) = intersection( posxy, posxyall(*,*,i), $
sex_to, sex_from )
sec_to(0,0,i) = round_off( sex_to )
sec_from(0,0,i) = sex_from + ( sec_to(*,*,i) - sex_to )
endfor
if max( intersecf ) NE 1 then begin
print,LF + $
" region does not intersect any images!" + BELL
print," try again..."
goto,BOXC
endif
printw, replicate( " ", 3 )
endif
if max( tag_names( image_List ) EQ "WAVELENGTH" ) then begin
if (intersecop EQ "VARIABLE") then begin
w = where( intersecf, Nmath )
order = w( sort( image_List(w).waveLength ) )
endif else order = sort( image_List.waveLength )
waveLens = image_List(order).waveLength
wavprint = 1
endif else begin
if (intersecop EQ "VARIABLE") then begin
w = where( intersecf, Nmath )
order = w( sort( image_List(w).name ) )
endif else order = sort( image_List.name )
if max( tag_names( image_List ) EQ "TIME" ) then begin
times = image_List(order).time
endif else $
if max( tag_names( image_List ) EQ "INFO" ) then begin
info = image_List(order).info
times = conv_ascii_time( info.time, info.date )
endif
waveLens = intarr( N_elements( order ) )
wavprint = 0
endelse
if keyword_set( calc_op ) then begin ;check for predefined option.
calc = calc_op
goto, CALC
endif
mins = image_List(order).min
maxs = image_List(order).max
print,BELL
LIST: print, LF + " current images in stack:"
names = image_List(order).name
Lens = strlen( names )
Lenmax = max( Lens )
for i = 0, (Nmath-1)<25 do begin
Letter = string( byte( 65+i ) )
name = names(i)
if (Lens(i) LT Lenmax) then name = name + $
string( replicate( 32B, Lenmax-Lens(i) ) )
range = string( [mins(i),maxs(i)],$
FORM="(' min-max:[',2G9.2,'] ')" )
if (wavprint) then $
range = string( waveLens(i),FORM="(' (',F6.2,'um) ')" ) +range
print," " + Letter + " = " + name + range
endfor
print," enter algebraic expression"
print," using upper case A,B,C,... for the images,"
print," lower case for other IDL functions,"
print," use > and parentheses for threshold rejection, e.g. A/(B>1.2)"
print,"(for more info enter: options)
READ: calc = ""
read, calc
if strlen( calc ) LE 0 then return,[-1]
if (strpos( strupcase( calc ), "LIST" ) GE 0) then goto,LIST
if (strpos( strupcase( calc ), "OP" ) GE 0) then begin
print,LF+" stack commands (up/low case) using ALL images:"
print," WRITE = write stack overlap to a binary file"
print," WRITE F77 = write to a Fortran-77 binary file""
print," WRITE FORM = formatted write to a file""
print," WRITE FITS = write stack to a FITS file""
print," SAVE = extract stack overlap and save to IDL/XDR file"
print," STATS = print SNR, sky, FWHM, of stack overlap"
print," LIST = print list of images again"
goto,READ
endif
CALC:
request = strupcase( calc )
if (strpos( request, "WRITE" ) GE 0) OR $
(strpos( request, "SAVE" ) GE 0) OR $
(strpos( request, "RETURN" ) GE 0) OR $
(strpos( request, "STATS" ) GE 0) then begin
if request EQ "RETURN" then begin
order = rotate( order( sort( image_List(order).Level ) ), 2 )
Nmath = Nmath < 6
order = order(0:Nmath-1)
imsel = order
names = image_List(order).name
endif
if (intersecop EQ "VARIABLE") then begin
sx = reform( sec_to(0,0,order) )
sy = reform( sec_to(0,1,order) )
xyb = reform( sec_from(0,*,order) )
xys = reform( sec_from(1,*,order) ) - xyb
images = make_array( sizreg(0), sizreg(1), Nmath, VAL=-999999, $
TYPE = sizims(sizims(0)+1) )
endif else begin
if (intersecop EQ "SIMPLE") then begin
xb = intarr( Nmath )
yb = xb
xt = replicate( min( sizims(1,*) ),Nmath )-1
yt = replicate( min( sizims(2,*) ),Nmath )-1
endif else overlap_images, image_List, order, xb,xt,yb,yt,/BOX
sizex = round_off( xt-xb )
sizey = round_off( yt-yb )
xyb = transpose( [ [xb], [yb] ] )
xys = transpose( [ [sizex<sizex(0)], [sizey<sizey(0)] ] )
sx = intarr( Nmath )
sy = intarr( Nmath )
images = make_array( sizex(0), sizey(0), Nmath, $
TYPE = sizims(sizims(0)+1) )
endelse
if keyword_set( verbose ) OR !DEBUG then begin
print," extracting..."
help,images
endif
for i=0,Nmath-1 do begin
if !DEBUG then print," " + names(i)," (x0,y0):",xyb(*,i)
images(sx(i),sy(i),i) = frac_pix_extrac( BEG=xyb(*,i), $
SIZE=xys(*,i), $
get_image( order(i), image_List, image_data ) )
endfor
CASE next_word( request ) OF
"RETURN": return, images
"WRITE": BEGIN
CASE next_word( request ) OF
"FITS": write_images, images, waveLens, /FITS, NAMES=names
"FORM": write_images, images, waveLens, /FORMATTED
"F77": write_images, images, waveLens, /F77
else: write_images, images, waveLens
ENDCASE
return,[-1]
END
"SAVE": BEGIN
filename = ""
read," filename to save image stack? ",filename
if (filename NE "") then begin
filesave = filename + ".images"
if N_elements( times ) GT 0 then $
save, images, times, names, FILE=filesave,/VERB,/XDR $
else $
save, images, waveLens, names,FILE=filesave,/VERB,/XDR
print," saved stack of images to: ",filesave
endif else print," nothing saved"
return,[-1]
END
"STATS": BEGIN
CASE next_word( request ) OF
"LOR": getfwhm=3
"G": getfwhm=2
"LIN": getfwhm=1
else: getfwhm=0
ENDCASE
stats = im_stats( N_STRUCT=Nmath )
for i=0,Nmath-1 do begin
stats(i) = im_stats( images(*,*,i), $
NAME=names(i), GET_FWHM=getfwhm )
print_stats, stats(i), TITLE=(i EQ 0)
endfor
return, stats
END
else: return,[-1]
ENDCASE
endif
equation = "calc_image = " + calc
inums = replicate( -1, Nmath )
Th_exp = ""
parens = byte( "()" )
Ncalc = 0
for i=0,Nmath-1 do begin
Letter = string( byte( 65+i ) )
pos = strpos( equation, Letter )
if (pos GE 0) then begin
Letsubs = "images(*,*," + strtrim( Ncalc, 2 ) + ")"
inums(Ncalc) = i
Ncalc = Ncalc+1
if strpos( equation, ">", pos ) EQ (pos+1) then begin
p2 = pos+2
thex = "(" + Letsubs + " LE " + $
strmid( equation, p2, $
strpos( equation,")",pos )-p2 ) + ")"
if total( byte( thex ) EQ parens(0) ) GT $
total( byte( thex ) EQ parens(1) ) then $
thex = thex + ")"
if strlen( Th_exp ) GT 0 then $
Th_exp = Th_exp + " OR " + thex $
else Th_exp = thex
endif
while (pos GE 0) do begin
equation = strmid( equation, 0, pos ) + $
Letsubs + $
strmid( equation, pos+1, 999 )
pos = strpos( equation, Letter )
endwhile
endif
endfor
if (Ncalc LE 0) then begin
print,LF + " expression does not refer to any images!" + BELL
print," remember to use upper case for images, try again..."
wait,1
goto,LIST
endif else inums = inums(0:Ncalc-1)
imsel = order(inums)
if (intersecop EQ "VARIABLE") then begin
sx = reform( sec_to(0,0,imsel) )
sy = reform( sec_to(0,1,imsel) )
xyb = reform( sec_from(0,*,imsel) )
xys = reform( sec_from(1,*,imsel) ) - xyb
images = fltarr( sizreg(0), sizreg(1), Ncalc )
endif else begin
if (intersecop EQ "SPECIFIC") then begin
print," marks are NOT valid in SPECIFIC intersection"
overlap_images, image_List, imsel, xb,xt, yb,yt,/BOX
endif else if (intersecop EQ "SIMPLE") then begin
xb = intarr( Ncalc )
yb = xb
xt = replicate( min( sizims(1,inums) ), Ncalc ) -1
yt = replicate( min( sizims(2,inums) ), Ncalc ) -1
endif else begin
overlap_images, image_List, order, xb,xt, yb,yt, /BOX
xb = xb(inums)
xt = xt(inums)
yb = yb(inums)
yt = yt(inums)
endelse
sizex = round_off( xt-xb+1 )
sizey = round_off( yt-yb+1 )
xyb = transpose( [ [xb], [yb] ] )
xys = transpose( [ [sizex<sizex(0)], [sizey<sizey(0)] ] )
images = fltarr( sizex(0), sizey(0), Ncalc )
sx = intarr( Ncalc )
sy = intarr( Ncalc )
endelse
print," extracting..."
names = strarr( Ncalc )
for i=0,Ncalc-1 do begin
k = imsel(i)
names(i) =string( byte( 65+inums(i) ) )+" = "+image_List(k).name
print," " + names(i)
images(sx(i),sy(i),i) = frac_pix_extrac( BEG=xyb(*,i), $
SIZE=xys(*,i), $
get_image( k, image_List, image_data ) )
endfor
if !DEBUG then stop
if (NOT execute( equation )) then begin
message,"error executing: " + equation + BELL,/INFO
goto,LIST
endif
if strlen( Th_exp ) GT 0 then begin
Th_exp = "wth = where( " + Th_exp + ", Nth )"
if (NOT execute( Th_exp )) then begin
message,"error executing: " + Th_exp + BELL,/INFO
goto,LIST
endif
if (Nth GT 0) then begin
maxim = max( calc_image, MIN=minim )
fill = minim - (maxim - minim)/100.
print," data rejection threshold was specified"
print," computation min, max : ", minim, maxim
print," fill value for pixels with no data =", fill
calc_image( wth ) = fill
endif
endif
print," done computing"
print," " + calc
wait,0.5
box_erase
return, calc_image
end