;------------------------------------------------------------------------
;  (C) Alexander Heger, August 2000 
;
;  IDEA: 
;  remove datapoints that result into ploting below the
;  specified resolution.
;  
;  HOW IT WORKS: 
;  From a starting point determine opening angle under which an
;  ellipse with major axis of the given resolution appears.  Compute
;  Intersection of this opening angle and the previous.  Keep only
;  most remote point as long as there is an overlap of the opening
;  angles (i.e., remove intermediate points that lay on a streight
;  line within the given resolution).
;  
;  Modifications:
;  30 Dec 2000: rewritten
;  31 Dec 2000: debugged (variable d2 introduced)
;  01 Jan 2001: parameter PLOT added
;  14 Aug 2001: support for XY data in one array
;  04 Jun 2003: TRUNCATE parameter added
;------------------------------------------------------------------------

PRO reduce,xv,yv,XLOG=xlog,YLOG=ylog,XRES=xres,YRES=yres, $
           NREM=nrem, RES=res, $
           XRANGE=xrange, YRANGE=yrange, $
           RELRES=relres, XRELRES=xrelres, YRELRES=yrelres, $
           METHOD=method, MAP=map, PLOT=plot, $
           VERBOSE=verbose, TRUNCATE=truncate

IF N_ELEMENTS(verbose) EQ 0 THEN verbose=0


nnx0=N_ELEMENTS(xv)
nny0=N_ELEMENTS(yv)

; if [x|y]log then res==d[x|y]/[x|y] ELSE res==d[x|y]

array_mode=0
nndx=SIZE(xv,/N_DIMENSIONS)
ndx=SIZE(xv,/DIMENSIONS)
IF nndx EQ 2 THEN BEGIN
    IF N_ELEMENTS(yv) NE 0 THEN BEGIN
        PRINT,' [REDUCE] ERROR: accept only one muli-D array'
        RETURN
    ENDIF
    IF ndx[0] EQ 2 THEN BEGIN
        array_mode=1
        yv=xv[1,*]
        xv=xv[0,*]
    ENDIF ELSE IF ndx[1] EQ 2 THEN BEGIN
        array_mode=2
        yv=xv[*,1]
        xv=xv[*,0]
    ENDIF
ENDIF

IF verbose NE 0 THEN PRINT,' [REDUCE] array_mode:'+STRING(array_mode)

n=N_ELEMENTS(xv)
IF n_ELEMENTS(yv) NE n THEN BEGIN
    PRINT,' [REDUCE] dim error'
    RETURN
ENDIF

IF N_ELEMENTS(xlog) EQ 0 THEN xlog=0
IF N_ELEMENTS(ylog) EQ 0 THEN ylog=0

IF N_ELEMENTS(plot) EQ 0 THEN plot=0
IF plot EQ 1 THEN BEGIN
    IF !X.TYPE EQ 0 THEN BEGIN 
        xrange=!X.CRANGE 
    ENDIF ELSE BEGIN 
        xrange=10.D0^!X.CRANGE
        xlog=1
    ENDELSE
    IF !Y.TYPE EQ 0 THEN BEGIN
        yrange=!Y.CRANGE 
    ENDIF ELSE BEGIN
        yrange=10.D0^!Y.CRANGE
        ylog=1
    ENDELSE
ENDIF

IF (N_ELEMENTS(relres) NE 0) THEN BEGIN
    relflagx=1
    relflagy=1
ENDIF

relflag=0
IF N_ELEMENTS(relres)  EQ 0 THEN relres =1.0D-4 ELSE relflag=1

relflagx=relflag
IF N_ELEMENTS(xrelres) EQ 0 THEN xrelres=relres ELSE relflagx=1

relflagy=relflag
IF N_ELEMENTS(yrelres) EQ 0 THEN yrelres=relres ELSE relflagy=1

IF (relflagx EQ 1) AND N_ELEMENTS(xrange) NE 2 THEN BEGIN
    xrange0=MIN(xv,MAX=xrange1)
    xrange=[xrange0,xrange1]
ENDIF 
IF (relflagy EQ 1) AND N_ELEMENTS(yrange) NE 2 THEN BEGIN
    yrange0=MIN(yv,MAX=yrange1)
    yrange=[yrange0,yrange1]
ENDIF 

IF N_ELEMENTS(truncate) EQ 0 THEN truncate=1

IF N_ELEMENTS(xrange) EQ 2 THEN BEGIN
    xrange0=MIN(xrange,MAX=xrange1)
    IF truncate THEN xv=(TEMPORARY(xv)>xrange0)<xrange1
    IF (N_ELEMENTS(res) EQ 0) AND (N_ELEMENTS(xres) EQ 0) THEN BEGIN
        IF xlog NE 0 THEN BEGIN
            xres=(ALOG(xrange1)-ALOG(xrange0))*xrelres
        ENDIF ELSE BEGIN
            xres=(xrange1-xrange0)*xrelres            
        ENDELSE
    ENDIF
ENDIF
IF N_ELEMENTS(yrange) EQ 2 THEN BEGIN 
    yrange0=MIN(yrange,MAX=yrange1)
    IF truncate THEN yv=(TEMPORARY(yv)>yrange0)<yrange1
    IF (N_ELEMENTS(res) EQ 0) AND (N_ELEMENTS(yres) EQ 0) THEN BEGIN
        IF ylog NE 0 THEN BEGIN
            yres=(ALOG(yrange1)-ALOG(yrange0))*yrelres
        ENDIF ELSE BEGIN
            yres=(yrange1-yrange0)*yrelres            
        ENDELSE
    ENDIF
ENDIF

IF N_ELEMENTS(res)  EQ 0 THEN res =1.0D0
IF N_ELEMENTS(xres) EQ 0 THEN xres=res
IF N_ELEMENTS(yres) EQ 0 THEN yres=res

nrem=0L

xresi=1.D0/xres
yresi=1.D0/yres

; some specials for the log treatment
; use working array with log values and scaled by resolution
; use map to select final choice of elements
IF xlog EQ 1 THEN BEGIN
    x=alog(xv)*xresi
ENDIF ELSE BEGIN
    x=xv*xresi
ENDELSE
IF ylog EQ 1 THEN BEGIN
    y=alog(yv)*yresi
ENDIF ELSE BEGIN
    y=yv*yresi
ENDELSE
map=lindgen(n)

IF N_ELEMENTS(method) EQ 0 THEN method=2

IF (method mod 2) EQ 0 THEN GOTO,PART_II
;--------------------
PART_I:
;--------------------

; the obvious: remove dense points
; keep, however, first & last point
; remove point when within a ellipse with major axis xres and yres
i1=0L
n1=0L
map[n1]=map[i1]
FOR i=1L,n-2 DO BEGIN
    dx=x[i]-x[i1]
    dy=y[i]-y[i1]
    IF (dx^2+dy^2) GT 1.D0 THEN BEGIN
        i1=i
        n1=n1+1L
        map[n1]=map[i1]
    ENDIF
ENDFOR
i1=n-1L
n1=n1+1L
map[n1]=map[i1]
nrem=n-n1-1L

; reduce vector length
n=n1+1L
map=map[0:n1]
x=x[map]
y=y[map]

IF (method mod 4)/2 EQ 0 THEN GOTO,FINAL
;--------------------
PART_II:
;--------------------

; the difficult: non-horz./vert. lines
; keep first and last point 
pi=!DPI
pi2=2.D0*pi

n1=0L
i0=0L

i1=1L
i2=1L
d1=0.D0
d2=0.D0
dchi=-1.D0
FOR i=1L,n-1L DO BEGIN
begin_of_loop:
    dx=x[i]-x[i0]
    dy=y[i]-y[i0]
    d=SQRT(dx^2+dy^2)
; d < 1?
    IF (d LE 1.D0) THEN GOTO,remote_point
; d > 1
    psi=ASIN(1.D0/d) ; ==> 0 < psi < pi/2
    phi=ATAN(dy,dx) ; -pi < phi < pi
    phi=(phi+psi+pi2) MOD pi2 ; make sure all angles are positive
    dphi=2.D0*psi ; < pi
    IF (dchi LT 0.0D0) THEN BEGIN
; setup of first point (with d > 1)
        chi=phi
        dchi=dphi
        d1=d
        GOTO,remote_point
    ENDIF ELSE BEGIN
; compute overlap
        xi=(phi-chi+pi2) mod pi2
        IF (xi GT pi) then xi=xi-pi2
        IF xi LT 0.0D0 THEN BEGIN
            chi=phi
            dchi=MIN([dchi+xi,dphi])
        ENDIF ELSE BEGIN
            dchi=MIN([dphi-xi,dchi])
        ENDELSE
    ENDELSE
    IF dchi LT 0.0D0 THEN BEGIN
; point is out of line
        n1=n1+1
        map[n1]=map[i1]
        IF (d1 GT (d2 + 1.0D0)) THEN BEGIN
;        IF (i2 GT i1) THEN BEGIN
            n1=n1+1
            map[n1]=map[i2]
        ENDIF 
; start new segment
        d1=0.0D0
        i0=i-1
        i1=i
        i2=i
        d2=d
        dchi=-1.D0
; redo current point
        GOTO,begin_of_loop
    ENDIF
remote_point:
; save current coordiantes
    i2=i
    d2=d
; check if most remote point in that direction
    IF d GE d1 THEN BEGIN
        i1=i
        d1=d
    ENDIF
ENDFOR
n1=n1+1
map[n1]=map[i1]
IF (i2 GT i1) THEN BEGIN
    n1=n1+1
    map[n1]=map[i2]
ENDIF 
nrem=nrem+n-n1-1

; reduce vector length
map=map[0:n1]
n=n1+1
x=0
y=0

;--------------------
FINAL:
;--------------------

; assign values
xv=xv[map]
yv=yv[map]
; map=0
IF array_mode EQ 1 THEN BEGIN
   xx=xv
   xv=DBLARR(2,n)
   xv[0,*]=TEMPORARY(xx)
   xv[1,*]=TEMPORARY(yv)
ENDIF ELSE IF array_mode EQ 2 THEN BEGIN
   xv=[[TEMPORARY(xv)],[TEMPORARY(yv)]]
ENDIF

IF verbose THEN BEGIN
    nnx1=N_ELEMENTS(xv)
    nny1=N_ELEMENTS(yv)
    PRINT,' [REDUCE] reduction from '+STRING(nnx0+nny0)+' to'+$
      STRING(nnx1+nny1)+' points ('+STRING((1.D0*nnx1+nny1)/(nnx0+nny0))+')'
ENDIF

END




PRO test_reduce

n=10000000L
t=dindgen(n+1)/n

;t=randomu(1,n+1)

;t=t^100

t=t*2*!DPi

;x=exp(cos(3*t)*2)
;y=sin(exp(2.5*t/MAX(t)))

;y=cos(t)
;x=sin(5*t)*exp(5*t/MAX(t))

;y=cos(t*!DPi)
;x=sin(5*t)*exp(5*t/MAX(t))

y=cos(t)^5
x=sin(t)^5

; rotate
phi=!DPI/2
c=COS(phi)
s=SIN(phi)
x1=c*x-s*y
y=s*x+c*y
x=x1

plot,x,y,/NODATA,XMARGIN=[1,1],YMARGIN=[1,1]
oplot,x,y,color=255
xx=[[x],[y]]

;xx=DBLARR(2,n+1)
;xx[0,*]=x
;xx[1,*]=y

;reduce,x,y,relres=1.D-3,NREM=nn,/plot
reduce,xx,relres=1.D-3,NREM=nn,/plot,/verbose
x=xx(*,0)
y=xx(*,1)
;x=xx[0,*]
;y=xx[1,*]

oplot,x,y
psymdot,FILL=0
oplot,x,y,color=255*256L*256L,psym=8
;help,x

PRINT,' reduction from ',n+1,' to ',n-nn+1,' (factor ',n/(1.d0*(n-nn+1)),')'

END