; Date: April 27, 1993 
; 
PRO SINUS,ANGLE,IMAGE,PHAS,AMP,OFFSET 
;PHAS is phase in degrees of MINIMUM of sinusoid, not maximum!  
dr=!pi/180 
dim=size(image) 
n_x = dim(1) 
n_y = dim(2) 
n   = dim(3) 
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) 
thr=dr*transpose(angle) ;convert angles to radians 
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,*) 
                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) 
                B=gamma-A 
                phir=atan(-beta/A,-alpha/A) 
                phi=phir/dr 
                x=angle 
                y=A*(1.-cos(2.*(x*dr-phir)))+B 
                xmod=A/(A+B) 
                amp(i,j)=A 
                offset(i,j)=gamma 
                phas(i,j)=phi 
        endfor 
endfor 
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) 
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,*)) 
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)) 
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,HIFWHM 
;    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!

IF (N_PARAMS(0) LT 3) THEN 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_PARAMS(0) LT 4) THEN 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) 
plot_io,w(*,0) & oplot,w(0,*)
tvscl,w,2 

FOR I=0,N-1 DO BEGIN 
  T=IMAGE(*,*,I) 
  tvscl,image(*,*,i),0 
  AVER=TOTAL(T)/NP 
  TT=FFT(T-AVER,1) 
  T=FFT(TT*W,-1) 
  IMAGE(*,*,I)=AVER+FLOAT(T) 
  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,RAD,FRINGE,FWHM,HIFWHM
; Fits sinusoid to each pixel of image(x,y,*) 
; PHAS is phase in degrees of MINIMUM of sinusoid, not maximum! 
; Returns RAD(0) = radius of circular image, in pixels 
;         RAD(1) = x coordinate of center, RAD(2) = y coord. 
; Uses quick & dirty centering routine 
; If FRINGE is given, defringes each image  
 
CIRCLEFIT,IMAGE,RAD 
IF (N_PARAMS(0) GT 8) THEN DEFRINGE,IMAGE,FRINGE,FWHM,HIFWHM
IF (N_PARAMS(0) GT 7) THEN DEFRINGE,IMAGE,FRINGE,FWHM 
IF (N_PARAMS(0) GT 6) THEN DEFRINGE,IMAGE,FRINGE 
SINUS,ANGLE,IMAGE,PHAS,AMP,OFFSET 
RETURN 
END