; Date: April 27, 1993 
; 
PRO SINUS,ANGLE0,IMAGE,PHAS,AMP,OFFSET,SIGMA,img_avoid=img_avoid
;PHAS is phase in degrees of MINIMUM of sinusoid, not maximum!  
;	Feb-94 (MDM) - Various changes
;	 4-Mar-94 (TDT) - added SIGMA = rms deviation of fit from data (optional)
;			- deleted some unused & inefficient stuff in inner loop
;	 7-Mar-94 (MDM) - Corrected SIGMA definition
;			- Added IMG_AVOID option
;
dr=!pi/180 
dim=size(image) 
n_x = dim(1) 
n_y = dim(2) 
if (keyword_set(img_avoid)) then begin		;MDM 7-Mar-94
    barr = bytarr(n_elements(angle0)) + 1
    barr(img_avoid-1) = 0
    ss = where(barr)
end else begin
    ss = indgen(n_elements(angle0))
end
angle = angle0(ss)	;subset posibily
;
image=float(image)      ;convert to floating point array 
phas=fltarr(n_x,n_y) 
amp=fltarr(n_x,n_y) 
offset=fltarr(n_x,n_y) 
sig2=fltarr(n_x,n_y)
thr=dr*transpose(angle) ;convert angles to radians 
x=dr*angle
c1=cos(thr) 
c2=c1*c1 
s1=sin(thr) 
s2=s1*s1 
mat=dblarr(3,3) 
mat(0,0)=total(c2) 
mat(0,1)=total(s1*c1) 
mat(1,0)=mat(0,1) 
mat(0,2)=total(c1) 
mat(1,1)=total(s2) 
mat(2,0)=mat(0,2) 
mat(1,2)=total(s1) 
mat(2,1)=mat(1,2) 
mat(2,2)=n_elements(angle) 
;mat=double(mat)        ;if double precision needed 
mati=invert(mat) 
for i=0,n_x-1 do begin 
        for j=0,n_y-1 do begin 
                value=image(i,j,ss) 
                rhs=dblarr(1,3) 
                rhs(0,0)=total(value*c1) 
                rhs(0,1)=total(value*s1) 
                rhs(0,2)=total(value) 
                res=rhs#mati 
                alpha=res(0,0) 
                beta=res(0,1) 
                gamma=res(0,2) 
                A=sqrt(alpha^2+beta^2) > 0.01	;MDM added > 0.01 on 8-Feb-94
                B=gamma-A 
                ;phir=atan(-beta/A,-alpha/A) 
                if (abs(alpha)+abs(beta) eq 0) then phir = 0 else phir=atan(-beta/A,-alpha/A) 	;MDM added check 8-Feb-94
                ;;phi=phir/dr 
		;;if (a lt 1) then print,i,j,a,b,alpha,beta,phi
                 y=A*(1.-cos(x-phir))+B 		;MDM redefined 7-Mar-94
                ;;y=A*(1.-cos(2.*(x-phir)))+B
                 amp(i,j)=A 
                offset(i,j)=gamma 
                phas(i,j)=phir 
                sig2(i,j)=total((value-y)^2)
;;		if (i eq 32) and (j eq 32) then stop
        endfor 
endfor 
phas=phas/dr
if (n_params(0) gt 5) then sigma=sqrt(sig2/n_elements(angle))
return 
end 
 
PRO CIRCLEFIT,IMAGE,RAD,RTHRESH 
; Finds center and radius of circular footprint on a set of images 
; Looks at diagonals, using threshold of 5% of central intensity 
; Assumes corners are dark, footprint is bright, x_scale=y_scale 
; Finds median center and crops all images to RTHRESH % radius (def = .95) 
;HISTORY:
;	Written by Ted Tarbell
;	 8-Mar-94 (MDM) - Modified to take the average of all images and
;			  find the edges from that image, rather than 
;			  finding the edges of all images and then
;			  returning the mean of the edges.
THRESH=0.05 
IF (N_PARAMS(0) LT 3) THEN RTHRESH=0.95 
Dim=size(image) 
nx = dim(1) 
ny = dim(2) 
n  = dim(3) 
ND=MIN([NX,NY]) 
NDL=FIX(.25*ND) & NDU=FIX(.75*ND) & NDS=FLOAT(NDU-NDL+1) 
T=FLTARR(ND) & TT=T 

;;EDGE=FLTARR(4,N) 
;;FOR I=0,N-1 DO BEGIN 
;;  FOR J=0,ND-1 DO BEGIN T(J)=IMAGE(J,J,I) & TT(J)=IMAGE(NX-1-J,J,I) & ENDFOR 
;;  TH=FIX(THRESH*TOTAL(T(NDL:NDU))/NDS) 
;;  W=WHERE(T GE TH) 
;;  EDGE(0,I)=MIN(W)  &  EDGE(1,I)=MAX(W)  
;;  TH=FIX(THRESH*TOTAL(TT(NDL:NDU))/NDS) 
;;  W=WHERE(TT GE TH) 
;;  EDGE(2,I)=MIN(W)  &  EDGE(3,I)=MAX(W)  
;;ENDFOR 
;;; for i=0,n-1 do print,edge(*,i) 
;;E0=MEDIAN(EDGE(0,*)) & E1=MEDIAN(EDGE(1,*))  
;;E2=MEDIAN(EDGE(2,*)) & E3=MEDIAN(EDGE(3,*)) 

temp = total(image,3)/n
FOR J=0,ND-1 DO BEGIN T(J)=temp(J,J) & TT(J)=temp(NX-1-J,J) & ENDFOR
bkg = total(t(0:2) + t(nd-4:nd-1) + tt(0:2) + tt(nd-4:nd-1))/4/3
t = t - bkg
tt = tt - bkg
TH=FIX(THRESH*TOTAL(T(NDL:NDU))/NDS) 
W=WHERE(T GE TH) 
e0 = MIN(W)  &  e1=MAX(W)  
TH=FIX(THRESH*TOTAL(TT(NDL:NDU))/NDS) 
W=WHERE(TT GE TH) 
e2 = MIN(W)  &  e3=MAX(W)  

UC=(E0+E1)/2.  &  VC=(E2+E3)/2. 
RAD=FLTARR(3) 
RAD(0)=1.414*(E1-E0+E3-E2)/4. 
RAD(1)=UC+NX/2.-VC  &  RAD(2)=UC-NX/2.+VC 
; print,'E0 - E3, RAD',e0,e1,e2,e3,rad 
TH=(RTHRESH*RAD(0))^2 
X2=(FINDGEN(NX)-RAD(1))^2  &  Y2=(FINDGEN(NY)-RAD(2))^2   
MASK=FLTARR(NX,NY) 
FOR I=0,NY-1 DO MASK(*,I)=Y2(I)+X2 
TH=(RTHRESH*RAD(0))^2 
W=WHERE(MASK GT TH) &  WW=WHERE(MASK LT TH) 
NW=FLOAT(N_ELEMENTS(WW)) 
if (ww(0) ne -1) then FOR I=0,N-1 DO BEGIN  
  MASK=IMAGE(*,*,I) 
  AVERAGE=TOTAL(MASK(WW))/NW 
  MASK(W)=AVERAGE 
  IMAGE(*,*,I)=MASK 
  IMAGE(0,*,I)=AVERAGE 
  IMAGE(NX-1,*,I)=AVERAGE 
  IMAGE(*,0,I)=AVERAGE 
  IMAGE(*,NY-1,I)=AVERAGE 
ENDFOR 
RETURN 
END 
  
PRO DEFRINGE,IMAGE,FRINGE,FWHM=FWHM0,HIFWHM=HIFWHM0
;    Defringes each image at spatial frequency index 
;    Fx = FRINGE(0), Fy = FRINGE(1) (cycles per side of image) 
;    Uses 1 - Gaussian filter of width FWHM (default = 12)  
;    Also does low pass filter, using Gaussian
;     of width HIFWHM (default = NX/2)
;  STILL HAS TV DISPLAYS IN IT!
;HISTORY:
;	Written by Ted Tarbell
;	 7-Mar-94 (MDM) - Turned FWHM and HIFWHM into keywords
;

if (n_elements(fwhm0) ne 0) then fwhm=fwhm0 else FWHM=12.
Dim=size(image) 
nx = dim(1) 
ny = dim(2) 
n  = dim(3) 
NP=FLOAT(NX)*NY 
NX2 = NX/2  &  NY2=NY/2
if (n_elements(hifwhm0) ne 0) then hifwhm=hifwhm0 else HIFWHM=NX

; Defringing flter
WID=4.*ALOG(2.)/FWHM^2 
X2=FINDGEN(NX)-FRINGE(0)  &  X2=WID*X2^2 
Y2=FINDGEN(NY)-FRINGE(1)  &  Y2=WID*Y2^2 
W=0.*FLTARR(NX,NY) 
FOR I=0,NY2 DO W(*,I)=1.-EXP(-(X2+Y2(I))) 

;  Low pass filter
WID=4.*ALOG(2.)/HIFWHM^2 
X2=FINDGEN(NX2+2)  &  X2=WID*X2^2 & XG = EXP(-X2)
Y2=FINDGEN(NY2+2)  &  Y2=WID*Y2^2 & YG = EXP(-Y2) 
REVNX2 = NX-NX2-1-INDGEN(NX-NX2-1)
FOR I=0,NY2 DO BEGIN
      W(0:NX2,I) = W(0:NX2,I)*XG(0:NX2)*YG(I)
      W(NX2+1:*,I) = W(NX2+1:*,I)*XG(REVNX2)*YG(I)
      ENDFOR
W(*,NY2+1:NY-1)=ROTATE(W(*,1:NY-NY2-1),2) 
if (getenv('ys_batch') eq "") then begin
    plot_io,w(*,0) & oplot,w(0,*)
    tvscl,w,2 
end

FOR I=0,N-1 DO BEGIN 
  T=IMAGE(*,*,I) 
  if (getenv('ys_batch') eq "") then tvscl,image(*,*,i),0 
  AVER=TOTAL(T)/NP 
  TT=FFT(T-AVER,1) 
  T=FFT(TT*W,-1) 
  IMAGE(*,*,I)=AVER+FLOAT(T) 
  if (getenv('ys_batch') eq "") then tvscl,image(*,*,i),1 
ENDFOR 
RETURN 
END 
 
PRO LOWPASS,IMAGE,FCUTOFF 
; Fourier lowpass filter deleting power above FCUTOFF cycles 
;   across the images 
Dim=size(image) 
nx = dim(1) 
ny = dim(2) 
n  = dim(3) 
NP=FLOAT(NX)*NY 
X2=FINDGEN(NX) &  X2(NX/2:NX-1)=NX-X2(NX/2:NX-1)  &  X2=X2^2 
Y2=FINDGEN(NY) &  Y2(NY/2:NY-1)=NY-Y2(NY/2:NY-1)  &  Y2=Y2^2 
F2=FLOAT(FCUTOFF)^2 
W=FLTARR(NX,NY) 
FOR I=0,NY-1 DO W(*,I)=X2+Y2(I) 
WW=WHERE(W GT F2) 
FOR I=0,N-1 DO BEGIN 
  T=IMAGE(*,*,I) 
  AVER=TOTAL(T)/NP 
  TT=FFT(T-AVER,1) 
  TT(WW)=COMPLEX(0,0) 
  T=FFT(TT,-1) 
  IMAGE(*,*,I)=AVER+FLOAT(T) 
ENDFOR 
RETURN 
END 
 
PRO SINFIT,ANGLE,IMAGE,PHAS,AMP,OFFSET,SIGMA,RAD,FRINGE, $
		FWHM=FWHM,HIFWHM=HIFWHM,img_avoid=img_avoid
;
;NAME:
;	sinfit
;PURPOSE:
;	Fits sinusoid to each pixel of image(x,y,*) 
;INPUT:
;	angle	- The angle in degrees
;	image	- 3-D array
;OUTPUT:
;	PHAS 	- is phase in degrees of MINIMUM of sinusoid, not maximum! 
;		  (It is a 2-D array)
;	AMP	- amplitude of sine fit (It is a 2-D array)
;	OFFSET	- offset of sine fit (It is a 2-D array)
;	SIGMA	- sigma (fit quality) of sine fit (It is a 2-D array)
;	RAD	- RAD(0) = radius of circular image, in pixels 
;		  RAD(1) = x coordinate of center
;		  RAD(2) = y coord. 
;OPTIONAL KEYWORD INPUT:
;	FRINGE	- ??
;	FWHM	- ??
;	HIFWHM	- ??
;	img_avoid - An optional list of image numbers to avoid during the
;		    fit.  If not passed, use all images passed in.
;METHOD:
;	Uses quick & dirty centering routine 
;	If FRINGE is given, defringes each image  
;HISTORY:
;	Written by Ted Tarbell
;	 4-Mar-94 (TDT) - Added SIGMA output
;	 7-Mar-94 (MDM) - Converted to have keywords for
;			  FRINGE, FWHM and HIFWHM
;

CIRCLEFIT,IMAGE,RAD 
;;defringe, image, fringe, fwhm=fwhm, hifwhm=hifwhm
SINUS,ANGLE,IMAGE,PHAS,AMP,OFFSET,SIGMA,img_avoid=img_avoid
return 
END