pro tracing_auto2,image1,image2,loopfile,para,output,test
;+
; Project     : SOHO, TRACE, STEREO, HINODE, SDO
;
; Name        : TRACING_AUTO1
;
; Category    : Automated 2D Pattern Recognition
;
; Explanation : An input image (IMAGE1) is highpass-filtered (IMAGE2)
;		and automatically traced for curvi-linear structures,
;		whose [x,y] pixel coordinates are stored in LOOPFILE.
;		The algorithm iteratively starts to trace a loop from
;		the position of the absolute flux maximum in the image,
;		tracks in a bi-directional way by oriented directivity,
;		based on the direction of the ridge with maximum flux.
;
; Syntax      : IDL> tracing_auto2,image1,image2,loopfile,para,test
;
; Inputs      : image1   - input image 
;               para     - [nsm1,nsm2,nstruc,rmin,qfill,lmin,qmed]
;			   contains algorithm control parameters, such as:
;			     nsm1  = lowpass filter [pixels]
;			     nsm2  = highpass filter [pixels]
;			     nstruc= maximum number of analyzed structures
;			     rmin  = minimum curvature radius [pixels]
;		             qfill =minimum filling ratio (flux_avg/flux_peak)
;			     lmin  = minimum length of loop segments
;			     qmed  = ratio of image noise level
;               test     - flag for no test (0) or interactive test output (1)
;
; Outputs     : image2	 - output image of highpass-filtered image
;               loopfile - output file containing loop coordinates
;                          with consecutive loop numbering
;		output	 - [wid,nlen,fluxmin,fluxmax,fmin,qarea]
;	                      wid  = half widht of loop [pixel]	
;			      nlen = length of moving box [pixel]
;			      fluxmin = minimum flux in image 
;			      fluxmax = maximum flux in image 
;			      fmin    = noise flux level in image 
;			      qarea   = effective image area above base flux
;
; Procedures called inside this routine:
;		TRACING_DIRECTION.PRO
;		TRACING_STEP.PRO
;
; History     : 16-Oct-2009, TRACING_AUTO.PRO written by Markus J. Aschwanden
;	         1-May-2013, noise threshold statistics in 8x8 macropixels
;	         1-May-2013, minimum flux correction (QMED)
;		 8-May-2013, code OCCULT-2 released 
;
; Contact     : aschwanden@sag.lmsal.com
;-

;UNPACK CONTROL PARAMETERS________________________________________________
nsm1	=para(0)
nsm2	=para(1)
nstruc	=para(2)
rmin	=para(3)
qfill 	=para(4)
lmin	=para[5]
qmed	=para(6)
wid 	=long((nsm2-3)/2)
nlen	=long(sqrt(8*rmin*wid))
thresh0 =0.001
reso	=1
ang	=30.	;deg

;FLUX THRESHOLD___________________________________________________________
dim	=size(image1)
nx	=dim(1)
ny	=dim(2)
fluxmin =min(image1)
fluxmax =max(image1)
fluxmed =median(image1)
fmin	=qmed*fluxmed 
ind	=where(image1 lt fmin,nind)
if (nind ge 1) then image1(ind)=fmin
qarea	=(1.-float(nind)/float(nx*ny))
print,'Min/max flux in image = ',fluxmin,fluxmax
print,'Noise level  in image = ',qmed,fmin
print,'Image area      qarea = ',100.*qarea,' %'

;HIGHPASS FILTER___________________________________________________________
if (nsm1 le 2) then image2=image1-smooth(image1,nsm2) 
if (nsm1 ge 3) then image2=smooth(image1,nsm1)-smooth(image1,nsm2) 

;LOOP TRACING______________________________________________________________
openw,2,loopfile
close,2
count	=0
residual=image2
for nloop=0,nstruc-1 do begin
 zmax	=max(residual,im) 
 if (nloop mod 100) eq 0 then begin
  print,'Struc#'+string(nloop,'(I7)')+' Loop#'+string(count,'(I5)')+$
   ' Flux='+string(zmax,'(f10.2)')
 endif
 jstart	=long(im/nx)
 istart	=long(im mod nx)
 tracing_step,residual,istart,jstart,nlen,ang,thresh0,qfill,test,xloop,yloop,zloop,wid,count
 np	=n_elements(xloop)
 s	=fltarr(np)
 looplen=0
 if (np ge 2) then begin
  for ip=1,np-1 do s(ip)=s(ip-1)+sqrt((xloop(ip)-xloop(ip-1))^2+(yloop(ip)-yloop(ip-1))^2)
  looplen=s(np-1)
 endif

;ELIMINATE COLUMN AND ROW ARTIFACTS________________________________________
  statistic,xloop,xav,xdev
  statistic,yloop,yav,ydev
  if (xdev lt 1) then looplen=0    ;vertical column artifact
  if (ydev lt 1) then looplen=0    ;horizontal row  artifact

;STORE LOOP COORDINATES____________________________________________________
  if (looplen ge lmin) then begin
   ns	 =long(looplen)>3
   ss	 =findgen(ns)
   nn    =long(ns/reso+0.5)
   ii    =findgen(nn)*reso
   xx    =interpol(xloop,s,ii)                           ;interpolate
   yy    =interpol(yloop,s,ii)
   ff    =interpol(zloop,s,ii)                           ;flux average
   openw,2,loopfile,/append
   for ip=0,nn-1 do printf,2,count,xx(ip),yy(ip),ff(ip),ss(ip)
   close,2
   if (zmax le thresh0) then goto,end_trace
   count=count+1
  endif				  ;if (looplen ge lmin) then begin

;ERASE LOOP IN RESIDUAL IMAGE____________________________________________
 ns	=n_elements(xloop) > 1
 for is=0,ns-1 do begin
  i0    =(long(xloop(is))>0)<(nx-1)
  j0    =(long(yloop(is))>0)<(ny-1)
  i3    =long(i0-wid)>0
  j3    =long(j0-wid)>0
  i4    =long(i0+wid)<(nx-1)
  j4    =long(j0+wid)<(ny-1)
  residual(i3:i4,j3:j4)=0.
 endfor

endfor
end_trace:
output	=[wid,nlen,fluxmin,fluxmax,fmin,qarea]
end