pro pgisol_s,movin,movout,msk

; Procedure to isolate penumbral grains in a 3-D sunspot movie.
; "Everything" except PGs is set to zero.
;
;Version of PGISOL.PRO - processing of P.Suetterlin's DOT data.
;General mask is set individually
;for each frame, trying to outline the actual penumbra boundary.
;Then a mask to remove too dark BFs (i.e. UDs) is applied.
;Moreover, an external mask MSK can be used.
;The threshold for segmentation THR is a linear function of
;the mean MEA of the difference penumbral image.
;
; Method: Differentiation based on unsharp masking,
;         like in SBV, ApJ 415 (1993), 832
; MOVIN - input movie
; MOVOUT- output movie with isolated PGs (544 x 368 pix)
; MSK   - external mask (optional)
; Adapted from UDISOL2 and PGISOL2, 24 March 2000, Michal

on_error,1

;*************** PARAMETERS **********************************
; A,B   - segmentation threshold coefficients (THR=A+B*MEA)
; SMO   - unsharp masking smoothing parameter (in pixels)
; PETHR - intensity below which the bright BFs are removed
; GMSMO - smoothing parameter to define general mask
; GMTHO - threshold to define outer general mask
; GMTHI - threshold to define inner general mask (used only
;	  to calculate MEA)
 a=21
 b=0.45
 smo=7 		; 0.57"
 pethr=600	; minimum allowed intensity (here Iph = 1219)
 gmsmo=31 	; 2.5" to smooth the penumbra
 gmtho=1070	; max. penumbral intensity in smoothed frame
 gmthi=550
;*************************************************************

if n_params() lt 3 then msk=1	; external mask, default

smov=size(movin)
if smov(0) ne 3 then message,'Input must be a 3-D array!'
movout=intarr(544,368,smov(3))

for i=0,smov(3)-1 do begin

   oima=movin(*,*,i)		; original image
   ima=SMOOTH(float(oima),smo)	; smoothed image (unsharp masking)
   ima=oima-ima			; difference image (unsharp masking)

   ois=SMOOTHE(oima,gmsmo)
   mge=ois lt gmtho		;outer general mask
   mpe=ois gt gmthi		;inner auxiliar mask
   mima=(ima>0)*mge*mpe		;intermediate masked image
   mea=mean(mima(where(mima)))	;MEA - mean of difference penumbral image

   thr=a+b*mea			; segmentation threshold
   print,i,thr
   ima=float(ima) gt thr	; thresholding mask
   ima=ima*oima			; masked original image

   mpe=ima gt pethr	; mask to remove too faint BFs
   ima=ima*mpe		; too faint features eliminated
   ima=ima*mge*msk			; resulting masked image

   movout(0,0,i)=ima(32:575,66:433)	; compose resulting movie

endfor

end