;+ 
;NAME:
;
;  fancompress
;
;PURPOSE:
;  Decimates polylines in an aesthetically pleasing fashion.
;
;CALLING SEQUENCE:
;  outidx = fancompress(inpts,err)
;
;INPUT:
;  inpts: N x 2 dimension array, where inpts[*,0] are the x components of the polyline and inpts[*,1] are the y components of the polyline
;  err: The amount of error allowed before including a point 
; 
;
;OUTPUT:
;  An array of indexes into inpts.  Indices will range from 0 to N-1.  First and Last points are always included.
;
;NOTES:
;  1. Based almost entirely on the paper: 
;  Fowell, Richard A. and McNeil, David D. , “Faster Plots by Fan Data-Compression,” 
;  IEEE Computer Graphics & Applications, Vol. 9, No. 2,Mar. 1989, pp. 58-66.
;  
;  2. One modification from published algorithm, handles NaNs by always including the point
;  before a group of NaNs, 1 NaN and the point after the NaNs.  This ensures that gaps will
;  be drawn accurately.
;  
;  3. Algorithm is fairly slow, because it requires 1 pass over all data points.
;  Optimizing this algorithm by divide and conquer, vectorization, or dlm may be
;  a worthwhile use of time in the future. 
;
;$LastChangedBy: jimmpc $
;$LastChangedDate: 2009-05-29 15:09:16 -0700 (Fri, 29 May 2009) $
;$LastChangedRevision: 6003 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_5_1/misc/fancompress.pro $
;-----------------------------------------------------------------------------------


pro fn_save_pt,n_out,outpts,o,k,i,n_in

  compile_opt idl2,hidden

  if n_out eq 0 then begin
    outpts = i
  endif else begin
    outpts = [outpts,i]
  endelse 

  n_out++
  o = i
  k = n_in

end

function fn_dot_p,u,v

  compile_opt idl2,hidden
  
  return,total(u*v)
  
end

function fn_norm,v

  compile_opt idl2,hidden

  return,sqrt(fn_dot_p(v,v))

end

function fancompress,inpts,err

  compile_opt idl2
  
  n_in = (dimen(inpts))[0]
  n_out = 0
  keep = 1
  
  i = 1
  
  fn_save_pt,n_out,outpts,o,k,i,n_in
  
  if n_in eq 1 then begin
    return,inpts
  endif
  
  err_1 = max([0,err])

  while i lt n_in - 1 do begin
    i++
    ;properly mark NaNs in data, this addition is currently the only significant
    ;modification from the published algorithm
    if ~finite(inpts[i-1,0]) || ~finite(inpts[i-1,1]) then begin
      if keep eq 0 then begin
        outpts = [outpts,i-1]
        n_out++
      endif
      outpts = [outpts,i]
      n_out++
      for k = i,n_in-1 do begin
        if finite(inpts[k-1,0]) && finite(inpts[k-1,1]) then begin
          break
        endif
      endfor
      i = k
      if k eq n_in-1 then begin
        fn_save_pt,n_out,outpts,o,k,i,n_in
        break
      endif else if k eq n_in then begin
        break
      endif
      fn_save_pt,n_out,outpts,o,k,i,n_in
    endif
    
    if i lt k then begin
      distance = fn_norm(inpts[i-1,*] - inpts[o-1,*])
      if distance gt err_1 then begin
        k = i
        p_i_u = distance
        u_hat = (inpts[i-1,*] - inpts[o-1,*]) / p_i_u
        v_hat = [-u_hat[1],u_hat[0]]
        f = [p_i_u,err_1/p_i_u,-err_1/p_i_u]
      endif
    endif
    if i ge k then begin
      keep = 1
      p_i1_u = fn_dot_p((inpts[(i-1)+1,*] - inpts[(o-1),*]),u_hat)
      p_i1_v = fn_dot_p((inpts[(i-1)+1,*] - inpts[(o-1),*]),v_hat)
      if p_i1_u ge f[0] then begin
        m_i = p_i1_v/p_i1_u
        if m_i le f[1] && m_i ge f[2] then begin
          delta_m_i = err_1 / p_i1_u
          keep = 0
          f[0] = p_i1_u
          f[1] = min([f[1],m_i+delta_m_i])
          f[2] = max([f[2],m_i-delta_m_i])
        endif
      endif
      if keep then begin
        fn_save_pt,n_out,outpts,o,k,i,n_in
      endif
    endif
  endwhile
  
  outpts = [outpts,n_in]
  
  return,outpts-1

end