;
;Copyright (C) 2005 Patricio Rojo
;
;This program is free software; you can redistribute it and/or
;modify it under the terms of the GNU General Public License
;as published by the Free Software Foundation; either version 2
;of the License, or (at your option) any later version.
;
;This program is distributed in the hope that it will be useful,
;but WITHOUT ANY WARRANTY; without even the implied warranty of
;MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;GNU General Public License for more details.
;
;You should have received a copy of the GNU General Public License
;along with this program; if not, write to the Free Software
;Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
;

;+
; NAME:
;	MAXTRACE
;
; PURPOSE:
;	This function fits the trace of some predominant features
;	along columns in a two dimensional array.
;
; CATEGORY:
;	Pato's Array Data Reduction Potpourri.
;
; CALLING SEQUENCE:
;
;	trace = MAXTRACE(Data, Y, Refy, Refx)
;
; INPUTS:
;       Data:   2 dimensional array in where to find the trace along
;               the columns
;       Y:      Array containing the value of the y axis
;       Refy:   Index of row where to start looking for the trace
;       Refx:   Index of column where to start looking for the trace
;
;
; KEYWORD PARAMETERS:
;       OSAMP:    Oversampling for peak finding. If using this option
;                 and a wavelet as a input you will want to use expy
;                 and be careful about converting back the output. If
;                 you want to use this switch make sure that the y
;                 scale is linear and equispaced.
;       EXPY:     Process through the natural logarithm of y
;                 coordinate.
;       FITNAME:  Name of the fit you want to use. Currently
;                 'true', 'gaussest', and 'gauss'
;       NOVAL:    Value of parameters when there was no succesfull fit
;       IGNORE:   Ignore data beyond this columns
;       TRACELIM: Force trace to be within these values, if it is
;                 going beyond this value, it goes back to refy value
;                 and start fitting in the opposite direction.
;
;
; OUTPUTS:
;	This function returns the parameters of the fitted trace
;
; OPTIONAL OUTPUTS:
;       TRUEBACK: Array containing the background when using trueshape
;                 fitting. 
;       PEAKS:    Array containing the Y-value of the highest peaks in
;                 the trace along the columns.
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;	Algorithm description.
;
; EXAMPLE:
;      Find the trace of array 'arr' using default gaussian fitting,
;   processing the background, given axis y is monospaced in
;   logarithmic space. Use position (cntx, xnty) as a reference to
;   start looking for the peak.
;
;       cntx = 220
;       cnty = 430
;	gaussprm = MAXTRACE(arr, y, cntx, cnty, /expy)
;
;
; MODIFICATION HISTORY:
; 	Written by:	Pato Rojo, Cornell.  Aug 2005
;			pato@astro.cornell.edu
;-

function peakborders, y, data, cy
compile_opt idl2, hidden

n = n_elements(y)
der = dblarr(n)
sder = dblarr(n)

der[2:n-3] = (data[3:n-2] - data[1:n-4]) / 2.0 + $
  (data[4:n-1] - data[0:n-5]) / 4.0

iyy = (where(y ge cy))[0]

;if it is at the right side of the peak
if der[iyy] lt 0 then begin
    lind = where(der[iyy:*] ge 0)
    if n_elements(lind) eq 1 then u = n - 1 $
    else begin
        pos = iyy + lind
;        pos = where(pos - shift(pos, 1) eq 2)
        while(pos[1] eq pos[0] + 2) do               $  
          if n_elements(pos) gt 2 then pos = pos[1:*] $
          else pos = [n-1, n-1]
        u = pos[0]
    endelse
    lind = where(der[0:iyy] ge 0)
    if n_elements(lind) eq 1 then return, [0, 1]  $ ;end of search
    else begin
        pos = rotate(lind, 2)
        while(pos[1] eq pos[0] - 2) do               $
          if n_elements(pos) gt 2 then pos = pos[1:*] $
          else return, [0, 1]
        lind = where(der[0:pos[0]-1] lt 0)
        if n_elements(lind) eq 1 then l = 0 $
        else begin
            pos = rotate(lind, 2)
            while(pos[1] eq pos[0] - 2) do               $
              if n_elements(pos) gt 2 then pos = pos[1:*] $
              else pos = [0, 0]
            l = pos[0]
        endelse
    endelse
;if it is at the left
endif else begin
    lind = where(der[0:iyy] lt 0)
    if n_elements(lind) eq 1 then l = 0 $
    else begin
        pos = rotate(lind, 2)
        while(pos[1] eq pos[0] - 2) do               $
          if n_elements(pos) gt 2 then pos = pos[1:*] $
          else pos = [0, 0]
        l = pos[0]
    endelse
    lind = where(der[iyy:*] lt 0)
    if n_elements(lind) eq 1 then return, [0, 1] $ ;end of search
    else begin
        pos = iyy + lind
        while(pos[1] eq pos[0] + 2) do               $
          if n_elements(pos) gt 2 then pos = pos[1:*] $
          else return, [0, 1]
        lind = where(der[pos[0]+1:*] ge 0)
        if n_elements(lind) eq 1 then u = n - 1 $
        else begin
            pos = pos[0] + 1 + lind
            while(pos[1] eq pos[0] + 2) do               $
              if n_elements(pos) gt 2 then pos = pos[1:*] $
              else pos = [n-1, n-1]
            u = pos[0]
        endelse
    endelse
endelse

return, [l, u]

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function findtrace, refy, inix, finx, stepx, cyy, tracelim, ignore, $
                    y, data, yfine, n, trueback, osamp, dogaussfit, $
                    traceprm

cy = refy
l  = -1

for x = inix, finx, stepx do begin
    nextiter = 0
    if(x gt ignore[1] or x lt ignore[0]) then break

    spdata = data[x, *]
    if(osamp gt 1) then begin
        sp     = spl_init(y, spdata)
        spdata = spl_interp(y, spdata, sp, yfine)
    endif

    ol = l
    l = peakborders(yfine, spdata, cy)

    if(l[0] ge l[1]) then break

    tmpd = max(spdata[l[0]:l[1]], cntidx)
    cntidx += l[0]
    cy   = yfine[cntidx]

;stop looking for trace if new center is outside previous borders
    if ol[0] ne -1 then $
      if(cntidx gt ol[1] or cntidx lt ol[0]) then break

;if I'm backsearching then finish doing so when locating previous
;point
    if n_elements(tracelim) eq 1 then begin
        if tracelim eq 1 then $
          if float(cyy[x]) eq float(cy) then return,x

;check whether I need to backtrace
    endif else if (n_elements(tracelim) eq 2) then begin
        if ((cy gt tracelim[1]) or (cy lt tracelim[0])) then begin
            b = findtrace(refy, x, inix, -stepx, cyy, 1, ignore, $
                          y, data, yfine, n,                     $
                          trueback, osamp, dogaussfit, traceprm)
        endif
    endif else $
      message, "TRACELIM must have either 1 or two elements"

    if dogaussfit gt 0 then begin
;prepare estimates
        wdt = yfine[l[1]] - yfine[l[0]]
        cnt = cy
        cidx = total(l) / 2     ;intdiv
        gausest = [spdata[cidx], cnt, wdt*0.2, 0]
        s = 0

                                ;process with the fitting
        switch dogaussfit of
            1: begin            ;Using fixed center
;Gaussian fit with fixed center is not working!! (curvefit() with fita
;doesn't work (see CVS 1.6))
                lgaussprm = gausest
                break
            end
            2: begin            ;Using fixed center v2
                vx = yfine[l[0]:l[1]]
                a = [gausest[0], gausest[2], gausest[3]]
                func = 'evengauss'

                gf = curvefit(vx-cnt, spdata[l[0]:l[1]],               $
                              replicate(1.0, l[1] - l[0] + 1), a,      $
                              status=s, functi=func, iter=iter)
                  
                if (s eq 1) or (a[0] le 0) or $
                  (a[1] le 0) or (a[1] ge 50 * gausest[2]) then begin
                    if ~ keyword_set(lbase) then begin
                        cy = refy
                        nextiter = 1
                        break
                    endif

                    lgaussprm = [(traceprm[x-stepx, 0] +        $
                                  (gausest[0] - lbase)) / 2.0,  $
                                 cnt,                           $
                                 (traceprm[x-stepx, 2] +        $
                                  gausest[2]) / 2.0,            $
                                 lbase]
                endif else $
                  lgaussprm = [a[0], cnt, a[1], a[2]]
                break
            end
            else: message, "Option not supported for gauss fitting" + $
              " 'dogaussfit'(" + string(dogaussfit) + ")"
        endswitch
        if (s ne 0) then break
        if nextiter then continue

        traceprm[x,*] = lgaussprm
        lbase = lgaussprm[3]
    endif else begin
                                ;Using real shape with borders where
                                ;the neighboring peak starts
        bslope = (spdata[l[1]] - spdata[l[0]]) / (yfine[l[1]] - yfine[l[0]])
        binter = spdata[l[1]] - bslope * yfine[l[1]] 
        c = ((long((l+1) / osamp + [0, 0.9999999])) < (n_elements(y)-1)) > 0
        traceprm[x, c[0]:c[1]] = data[x, c[0]:c[1]]
        trueback[x, *] = [binter, bslope]
    endelse

    cyy[x] = cy

endfor

return, x

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro evengauss, x, A, F, pder

nx = x / A[1]
ex = exp(-0.5 * nx * nx)
F = A[0] * ex + a[2]

if n_params() ge 4 then $
  pder = [[ex], [A[0] * ex * nx * nx / A[1]], $
          [replicate(1.0, n_elements(x))]]

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function maxtrace, data, y, refx, refy,                            $
                   osamp=osamp, fitname=fitname, expy=expy,        $
                   noval=noval, ignore=ignore, trueback=trueback,  $
                   tracelim=tracelim, peaks=peaks
;
;Finds trace of maximum point accross a continuous image
;
;Example:
;trace=maxtrace(data,y,refy,refx,/derivcenter,osamp=100,noval=0
;               ,/gaussianfit,/expy,ignore=[1,1000])
;

nx = (size(data,/dim))[0]
peaks = fltarr(nx)

goback = -1

if n_elements(noval) eq 0    then noval    = 0.0
if n_elements(tracelim) ne 2 then tracelim = 0
if ~ n_elements(ignore)      then ignore   = [0, nx]
if ~ keyword_Set(fitname)    then fitname  = 'gauss'


n = (size(data,/dim))[1]
if(n_elements(osamp)) then begin
    if keyword_set(expy) then yy = alog(y) $
    else yy = y
    span    = yy[n-1] - yy[0]
    osamp = long(osamp)
    yfine   = yy[0] + span *             $
      dindgen(long(osamp * (n-1) + 1)) / double(osamp * (n-1))
    if keyword_Set(expy) then yfine = exp(yfine)
    n       = n_elements(yfine)
endif else begin
    yfine   = y
    osamp = 1
endelse
lim = tracelim

case fitname of
    'gauss': begin
        dogaussfit = 2
        traceprm   = fltarr(nx, 4)
        nill       = replicate(noval, 4)
    end
    'gaussest': begin
        dogaussfit = 1
        traceprm   = fltarr(nx, 4)
        nill       = replicate(noval, 4)
    end
    'true' : begin
        dogaussfit = 0
        traceprm = dblarr(size(data,/dim))
        trueback = dblarr(nx, 2)
        nill     = replicate(noval, (size(data,/dim))[1])
    end
    else: message, "Unrecognized fit method '" + fitname + "'"
endcase

;Find the trace going forward and complete with nill what couldn't be
;fitted
x = findtrace(refy, refx, nx-1, 1, peaks, lim, ignore, y, data,  $
              yfine, n, trueback, osamp, dogaussfit, traceprm)
for x = x, nx-1, 1 do traceprm[x, *] = nill

;Find the trace going backward and complete with nill what couldn't be
;fitted
x = findtrace(refy, refx, 0, -1, peaks, lim, ignore, y, data,   $
              yfine, n, trueback, osamp, dogaussfit, traceprm)
for x = x, 0, -1 do traceprm[x, *] = nill

return, traceprm

end