pro nlfff_tracing_test,dir_tmp,dir,filename,vers,para,test,wave_,nr,diffract_
;+
; Project     : SDO, STEREO, SOHO
;
; Name        : NLFFF_TRACING
;
; Category    : data analysis and automated feature recognition
;
; Explanation : Reads and EUV image, enables automated loop tracing,
;               and writes loop coordinates into loopfile *loop*.dat.
;
; Syntax      : IDL>nlfff_tracing,dir_tmp,dir,filename,vers,para,test,wave_,nr
;
; Inputs      : dir_tmp,dir,filename,vers,para,test,wave_,nr
;
; Outputs     ; nr	;total number of loops in all wavelengths
;               saves loop coordinates in loopfile : *loop*.dat
;
; History     : 21-Sep-2012, Version 1 written by Markus J. Aschwanden
;	         9-May-2014, add mask of forbidden areas (LOOPTRACING_AUTO4.PRO)
;
; Contact     : aschwanden@lmsal.com
;-

print,'__________________NLFFF TRACING_____________________________
nmag0   =long(para[1])
nsm1	=long(para[4])
nsmax	=long(para[7])    
dfoot	=para[ 2]
dsmag	=para[ 3]  	  
qsig	=para[ 5]
nsig_d	=para[ 6]  	  ;significance level of diffraction pattern
ds	=para[13]
rmin	=para[15]
lmin	=para[16]	  ;in pixels 
x1      =para[17]	  ;FOV in solar radii
y1      =para[18]
x2      =para[19]
y2      =para[20]
nloopw	=para[21]
ngap	=para[22]
thresh2	=para[27]
dprox	=para[28]
qripple	=para[29]
nsig	=0.5    	  ;display contrast of highpass filter
nsm2	=nsm1+2
fsat	=2.^14-1 
black	=0
blue	=50
purple	=75
red	=110
orange	=150
yellow	=200
white	=255
char	=1
autofile  =dir_tmp+filename+'_auto_'+vers+'.tmp'
coeff_file=dir+filename+'_coeff_'+vers+'.dat'
loopfile  ='loop_dump.dat'
savefile  =dir+filename+'_field_'+vers+'.sav'
restore,savefile

;________________WAVELENGTH LOOP_____________________________
nr	=0
nwave	=n_elements(wave_)
openw,1,loopfile
for iwave=0,nwave-1 do begin
 wave	=wave_[iwave]
 filename_w=filename+'_'+wave_(iwave)
 imagefile =dir_tmp+filename_w+'_euv.fits'
 filterfile=dir+filename_w+'_euv_fov.fits'
 n_sat	=0
 n_bleed=0
 n_prox =0
 n_ripple=0

;print,'__________READ EUV IMAGE__________________________________
 data	=readfits(imagefile,h,/silent)
 dateobs=sxpar(h,'DATE_OBS')
 naxis1	=sxpar(h,'NAXIS1')
 naxis2	=sxpar(h,'NAXIS2')
 cdelt1 =sxpar(h,'CDELT1')
 cdelt2	=sxpar(h,'CDELT2')
 crpix1	=sxpar(h,'CRPIX1')
 crpix2	=sxpar(h,'CRPIX2')
 eph	=get_sun(dateobs)
 rsun 	=eph(1)			;rsun=index.rsun_obs in arcseconds 
 rpix    =rsun/cdelt1
 dpix	=1./rpix
 dx_euv  =cdelt1/rsun		;pixel size in solar radii
 i1      =long(crpix1+x1/dpix+0.5)
 i2      =long(crpix1+x2/dpix+0.5)
 j1      =long(crpix2+y1/dpix+0.5)
 j2      =long(crpix2+y2/dpix+0.5)
 nx      =(i2-i1+1)
 ny      =(j2-j1+1)
 image1  =float(data(i1:i2,j1:j2))

;_________________BYPASS-FILTERED IMAGE____________________________
 if (nsm1 le 2) then image2=image1-smooth(image1,nsm2)
 if (nsm1 ge 3) then image2=smooth(image1,nsm1)-smooth(image1,nsm2)
 writefits,filterfile,image2 

;_________________DISPLAY EUV IMAGE________________________________
 zoom    =long(1024/nx)>1
 window,0,xsize=nx*zoom,ysize=ny*zoom
 !p.position=[0,0,1,1]
 !x.range=[x1,x2]		;FOV coordinates
 !y.range=[y1,y2]
 !x.style=1
 !y.style=1
 loadct,0,/silent
 plot,[0,0] 
 statistic,image2,zav,zsig
 c1	=zav-nsig*zsig
 c2	=zav+nsig*zsig
 tv,bytscl(rebin(image2,nx*zoom,ny*zoom),min=c1,max=c2)
 loadct,5,/silent

;__________________BACKGROUND SUBIMAGE WITH LEAST NOISE_____________
 thresh =max(image1)
 nbin	=10
 nxx	=long(nx/nbin)
 nyy	=long(ny/nbin)
 for ii=0,nbin-1 do begin
  ib1=ii*nbin
  ib2=(ii+1)*nbin-1
  for jj=0,nbin-1 do begin
   jb1=jj*nbin
   jb2=(jj+1)*nbin-1
   subimage=image2[ib1:ib2,jb1:jb2]
   statistic,subimage,zavg,zsig
   noise=zavg+qsig*zsig
   if (noise lt thresh) then thresh=noise
  endfor
 endfor

;_______________________MAGNETIC SOURCES___________________________
 readcol,coeff_file,mm1,xm1,ym1,zm1,am1,/silent
 coeff=[[mm1],[xm1],[ym1],[zm1],[am1]]
 nmag	=n_elements(mm1)

;_______________________AUTOMATED LOOP TRACING_____________________
 mask	 =fltarr(nx,ny)
 x	 =x1+(x2-x1)*findgen(nx)/float(nx-1)
 y	 =y1+(y2-y1)*findgen(ny)/float(ny-1)
 dx	 =x(1)-x(0)
 dy	 =y(1)-y(0)
 for im=0,nmag-1 do begin
  i3	 =(long((xm1[im]-dfoot-x1)/dx+0.5) > 0) < (nx-1)
  i4	 =(long((xm1[im]+dfoot-x1)/dx+0.5) > 0) < (nx-1)
  j3	 =(long((ym1[im]-dfoot-y1)/dx+0.5) > 0) < (ny-1)
  j4	 =(long((ym1[im]+dfoot-y1)/dx+0.5) > 0) < (ny-1)
  if (i4 gt i3) and (j4 gt j3) then $
  mask(i3:i4,j3:j4)=1			;valid area for loop initialization
 endfor
 nstruc	 =10000
 para2	 =[nsm1,rmin,lmin,nstruc,thresh,nloopw,ngap,thresh2]
 if (test eq -1) then test1=nstruc+1
 if (test eq  0) then test1=0
 looptracing_auto4,image1,image2,autofile,para2,output,test1,mask

;_______________________COORDINATE TRANSFORMATION__________________
 file	=findfile(autofile,count=nfiles)
 if (nfiles le 0) then goto,end_wave
 readcol,autofile,iloop,ix,iy,flux,is,/silent
 il1	=long(min(iloop))
 il2	=long(max(iloop))
 nloop	=il2-il1+1
 for il=il1,il2 do begin
  ind	=where(iloop eq il,nind)
  if (nind eq 0) then goto,skip_loop
  ii	=ix(ind)
  jj	=iy(ind)
  x	=(i1-crpix1+ii)*dpix
  y	=(j1-crpix2+jj)*dpix
  f	=flux(ind)
  s	=is(ind)/rpix
  len	=0.
  for is=1,nind-1 do len=len+sqrt((x(is)-x(is-1))^2+(y(is)-y(is-1))^2)
  bad	=0
;......................check proximity to magnetic source........
  z     =sqrt(1.-x^2-y^2)
  n	=nind-1
  dist1	=min(sqrt((xm1-x[0])^2+(ym1-y[0])^2+(zm1-z[0])^2),id1)
  dist2	=min(sqrt((xm1-x[n])^2+(ym1-y[n])^2+(zm1-z[n])^2),id2)
  if (dist1 gt dprox) and (dist2 gt dprox) then begin
   bad	=1
   n_prox=n_prox+1
  endif
;......................check saturated pixels ...................
  ii1	=(ii-nsm1) > 0
  ii2	=(ii+nsm1) < (nx-1)
  jj1	=(jj-nsm1) > 0
  jj2	=(jj+nsm1) < (ny-1)
  dn1   =max(image1(ii1,jj))	;nsm1 pixels to the north of loop
  dn2   =max(image1(ii2,jj))	;nsm1 pixels to the south of loop
  dn3   =max(image1(ii,jj1))	;nsm1 pixels to the east  of loop
  dn4   =max(image1(ii,jj2))	;nsm1 pixels to the west  of loop
  if (dn1 ge fsat) or (dn2 ge fsat) or (dn3 ge fsat) or (dn4 ge fsat) then begin
   bad  =2
   n_sat =n_sat+1
  endif
;......................check vertical pixel bleeding.............
  dxx	=abs(x[1:nind-1]-x[0:nind-2])
  ivert =where(dxx le ds*0.01,nvert)
  if (nvert ge nind/2) then begin
   bad=3
   n_bleed=n_bleed+1
  endif
;......................flux ripple................................
  f1	=f(0:nind-2)
  f2	=f(1:nind-1)
  ind0	=where((f1 gt 0) or (f2 gt 0),nind0)
  df	=abs(f2-f1)
  fm	=f1 > f2
  qf	=total(df(ind0)/fm(ind0))/float(nind0)
  if (qf ge qripple) and (len le (2.*lmin*dx_euv)) then begin
   bad =4
   n_ripple=n_ripple+1
  endif
;.................................................................
  for ii=0,nind-1 do printf,1,nr,x(ii),y(ii),f(ii),s(ii),wave,bad,$
     format='(I6,2F12.6,F12.2,F10.6,2I6)'
  if (bad eq 0) then oplot,x,y,color=yellow,thick=3
  if (bad eq 1) then oplot,x,y,color=purple,thick=3
  if (bad eq 2) then oplot,x,y,color=red,thick=3
  if (bad eq 3) then oplot,x,y,color=orange,thick=3
  if (bad eq 4) then oplot,x,y,color=blue,thick=3
  nr	=nr+1
SKIP_LOOP:

 endfor
 print,wave,' A : sat,bleed,prox,ripple,nloop =',$
  n_sat,n_bleed,n_prox,n_ripple,nloop,format='(2a,5i6)'
 wait,1

;read,'continue?',yes
;if (yes eq 0) then stop

 END_WAVE:
endfor			;wavelength loop
close,1

;________________READ LOOP FILE____________________________
readcol,loopfile,iloop,xl,yl,fl,sl,wl,bl,/silent
nf	=long(max(iloop)+1)

;________________MAXIMUM LENGTH OF FIELDLINES______________
ns_loopmax=0
ns_loop =lonarr(nf)
for k=0,nf-1 do begin
 ind	=where(iloop eq k,ns)
 ns_loop[k]=ns
 ns_loopmax=ns_loopmax > ns
endfor

;________________DISPLAY MAGNETOGAM__________________________
dim	=size(bzmap)
nx	=dim(1)
ny	=dim(2)
zoom	=long(1024/nx)>1
window,0,xsize=nx*zoom,ysize=ny*zoom
loadct,0,/silent
!p.position=[0,0,1,1]
plot,[0,0],[0,0],xrange=[x1,x2],yrange=[y1,y2],xstyle=1,ystyle=1
statistic,bzmap,bz_avg,bz_sig
c1	=bz_avg-3*bz_sig
c2	=bz_avg+3*bz_sig
tv,bytscl(rebin(bzmap,nx*zoom,ny*zoom),min=c1,max=c2)

;________________FIELD LINE ARRAY___________________________
loadct,5,/silent
field	=fltarr(nsmax,3,nf)
icount	=0
dire	=fltarr(nf)
dire1	=fltarr(nf)
dire2	=fltarr(nf)
curv	=fltarr(nf)
len	=fltarr(nf)
for k=0,nf-1 do begin
 ind	=where(iloop eq k,ns)
 ns	=ns < nsmax
 if (ns ge 1) then begin
  xs	=xl(ind[0:ns-1])
  ys	=yl(ind[0:ns-1])
  fs	=fl(ind[0:ns-1])
  field(0:ns-1,0,icount)=xs 		;x-coordinate loop
  field(0:ns-1,1,icount)=ys		;y-coordinate loop
  field(0:ns-1,2,icount)=fs		;flux of loop
;..............................................................
  c	=linfit(xs,ys)
  c1	=linfit(xs(0:ns/2),ys(0:ns/2))
  c2	=linfit(xs(ns/2:ns-1),ys(ns/2:ns-1))
  dire(k)=atan(c(1))*(180./!pi)		;direction in degrees
  dire1(k)=atan(c1(1))*(180./!pi)	;direction in degrees
  dire2(k)=atan(c2(1))*(180./!pi)	;direction in degrees
  curv(k)=abs(dire2(k)-dire1(k))
  for is=1,ns-1 do $
  len(k)=len(k)+sqrt((xs(is)-xs(is-1))^2+(ys(is)-ys(is-1))^2)  ;solar radii
;..............................................................
  ns_loop[icount]=ns				
  icount=icount+1
 endif
endfor

;_______________histogram peak of direction____________________
bin_angle=10.
nbin	=long(180./bin_angle)
yh     =histogram(dire,binsize=bin_angle,min=-90.,max=90.)
xh     =-90.+bin_angle*(findgen(nbin)+0.5)
nloop_peak=max(yh,ipeak)
angle_peak=xh(ipeak)
nloop_peak=yh(ipeak)
noise_peak=sqrt(nloop_peak)
nloop_med =median(yh)
nsig_peak =(nloop_peak-nloop_med)/noise_peak
d1	=angle_peak-bin_angle
d2	=angle_peak+bin_angle
if (angle_peak gt 0) then begin &d3=d1-90. &d4=d2-90. &endif
if (angle_peak le 0) then begin &d3=d1+90. &d4=d2+90. &endif
print,'Significance diffraction pattern = ',nsig_peak,nsig_d

;_______________Selecting good loops__________________________
badloop =fltarr(nf)
for k=0,nf-1 do begin
 ind	=where(iloop eq k,ns)
 ns	=ns < nsmax
 if (ns ge 1) then begin
  xs	=xl(ind[0:ns-1])
  ys	=yl(ind[0:ns-1])
  zf	=fl(ind[0:ns-1])
  badloop(k)=median(bl(ind))	;previous bad values bad=0,1,2,3,4
  if (badloop(k) eq 0) then oplot,xs,ys,color=yellow,thick=3
  if (badloop(k) eq 1) then oplot,xs,ys,color=purple,thick=3
  if (badloop(k) eq 2) then oplot,xs,ys,color=red,thick=3
  if (badloop(k) eq 3) then oplot,xs,ys,color=orange,thick=3
  if (badloop(k) eq 4) then oplot,xs,ys,color=blue,thick=3

  NEXT_SOURCE:
  dire_sel=0
  curv_sel=0
  if (dire(k) ge d1) and (dire(k) le d2) then dire_sel=1   ;loops parallel to peak dir
  if (dire(k) ge d3) and (dire(k) le d4) then dire_sel=1   ;loops orthogonal to peak dir
  if (curv(k) le bin_angle) then curv_sel=1		   ;straight loops 
  if (dire_sel eq 1) and (curv_sel eq 1) and (nsig_peak ge nsig_d) then begin
   bl(ind)=5
   badloop(k)=5
   oplot,xs,ys,color=black,thick=3
  endif

 endif
 oplot,xm1,ym1,psym=1,symsize=5
endfor

ind_bad =where(badloop ne  0,nind_bad)
ind_good=where(badloop eq  0,nind_good)
ind1	=where(badloop eq  1,n1)
ind2	=where(badloop eq  2,n2)
ind3	=where(badloop eq  3,n3)
ind4	=where(badloop eq  4,n4)
ind5	=where(badloop eq  5,n5)

;final field length array.....................................
field_loops=field(*,*,ind_good)
ns_loop    =ns_loop[ind_good]

print,'Number of loops            NF    = ',nf 
print,'Number of good loops       NGOOD = ',nind_good 
print,'Number of bad loops        NBAD  = ',nind_bad 
print,'- magnetic proximity       N1    = ',n1
print,'- saturated pixels         N2	= ',n2
print,'- pixel bleeding           N3	= ',n3
print,'- rippled flux profile     N4	= ',n4
print,'- diffraction pattern      N5	= ',n5
features=[nf,nind_good,nind_bad,n1,n2,n3,n4,n5]

;______________DISPLAY NLFFF SOURCES______________________________
for im=0,nmag-1 do begin
 mm     =coeff(im,0)
 xm     =coeff(im,1)
 ym     =coeff(im,2)
 zm     =coeff(im,3)
 rm	=sqrt(xm^2+ym^2+zm^2)
 dm 	=1.-rm
 phi    =2.*!pi*findgen(101)/100.
 if (mm lt 0) then col=255
 if (mm gt 0) then col=0
 size	=char*3
 oplot,[xm,xm],[ym,ym],psym=1,color=col
 oplot,xm+dm*cos(phi),ym+dm*sin(phi),thick=3,color=col
 xyouts,xm,ym,string(im,'(i4)'),size=size,color=col
endfor
end