;+
; Interactive procedure that allows user to use the plotter and mouse buttons
; to specify regions and gaussians for fitting.
; Left Mouse Button: to mark regions (left then right)
; Center Mouse Button: to mark gaussians (center then width)
; Right Mouse Button: to calculate and show fits
;
; @param fit {out}{optional}{type=array} results of fit. identical to !g.gauss.fit
; @param fitrms {out}{optional}{type=array} errors of fit. identical to !g.gauss.fitrms
; @keyword modelindex {in}{optional}{type=long}{default=-1} index of
; data container in !g to hold the resulting model.  If not set (the
; default) then no index is filled with the model.
;
; @examples
; <pre>
;   fitgauss,modelindex=10    ; fit gaussian and store model in DC 10
;   show
;   subtract,0,10,11          ; store residual in DC 11
;   oshow, 11                 ; overlay residual
;   oshow, 10, color=!yellow  ; overlay gaussian model
; </pre>
;-
pro fitgauss, fit, fitrms, modelindex=modelindex

    ; init output
    fit = -1
    fitrms = -1

    npts = !g.line ? data_valid(!g.s[0]) : data_valid(!g.c[0])
    if npts le 0 then begin
        message,'No valid data in primary data container',/info
        return
    endif
    data = !g.line ? *!g.s[0].data_ptr : *!g.c[0].data_ptr
    allChans = dindgen(npts)

    clearovers

    print, "**********************************************************"
    print, "Instructions for fitgauss procedure:"
    print, "Left Mouse Button: to mark regions to be fit."
    print, "Center Mouse Button: to mark initial guesses (center then width)"
    print, "Right Mouse Button: to calculate and show fits"
    print, "**********************************************************"
   
    start_state = 0
    region_state = 1
    ready_state = 2
    gauss_state = 3
    done_state = 4

    state = start_state

    new_region = dblarr(2)
    nregions = 0
    ngauss = 0
    
    new_gauss = dblarr(3)

    clearoplots
    
    while (state ne done_state) do begin
    
        ; wait for mouse press
        c = click(/nocrosshair)
        
        get_closest_index,c.chan,allChans,closest_index
        cx = closest_index 
        cy = c.y
        
;        if (nregions gt 0) then print, "Current regions: ", regions
        
        case state of
            ; just started
            start_state: begin
                case c.button of
                    ; left - mark out left region and transition
                    1: begin
                        new_region[0] = cx
                        state = region_state
                    end    
                    ; middle - not allowed
                    2: print, 'must select a region first'
                    ; right - let the user quite
                    4: state = done_state
                    else: ; do nothing
                endcase    
            end    
            ; must finish marking region 
            region_state: begin
                case c.button of
                    ; left - finish region and transition
                    1: begin
                        ; new point cant appear in an existing region
                        overlap_check,nregions,regions,cx,no_overlap
                        if no_overlap then begin
                            if (cx lt new_region[0]) then begin
                                ; region entered backwards, thats okay
                                new_region[1] = new_region[0]
                                new_region[0] = cx
                            endif else begin
                                new_region[1] = cx
                            endelse
                            ; append this new region to the list of regions
                            if (nregions eq 0) then begin
                                regions=[new_region]
                            endif else begin
                                regions=[regions,new_region]
                            endelse    
                            nregions += 1
                            print, "Limits accepted for region ",nregions
                            gbtoplot, seq(new_region[0],new_region[1]), data[new_region[0]:new_region[1]], color=!cyan, /chan
                            state = ready_state
                        endif else begin
                            print, 'Ending overlaps with other regions, try again.'
                        endelse
                    end
                    ; middle - not allowed
                    2: print, 'Must finish or quit region first.'
                    ; right - quit region
                    4: begin
                        print, 'Quiting region marking, right click again to quit procedure'
                        state = ready_state
                    end
                    else: ; do nothing
                endcase 
            end
            ; can do anything from here 
            ready_state: begin
                case c.button of
                    ; left - start a new region
                    1: begin
                        ;start of new region cant be within a pre-existing region
                        overlap_check,nregions,regions,cx,no_overlap
                        if no_overlap then begin
                            new_region[0] = cx
                            state = region_state
                        endif else begin
                            print, 'Cannot start a new region inside pre-existing one'
                        endelse    
                    end    
                    ; middle - start a new gauss
                    2: begin
                        ; is this inside a region, and which one?
                        find_region,nregions,regions,cx,r
                        if (r eq -1) then begin
                            print, "Warning: this Gaussian is not inside any of the regions."
                        endif
                        new_gauss[0] = cy
                        new_gauss[1] = cx
                        state = gauss_state
                    end    
                    ; right - exit
                    4: begin
                       state = done_state
                    end   
                    else: ; do nothing
                endcase 
            end
            ; must finish marking gauss 
            gauss_state: begin
                case c.button of
                    ; left - not allowed
                    1: print, 'Must finish marking gauss width or quit Gauss marking first.'
                    ; middle
                    2: begin
                        ; is this inside the right region?
                        find_region,nregions,regions,cx,r
                        if (r eq -1) then begin
                            print, "Warning: this Gaussian is not inside any of the regions."
                        endif  
                        ; need special code for computing width at this height
                        get_full_width_half_max,new_gauss[1],cy,data,fwhm
                        new_gauss[2] = fwhm 
                        ; add a new guass
                        if (ngauss eq 0) then gauss=[new_gauss] else gauss=[gauss,new_gauss]
                        ngauss += 1
                        print, 'Guesses accepted for gaussian ',ngauss
                        state = ready_state
                    end
                    ; right - quit gauss marking
                    4: begin
                        print, 'Quitting Gauss marking, right click again to quit procedure'
                        state = ready_state
                    end
                    else: ; do nothing
                endcase    
            end
            else: ; do nothing
        endcase
        
    endwhile
    
    ; set guide structure using the local variables we just set up
    !g.gauss.nregion = nregions
    !g.gauss.ngauss = ngauss
    if (nregions ne 0) then gregion, regions
    if (ngauss ne 0) then begin
        for i=0,ngauss-1 do begin
            gparamvalues,i, gauss[(i*3):(i*3)+2]
        endfor
    endif 

    if (ngauss gt 0) then begin
        ; fit these gaussians!
        gauss, fit, fitrms, modelindex=modelindex
        ; plot results
        gshow, modelindex=modelindex
    endif
    
end
    
;+
; @private
;-
pro find_region,nregions,regions,x,region

    region = -1
    if (nregions eq 0) then return

    for i=0,nregions-1 do begin
        b = regions[i*2]
        e = regions[(i*2)+1]
        if (x gt b) and (x lt e) then begin
            region = i
            return
        endif    
    endfor

    return
    
end

;+
; @private
;-
pro overlap_check,nregions,regions,checkpoint,status

    
    status = 1

    if (nregions eq 0) then begin
        return
    endif    

    for i=0,nregions-1 do begin
        b = regions[i*2]
        e = regions[(i*2)+1]
        if (checkpoint gt b) and (checkpoint lt e) then begin
            status = 0
        endif    
    endfor

    return

end   


;+
; @private
;-
pro get_closest_index, x, y, result

    ; what is this value between?
    a = where(x gt y)
    lhs_index = a[n_elements(a)-1]
    b = where(x lt y)
    rhs_index = b[0]

    ; which value is x closest to
    lhs_diff = x - y[lhs_index]
    rhs_diff = y[rhs_index] - x
    if lhs_diff ge rhs_diff then result = rhs_index else result = lhs_index
    
end

;+
; @private
;-
pro get_full_width_half_max, cx, cy, y, fwhm

    cy_gt_y = where(cy ge y)
    tail = cy_gt_y[where(cx lt cy_gt_y)]
    fwhm = (tail[0] - cx) * 2.0

end