pro ift2,movin,outp
;
; INDIVIDUAL FEATURE TRACKING ver.2
;  The procedure tracks the motion of a bright or dark feature
;  (selected in the first frame by a cursor) during a time-series
;  of frames. Position (w.r. to the whole frame) and intensity of 
;  the maximum (minimum) are measured in each frame inside a subimage
;  of the tracking field (position interpolated with fivepoin).
;  The tracking offsets are computed with subpixel accuracy (ALIGN_C2)
;  and position of the tracking field and the reference frame are
;  refreshed when twice happens that modulus[offset] > 0.75 pix.
;
; Input : 3-d array MOVIN containing a time-series of frames
; Output: 2-d array OUTP: Frame.no.; X-position; Y-position; Intensity
; Calls : align_c2, fivepoin, mean, maxloc, nint, animate
;
; 2 June 1995, Michal
;
;------------------ TUNABLE PARAMETERS --------------------------------------
fld= 48 ; 24 for 0.06"/pix, 12 for 0.12"/pix, maximum 48
	;Size of the tracking field around the feature (pix). 1.5" recommended
ww = 19  ;  7 for 0.06"/pix,  5 for 0.12"/pix  - MUST BE ODD
        ;Size of the field where to find the extremum (pix).  0.5" recommended
;----------------------------------------------------------------------------

dim=size(movin)
if (dim(0) ne 3) then begin
   print,'Input must be a 3-D array (movie). STOPPED, type RETALL'
   stop
endif

outp=fltarr(4,dim(3))		; output array
chmov=intarr(fld,fld,dim(3))	; check movie
ld=fld/2			; defines the tracking field
ud=fld/2-1		        ; (ud=fld/2-1  for even value of fld)
w=ww/2
fact=float(fld)/48.		; resampling factor (align_c2)

window,0,xsize=dim(1),ysize=dim(2)
tvscl,movin(*,*,0)				; first frame
curs: print,'Select single feature with cursor'
      cursor,x,y,/device,/down			; get initial position

      read,'The selected feature is DARK (0), BRIGHT (1) or wrong (2) ?',mami
      if (mami gt 1 or mami lt 0) then goto,curs           ; repeat selection

if (mami eq 0) then begin			; dark - looking for minimum
      val=min(movin(x-w:x+w,y-w:y+w,0),pos)	;  in ww x ww environment
endif else begin				; bright-looking for maximum
      val=max(movin(x-w:x+w,y-w:y+w,0),pos)	;  in ww x ww environment
endelse
xm=(pos mod ww)+x-w				; rough position-
ym=pos/ww+y-w					;  -in the whole frame
cc=movin(xm-1:xm+1,ym-1:ym+1,0)			; 3x3 environment
fivepoin,cc,dx,dy
xm=xm+dx   &   ym=ym+dy				; interpolated position

outp(0,0)=[0,xm,ym,val]				; 1st row of output
print,'Initial position:',x,y,'      Intensity:',movin(x,y,0)
print,'Feature position:',xm,ym,'    Intensity:',val

x=nint(xm)   &   y=nint(ym)			; automatic centering
chmov(0,0,0)=movin(x-ld:x+ud,y-ld:y+ud,0)	; initial subframe

   iref=0				 ; 1st frame for reference
   ya=0					 ; counter of cases

for i=1,dim(3)-1 do begin  ;---------------------loop-------------------------

   a=movin(x-ld:x+ud,y-ld:y+ud,iref)	 ; reference
   b=movin(x-ld:x+ud,y-ld:y+ud,i)	 ; current tracking field
   ref=congrid(a,48,48,/cubic)		 ; resampling for align_c2
   ima=congrid(b,48,48,/cubic)

					    ; tracking window 20x20 pix,
   offset=ALIGN_C2(ref,ima,flg,20,20)*fact  ;   max.offset (pixels) =
					    ; = (32-12)/2*fact = 10*fact (old)
					    ; = (48-20)/2*fact = 14*fact
   
   print,'frame',i,'   ref',iref,'   offset:',offset
   if(flg eq 1) then begin
      print,'*** TRACKING LOST ***'
      goto,fin
   endif

   ofdst=sqrt(offset(0)*offset(0) + offset(1)*offset(1))
   if(ofdst gt 0.75) then begin

     if(ya lt 1) then begin			;one case could be random-
       ya=ya+1					; -wait for the 2nd case
     endif else begin				;significant case already
       x=x-nint(offset(0)) & y=y-nint(offset(1)) 	; new position
       iref=i						; new reference
       ya=0
     endelse

   endif

   chmov(0,0,i)=movin(x-ld:x+ud,y-ld:y+ud,i)

   if (mami eq 0) then begin			; dark - looking for minimum
      val=min(movin(x-w:x+w,y-w:y+w,i),pos)	;  in ww x ww environment
   endif else begin				; bright-looking for maximum
      val=max(movin(x-w:x+w,y-w:y+w,i),pos)	;  in ww x ww environment
   endelse
   xm=(pos mod ww)+x-w				; rough position-
   ym=pos/ww+y-w				;  -in the whole frame
   cc=movin(xm-1:xm+1,ym-1:ym+1,i)		; 3x3 environment
   fivepoin,cc,dx,dy
   xm=xm+dx   &   ym=ym+dy			; interpolated position

outp(0,i)=[i,xm,ym,val]				; output (val not interpolated)

endfor	;---------------------end of loop--------------------------------------

fin:
outp=outp(*,0:i-1)			; truncation of empty elements
chmov=chmov(*,*,0:i-1)			;(i-1 is the last processed frame)

; presentation

window,0,xsiz=740,ysiz=720
x1=min(outp(1,*))-5 & x2=max(outp(1,*))+5	; plot the path
y1=min(outp(2,*))-5 & y2=max(outp(2,*))+5
dxy=(x2-x1) > (y2-y1)  ;to preserve the aspect ratio
plot,outp(1,*),outp(2,*),/nodata,back=255,col=0,xra=[x1,x1+dxy],yra=[y1,y1+dxy],/xst,/yst
for j=1,i-1 do arrow,outp(1,j-1),outp(2,j-1),outp(1,j),outp(2,j),col=0,/data

chmov=bytscl(rebin(chmov,fld*5,fld*5,i))
fra=bytarr(fld*5,fld*5)				; mark the search subfield...
x1=5*(ld-w)  &  x2=5*(ld+w)
y1=x1        &  y2=x2
fra(x1,*)=130
fra(x2,*)=130
fra(*,y1)=130
fra(*,y2)=130
for k=0,i-1 do chmov(*,*,k)=byte(chmov(*,*,k)+fra)
animate,chmov					; show it in movie

end