;Program RECUR,datin,surr,gap,datout,newlf
;
;RECUR - a revision of the results of OTRACK2 to eliminate
;       gaps (due to bad seeing) in the tracking series of objects
;	(if the interruption is shorter than the interval specified in GAP,
;	the object is considered to be the same).
;NOTE:  Pre-processing of DATIN by OTRIM to set lower limit for lifetime
;       and upper limit to velocity is recommended.
;
;INPUTS:  DATIN - 3-D data volume created by OTRACK2
;	   ([intens., pos., size] * frames * objects)
;	  SURR -  linear size of surroundings where to look for
;		 the successive object (see comments in SUROUND)
;	  GAP  -  maximum allowed interval for missing
;		 frames in one object (= max.no. of missing frames)
;OUTPUTS: DATOUT - 3-D volume with revised data (objects rearranged)
;	  NEWLF  - 1-D array with new lifetimes
;
;This file contains two modules: RECUR and the function SUROUND
;Do not forget to set the 'sx' parameter in module RECUR.
;
; 9 July 1996, Michal.

;-----------------------------------------------------------------------

FUNCTION SUROUND,pos1,pos2,surr,sx

;A subroutine for RECUR.
;POS1, POS2 are positions in a 2-D frame with x-dim=SX
;in the 1-D representation.
;SURR defines a size of surroundings of POS1
;(if SURR=1 then surroundings is 3x3, if SURR=2 then 5x5 etc.)
;If POS2 is in the surroundings of POS1 then SURROUND returns 1,
;otherwise it returns 0 (byte).

rrus=surr*(-1)

for j=rrus,surr do begin
  for i=rrus,surr do begin
    if (pos1+j*sx+i) eq pos2 then RETURN,1B
  endfor
endfor

RETURN,0B
END

;-----------------------------------------------------------------------


PRO RECUR,datin,surr,gap,datout,newlf
;
;*********** PARAMETERS *******************************
sx = 70		;x-size of frames processed by OTRACK2
;******************************************************


s=size(datin)	;s(2) - number of frames, s(3) - number of objects

datout=datin
newlf=intarr(s(3))

hi=intarr(5,s(3))	;array of history of objects: 0 - no.,
			;1 - frame of birth, 2 - pos. of birth,
			;3 - frame of end, 4 - pos. of end.

for i=0,s(3)-1 do begin		;loop over objects
  w=where(datin(0,*,i),cnt)	;cnt frames where object exists
  hi(0,i)=i			;object number
  hi(1,i)=w(0)			;1st frame of existence
  hi(2,i)=datin(1,w(0),i)	;position in 1st frame
  hi(3,i)=w(0)+cnt-1		;last frame of existence
  hi(4,i)=datin(1,w(0)+cnt-1,i)	;position in last frame
endfor

;Elimination of gaps

pair=intarr(4,s(3)/2)		;list of objects to be linked
k=0

for i=0,s(3)-1 do begin		;main loop over objects
  wt=where((hi(1,*) gt hi(3,i)) and (hi(1,*)-hi(3,i) le gap+1),cnt)

  if cnt le 0 then goto,endloop		;check time difference
  for j=0,cnt-1 do begin
    if SUROUND(hi(4,i),hi(2,wt(j)),surr,sx) eq 1 then begin
      pair(0,k)=hi(0,i)		;1st object no.
      pair(1,k)=hi(0,wt(j))	;2nd object no.
      pair(2,k)=hi(3,i)		;last frame of 1st object
      pair(3,k)=hi(1,wt(j))	;first frame of 2nd object
      k=k+1
      goto,endloop	;WEAKPOINT - if these conditions are met
			;by more objects than 1 (ambiguous link -
			; - this should not happen), only the 1st
			;one is taken without checking the others.
    endif
  endfor

endloop:
endfor

pair=pair(*,0:k-1)		;length trim

;Linking the objects

openw,1,'links.dat'
for i=k-1,0,-1 do begin	;reverse order to perform multiple links

  datout(*,*,pair(0,i))=datout(*,*,pair(0,i))+datout(*,*,pair(1,i))
  ;in datout, the data are zero if the object does not exist -
  ; - i.e. we can link objects adding their datout sections

  datout(*,*,pair(1,i))=0	;destroy 2nd object

  ;fill the gap by replication of 1st object last position
  ; and by average values of brightness and size
  br=(datout(0,pair(2,i),pair(0,i))+datout(0,pair(3,i),pair(0,i)))/2
  sz=(datout(2,pair(2,i),pair(0,i))+datout(2,pair(3,i),pair(0,i)))/2
    for j=pair(2,i)+1,pair(3,i)-1 do begin
      datout(1,j,pair(0,i))=datout(1,pair(2,i),pair(0,i))
      datout(0,j,pair(0,i))=br
      datout(2,j,pair(0,i))=sz
    endfor

  print,'Link',i,':',pair(0,i),' with',pair(1,i),' fr.',pair(2,i),pair(3,i)
  printf,1,i,':',pair(0,i),' -',pair(1,i),' fr.',pair(2,i),pair(3,i)

endfor
close,1

;Get newlf and restrict datout to existing objects.
for i=0,s(3)-1 do begin
  w=where(datout(0,*,i),cnt)
  newlf(i)=cnt
endfor
datout=datout(*,*,where(newlf))
newlf=newlf(where(newlf))

print,k,' links'
sl=fix(n_elements(newlf))
print,fix(s(3)),' objects originally, after linking',sl

END