; Copyright (c) 1995, ALS SpectroMicroscopy Facility, jdd 
; 
;FILE:  structure.pro 
;+ 
;2D Sructure procedures for Microscopy images 
;Current Structure format: 
;   s = {D:data, X:xarr, Y:yarr, XL:xlabel, YL:ylabel, DL:title} 
; 
;  pro SETORIGIN, s [, x0, y0, /CENTER, /CORNER, /PLOT] 
;       -- reset origin of x & y axes; 
;          prompts for cursor input if location not specified 
;  pro SETSCALE, s [, opt,  /CENTER, /PLOT] 
;      -- reset x & y scale to image (assumes dx=dy) according to: 
;         WIDTH=width;   INCR=increment;  or 
;         LINE=length (place arbitrary line with cursors) 
;  fct SPROFILE, s		-- extract line profile from image structure 
;  fct SPROFILES, s [,opt]	-- interactive XYPROFILES 
;  fct SCROP, s [, xrange, yrange] 
;      -- crop image using box_cursor unless 
;         xrange=[x0,x1] and yrange=[y0,y1] are specified 
;  fct STAG, s, tag(s)	-- checks existence of structure field(s) 
; -- Misc -- 
;  fct ARR2STRUC, arr		-- convert 2D array to structure 
;  fct SCHECK, s 
;  fct SINTERP, s, xx, yy, MISSING=ms 
;     -- interpolate 2D structure to new X and Y arrays 
;  fct SCUTOUT, s, [x0,x1], [y0,y1], INC=inc, , _extra=e 
;     -- interpolate 2D structure to new X and Y ranges [min,max] 
;      default: uses original increment unless new INC specified 
;  pro CVTUNIT, s, TO='newunit' 
 
;  pro SUBST, s, new, /X, /Y, /UNIT 	-- substitute new item into structure 
;  fct SMOD, s, new, /X, /Y, /UNIT 
;  pro SSAVE, file, s		-- save 1D structure as columns in text file 
;
; (11-nov-00 cgz)
;	NOTE : Nothing other than STAG is used in AXIS.  
;		Created STAG.PRO.  STRUCTURE.PRO is obsolete and can be eliminated
;
;- 
 
 
function STAG, s, tags 
;====================== 
; check existence of structure field (or tag); 
; returns T/F (0/1) byte or byte array 
on_error,2 
ss=size(s) 
if ss(0)*ss(2) ne 8 then message, 'input not a structure' 
 
fields=strlowcase( tag_names(s) ) 
fstr='' 
for i=0, n_tags(s)-1 do fstr=fstr+fields(i)+';' 
;print, fstr 
 
tags=strlowcase(tags) 
nt=n_elements(tags) 
exist=bytarr(nt) 
for i=0, nt-1 do exist(i)=strpos( fstr, tags(i) ) gt -1 
;print, exist 
;exist= fix( (exist+.1)/abs(exist+.1) )>0 
if nt eq 1 then return, exist(0) $ 
else return, exist 
end 
 
pro subst, s, new, X=ix, Y=iy 
;******************* 
; structure substitute item PROCEDURE 
on_error,2 
sbackup=s 
mode='D' 
if keyword_set(ix) eq 1 then mode='X' 
if keyword_set(iy) eq 1 then mode='Y' 
case mode of 
'D':  s.d=new 
'X':  s.x=new 
'Y':  s.y=new 
endcase 
 
if scheck(s) gt 0 then s=sbackup	; no change if error 
end 
 
function smod, s, new, X=ix, Y=iy, DL=dl, XL=xl, YL=yl 
;==================================== 
; structure substitute item FUNCTION 
on_error,2 
snew=s 
mode='D' 
if keyword_set(ix) eq 1 then mode='X' 
if keyword_set(iy) eq 1 then mode='Y' 
if keyword_set(xl) eq 1 then mode='XL' 
if keyword_set(yl) eq 1 then mode='YL' 
if keyword_set(dl) eq 1 then mode='DL' 
case mode of 
'D':  snew.d=new 
'X':  snew.x=new 
'Y':  snew.y=new 
'XL':  snew.xl=new 
'YL':  snew.yl=new 
'DL':  snew.dl=new 
endcase 
 
if scheck(snew) gt 0 then snew=s 
return, snew 
end 
 
function scheck, s 
;************** 
;consistency check of dimensions ;and valid units 
on_error,2 
ierr=0 
siz=size(s.d)	;get num dimensions 
if siz(0) ne 2 then begin 
   print,'STRUCTURE ERROR: Data not 2-dimensional!' 
   ierr=1 
endif 
if ncols(s.x) ne ncols(s.d) then begin 
   print,'STRUCTURE ERROR: Incompatible X dimensions!' 
   ierr=2 
endif 
if ncols(s.y) ne nrows(s.d) then begin 
   print,'STRUCTURE ERROR: Incompatible Y dimensions!' 
   ierr=3 
endif 
;if strpos('pt;in;cm;mm;um;nm', s.u) eq -1 then begin 
;   print,'STRUCTURE ERROR: Invalid unit!' 
;   ierr=4 
;endif 
return, ierr 
end 
 
pro SETSCALE, s, WIDTH=width, INCR=incr, UNIT=unit, LINE=csrlen, $ 
      CENTER=ctr, PLOT=pl 
;---------------------------------------------- 
;assumes unit aspect ratio of pixel increment 
on_error,2 
sbackup=s 
nx=ncols(s.x) 
ny=ncols(s.y) 
 
mode=0    ; no value change 
if keyword_set(csrlen) eq 1 then begin 
; interactively obtain length in pixels from two cursors 
; assumes image is displayed with SPLOT2D 
  ; mag=floor(400/nx) 
  ; nwin,20 
  ; tvscln, s.d, mag, /WIN 
  ; pixlen=CLENGTH()/mag 
  ; increment=float(csrlen)/float(pixlen-1.) 
  ; wdelete,20 
   oldwidth=abs(max(s.x)-min(s.x)) 
   increment=float(oldwidth)/float(nx-1.) 
;   scale=float(csrlen)/CLENGTH() 
   increment=increment*scale 
  ; print, oldwidth, increment, scale 
   mode=1 
endif 
if keyword_set(width) eq 1 then begin 
   increment=float(width)/float(nx-1.) 
   mode=1 
endif 
if keyword_set(incr) eq 1 then begin 
   increment=float(incr) 
  mode=1 
endif 
if mode eq 1 then begin 
   ;print, increment 
   s.x=findgen(nx)*increment 
   s.y=findgen(ny)*increment 
endif 
 
;if keyword_set(unit) eq 1 then begin   ; rewrite unit 
;   s.u=STRTRIM(unit) 
;endif 
 
if scheck(s) gt 0 then s=sbackup     ; no change if error 
 
if keyword_set(ctr) eq 1 then  setorigin, s, /CENTER 
if keyword_set(pl) then splot2d, s 
end 
 
pro SETORIGIN, s, x0, y0, CENTER=ctr, CORNER=cnr,  PLOT=pl 
;--------------------------------------------- 
;assumes unit aspect ratio of pixel increment 
;uses cursor input if no location specified 
on_error,2 
sbackup=s 
 
ierr=0 
if keyword_set(ctr) eq 1 then  begin 
   x0=(max(s.x)+min(s.x))/2. 
   y0=(max(s.y)+min(s.y))/2. 
endif 
 
if keyword_set(cnr) eq 1 then  begin 
   x0=s.x(0) 
   y0=s.y(0) 
endif 
 
if n_elements(x0) eq 0 then  begin 
;assumes image is displayed via splotx2d; use data coordinates 
   print,'Select new origin within image using cursor; assumes use of SPLOT2D' 
   tvcrs, !d.x_size/2, !d.y_size/2 
   cursor, x0, y0, /data 
endif 
 
if ierr eq 0 then begin 
   s.x=s.x-x0 
   s.y=s.y-y0 
endif 
if scheck(s) gt 0 then s=sbackup     ; no change if error 
if keyword_set(pl) then splot2d, s 
end 
 
function SINTERP, s, xx, yy, MISSING=ms, _extra=e 
;------------------------------------- 
; bilinear interpolate 2D structure to new X and Y arrays 
; outputs new structure with same units 
; /GRID returns 2D structure 
on_error,2 
nx=ncols(xx) 
ny=ncols(yy) 
new=fltarr(nx, ny) 
; map new xy values to fractional array indices 
ix=INTERPOL( findgen(ncols(s.x)), s.x, xx) 
iy=INTERPOL( findgen(ncols(s.y)), s.y, yy) 
;print, ix,iy 
if keyword_set(ms) eq 1 then $ 
   data=INTERPOLATE( s.d, ix, iy, _extra=e, MISSING=ms) $ 
else $ 
   data=INTERPOLATE( s.d, ix, iy, _extra=e) 
 
snew={d:data, x:xx, y:yy, xl:'', yl:'', dl:''} 
return, snew 
end 
 
function SCUTOUT, s, xr, yr, INC=inc, _extra=e 
;============================================= 
;     ** interpolate 2D structure to new X and Y ranges [min,max] 
;      default: uses original increment unless new INC specified 
;      converts XY ranges to arrays and calls SINTERP 
if keyword_set(inc) eq 1 then begin		;INCREMENT SPECIFIED 
   inc=float(inc) 
; accept both single value and vector increment inputs 
   ni=ncols(inc) 
   if ni eq 2 then incr=inc else incr=[inc, inc] 
endif else begin 
; USE INCREMENT FROM ORIGINAL (assume linear) 
   incr=[ abs(s.x(1)-s.x(0)), abs(s.y(1)-s.y(0)) ] 
endelse 
nx=floor( abs(xr(1)-xr(0))/incr(0) ) 
ny=floor( abs(yr(1)-yr(0))/incr(1) ) 
xx=findgen(nx)*incr(0) + xr(0) 
yy=findgen(ny)*incr(1) + yr(0) 
;print,xx,yy 
 
snew=SINTERP( s, xx, yy, /grid, _extra=e) 
return, snew 
end 
 
 
function SCROP, s, xr, yr, MAG=mag, DEV=dev 
;================================== 
;  cutout rectangular section of structure 
; use box_cursor in DATA coordinates unless input 
; ranges [x0,x1], [y0,y1] are specified 
;  if using TVSCL and device coordinates, MAG=magnification is important 
on_error,2 
if n_elements(yr) eq 0 then begin  ; no input ranges specified 
if keyword_set(dev) then begin 
   print,'image must be displayed with TVSCL or TVSCLN 
   box_cursor, px0, py0, pnx, pny, /MESSAGE 
   if keyword_set(mag) eq 0 then mag=1   ; default unit magnification 
   px=floor( [px0, px0+pnx]/mag ) 
   py=floor( [py0, py0+pny]/mag ) 
endif else begin 
   box_cursor2, x0, y0, nx, ny, /MESSAGE   ;data coordinates 
   ;print, x0, y0, nx, ny 
   ; convert to pixels 
   sz=size(s.d)		; assume s.x, s.y agree in dimension 
   ;print, sz 
   px=round( interpol(indgen(sz(1)), s.x, [x0, x0+nx]) ) 
   py=round( interpol(indgen(sz(2)), s.y, [y0, y0+ny]) ) 
endelse 
endif else begin	;input range given 
   if n_elements(xr) ne 2 then message, 'input ranges as 2-vector' 
   if n_elements(xr) ne 2 then message, 'input ranges as 2-vector' 
if keyword_set(dev) then begin  ; DEVICE (pixel) range given 
   px=xr 
   py=yr 
endif else begin	; DATA coordinate range given 
   sz=size(s.d)		; assume s.x, s.y agree in dimension 
   px=round( interpol(indgen(sz(1)), s.x, xr) ) 
   py=round( interpol(indgen(sz(2)), s.y, yr) ) 
endelse 
endelse 
;print, px, py 
data=s.d( px(0):px(1), py(0):py(1) ) 
xaxis=s.x( px(0):px(1) ) 
yaxis=s.y( py(0):py(1) ) 
snew={d:data, x:xaxis, y:yaxis, xl:s.xl, yl:s.yl, dl:s.dl} 
print,'New image: ',ncols(snew.d),' x ',nrows(snew.d) 
return, snew 
end 
 
pro cvtunit, s, TO=nunit 
;------------------------------ 
; rescale X and Y values to new unit 
on_error,2 
 
if n_elements(s.u) eq 0 then message, 'No unit field (s.u) in structure' 
if keyword_set(nunit) eq 0 then begin 
   print,'Input new unit string (pt,in,cm,mm,um,nm):' 
   read, nunit 
endif 
 
if strpos('pt;in;cm;mm;um;nm', nunit) eq -1 then begin 
   print,'Invalid unit!  pt,in,cm,mm,um,nm supported' 
endif else begin 
 
case nunit of 
'pt':  nscf=1. 
'in':  nscf=1./25.4 
'cm':  nscf=0.1 
'mm':  nscf=1.	;reference unit 
'um':  nscf=1000. 
'nm':  nscf=1.E6 
endcase 
 
ounit=s.u 
case ounit of		;old unit 
'pt':  oscf=1. 
'in':  oscf=1./25.4 
'cm':  oscf=0.1 
'mm':  oscf=1. 
'um':  oscf=1000. 
'nm':  oscf=1.E6 
endcase 
 
scale=nscf/oscf 
s.x=s.x*scale 
s.y=s.y*scale 
s.u=nunit 
print, 'converted ',ounit,' to ', nunit,', scale factor=',scale 
 
endelse 
end 
 
function arr2struc, arr, XL=xl, YL=yl, DL=dl 
;=========================================== 
; convert 2D array to structure with default point scaling 
xarr=findgen(ncols(arr)) 
yarr=findgen(nrows(arr)) 
if not keyword_set(xl) then xl='' 
if not keyword_set(xl) then yl='' 
if not keyword_set(xl) then dl='' 
s={d:arr, x:xarr, y:yarr, xl:xl, yl:yl, dl:dl, u:'pt'} 
return, s 
end 
 
function SPROFILE, s 
;=========================== 
; extract line profile from structure 
; put into a 1-dimensional structure 
; 
on_error,2 
;print,'image must be displayed as TVSCL' 
;if keyword_set(mag) eq 0 then mag=1   ; default unit magnification 
;lin=profile(zm(s.d,mag), xx, yy) 
 
cline, r1, r2, len, ang, /SHOW, /data 
print, 'length', len 
 
p1=convert_coord( r1(0), r1(1), /data, /to_device ) 
p2=convert_coord( r2(0), r2(1), /data, /to_device ) 
np=round( sqrt( (p2(0)-p1(0))^2 + (p2(1)-p1(1))^2 ) ) 
 
xarr=rangen( np, r1(0), r2(0) ) 
yarr=rangen( np, r1(1), r2(1) ) 
 
sslin=sinterp( s, xarr, yarr ) 
lin=sslin.d 
 
;parr=convert_coord( xarr, yarr, /data, /to_device ) 
;print, parr 
;lin = interpolate( s.d, parr(0,*), parr(1,*) ) 
 
incr=len/float(np-1) 
;rarr=findgen(np)*incr 
rarr=rangen( np, 0, len ) 
 
;nx=ncols(xx) 
;xa=rebin(s.x, ncols(s.x)*mag) 
;ya=rebin(s.y, ncols(s.y)*mag) 
;length=sqrt( (xa(xx(0))-xa(xx(nx-1)))^2 + (ya(yy(0))-ya(yy(nx-1)))^2) 
 
; xaxis=findgen(plen)*incr 
; lin = interpol( s.d, ) 
;unit=s.u 
slin={d:lin, x:rarr, xl:'', dl:''} 
return, slin 
end 
 
pro SPROFILES, s, _extra=e 
;---------------------------- 
; pass along X and Y arrays to XYPROFILES 
; options: /PLOT, MAG=mag 
on_error,2 
XYPROFILES2, s.d, s.x, s.y, _extra=e 
end 
 
pro ssave, file, s, error=err 
;********************* 
;save 1-D structure {d:, x:, dl:, xl:} in 2 column text file 
on_error,2 
np=nrows(s.d) 
openw, iunit, file, /get_lun, error=err 
printf, iunit, 'd: '+s.dl 
printf, iunit, 'x: '+s.xl 
for i=0,np-1 do printf, iunit, s.x(i), s.d(i) 
close, iunit 
 
end 
 
 
 
 
 
 
