;+
; NAME:
;   viper_mesma
;
; AUTHOR:
;   Kerry Halligan
;   halligan@vipertools.org
;   www.vipertools.org
;
; LICENSE:
;   Copyright © 2005,2006,2007,2008 Kerry Halligan
;
;   This file is part of ViperTools.
;
;   ViperTools 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 3 of the License, or
;   (at your option) any later version.;
;
;   ViperTools 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 ViperTools.  If not, see <http://www.gnu.org/licenses/>.
;
; PURPOSE:
;   This routine performs Multiple Endmember Spectral Mixture Analysis (MESMA) on an input
;   image or images using endmembers containted in one or more ENVI spectral libraries. It
;   allows simple SMA or MESMA, photometric shade or non-photometric shade and a range of
;   constraints and outputs.
;
; CATEGORY:
;   MESMA
;   SMA
;
; CALLING SEQUENCE:
;
; INPUTS:
;   Reflectance image
;     Input image can be in any format (byte, integer, unsigned integer, floating point) or
;     interleave (bip, bil, bsq).  Integer or unsigned integer data should be in reflectance
;     times 10,000 (0 - 10,000). If data are in byte format, they should be in reflectance times
;     250 (0 - 250).  If data are in floating point they should be in reflectance (0 - 1).
;     Utilizes ENVI's band bands list (if present) to spectrally subset both image and spectral
;     library.
;   Spectral libraries
;     Up to 3 spectral libraries allowed.  Need to be the same number of bands and same data
;     type (e.g. floating point, integer, etc.) as image.  Libraries should not contain shade.
;     Most common libraries would be 1) all spectra (used for 2em mode), 2) vegetation (used
;     for 3em and 4em mode, and 3) npv + soil (used for 3em case) 4) npv (used for 4em case) and
;     5) soil (used for 4 em case).  Note that all combinations of spectra are used, so a
;     4em run with 10 green veg spectra, 5 npv and 6 soil spectra runs a total of 300 models.
;
; OPTIONAL INPUTS:
;   None
;
; KEYWORD PARAMETERS:
;   None
;
; OUTPUTS:
;   Minimum RMS image: Non-shade and shade fractions (1st bands) plus the rms and model number
;     of the minimum rms model.
;   Classification image: Classified image with the minimum RMS model for each pixel.
;
; OPTIONAL OUTPUTS:
;   None
;
; PROCEDURE:
;   The proceedure builds a lookup table for the endmembers for each model, reads in all spectral
;   libraries, then begins a line by line loop.  For each image line it reads in the data, and then
;   loops through each model.  For each model the spectra are selected and used to build an
;   endmember array, which is passed to viper_mesma_fracCalc.  fraccalc returns the fractional abundances
;   of all non-shade and shade endmembers, the model RMSE and optionally the residuals.  If
;   the current model produced a lower RMSE for any given pixel, the selected fraction, RMSE,
;   and residual constraints are tested.  All pixels that meet the constraints are considered
;   valid.  All pixels that produce a lower RMSE than the stored best value AND are valid are
;   updated with the new model number, fractions and RMSE.  After all models have been run,
;   results for that line of image data are saved to disk and the next line is read in.
;   This procedure selects the single model for each pixel that meets all constraints AND
;   has the lowest RMS error.  If no model meets the contraints then pixel is left with all
;   zero values and appears as 'unclassified' in the output image.  When the image is complete,
;   the file is closed, headers are written, and a classification image is produced. Calculates
;   fractions first, using Singular Value Decomposition to invert the endmember matrix, then
;   calculates shade as 1-sum of non-shade fractions.  If a non-photometric shade endmember is
;   used, then this endmember is subtracted first from each endmember, then from the image spectrum
;   then the fractions are calculated as above.
;
; REQUIRES:
;   sel_from_list
;   cmapply
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;   6-30-05 Written by Kerry Halligan from eMESMA.pro, with earlier versions of this code
;       dating back to 2001 and used for Master's thesis work
;   9-29-05 Kerry Halligan modified including
;       significant debugging effort to address a range of problems and to add
;       batch mode capabilities.  Also removed some of the cmapply calls which
;       were causing an unknown error during the constraint process
;   6-30-06 to fix problems with dropdown list for non-photometric shade selection
;   9-6-06 Changed counters from integer to float datatype to allow for greater
;       than 2^16/2 models (to accomidate Ted Eckman's fire temp mapping work)
;   11-1-06  Made various changes to fix two problems with the residual images
;       removed r_fid capture if/when residual image is opened
;       added a case statement for scaling residual images such that byte, integer and floating
;       point are handled explicitly, all others are treated as floating point
;   11-21-06 Kerry Halligan modified to do the following:
;       added error handler to _run routine to report error message
;       activated the save and restor control files after major re-writing of these routines
;   12-6-06 Kerry Halligan modified to first check all input files to make sure they
;     exist before trying to load them in ENVI - this prevents crashes when files not found.
;     Also added format statements to printf calls for filenames of image and spectral
;     libraries that were causing new lines when long filenames/paths were used.
;   12-20-06 Kerry Halligan odified add back in the thres_resids function that had been
;     inadvertantly removed from previous version when cleaning up code.  Also fixed bug
;     in the 'run' routine which was preventing the read of the text widgets for the
;     max RMSE and max residual text widgets.  Now it is no longer necessary to hit
;     enter to update these values.  Changed output residual image to now be:
;      same number of bands as input image, regardless of bad bands list
;      bands that are denoted as 'bad' in bad bands list will be zeros in residual image
;      floating point data with no scale factor (e.g. DN,radiance, reflectance)
;          regardless of scale factor of input image
;      renamed fracCalc2 to viper_mesma_fraccalc
;      renamed thresh_resids to viper_mesma_thresh_resids
;   12-21-06 Kerry Halligan modified to fix output residuals.  Bug had been resulting in the
;       output residuals being from just the last model, not the best model.
;   1-1-07 Kerry Halligan made the following changes
;     renamed to viper_mesma from eMESMA_viper_batch to maintain consistancy
;     fixed a bug in residual calculation that was causing it to crash when output residuals not selected
;     made significant modification to the handling of input datasets with regards to bad bands lists - copied code
;      from recently updated CRES module.  Now uses bad bands common to all spectral libraries and input image
;     warns user if bad bands are different between spectral libraries and image, and asks user if new bad bands
;      list should be output as an ASCII file.  This file can be imported in ENVI
;   1-7-07 Kerry Halligan fixed bug in the restore control file routine (viper_mesma_restore) which was causing
;     the index for the scale factor (0-3) to be stored rather than the scale factor itself (0,1,1000,10000)
;   2-27-07 Kerry Halligan fixed bug in the residual constraint.  Dropped code which autmatically set number
;     of bands to the number of good bands, and changed so that it only checks to see if number of good bands
;     is less than number of residual bands IF residual constraint is used.
;   7-3-07 Kerry Halligan made miner change to default file names in the 'browse' and 'save' routines to speed
;     selection of input and output filenames.
;-

; The file browse routine for file selection
pro viper_mesma_browse, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; Retrieve the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
inpath=(*pstate).inpath
case event.id of
    (*pstate).file0browse: begin
            ;* Get the file
            infile=dialog_pickfile(/must_exist, path=inpath, get_path=gp, /read, /multiple)
            ; open up the file in envi (if multiple, just open first)
            envi_open_file, infile[0], r_fid=fid, /no_realize
            ; query intensity image for dimensions
            envi_file_query, fid, ns=ns, file_type=filetype, nl=nl, nb=nb, sname=sname
            if fid ne -1 then begin
                if filetype le 1 then begin
                    (*pstate).fid0=fid
                    ; set the text widget to display the selected file
                    widget_control, (*pstate).file0text, set_value=sname
                    (*pstate).infile=ptr_new(infile)
                    (*pstate).inpath=gp
                    slidemax= 20 < nb
                    widget_control, (*pstate).residb_slider, SET_SLIDER_MAX=slidemax
                    ;(*pstate).constraints[5]=slidemax
                endif else begin
                    result=dialog_message('Not an ENVI standard file, please enter new file')
                    widget_control, (*pstate).file0text, set_value='Required'
                    (*pstate).fid0=-1
                endelse
            endif else begin
                result=dialog_message('No valid ENVI file selected, please select file')
                widget_control, (*pstate).file0text, set_value='Required'
                (*pstate).fid0=-1
            endelse

           ; if multiple images selected, set outfile automatically and grey out
           if n_elements(infile) gt 1 then begin
              ; strip off any extensions
              outfile=strarr(2,n_elements(infile))
              for i=0,n_elements(infile)-1 do outfile[*,i]=strsplit(infile[i], '.', /extract)
              outfile=outfile[0,*]+'_sma'
              ; add 'sma' to end of filenames
              (*pstate).outfile=ptr_new(outfile)
              ; grey out outfile widgets
              widget_control, (*pstate).file4browse, sensitive=0
              widget_control, (*pstate).sel4lab, sensitive=0
              widget_control, (*pstate).file4text, set_value=outfile, editable=0
              (*pstate).multiple=1
           endif else begin
              widget_control, (*pstate).file4browse, sensitive=1
              widget_control, (*pstate).sel4lab, sensitive=1
              widget_control, (*pstate).file4text, set_value='', editable=1
              (*pstate).outfile=ptr_new()
              (*pstate).multiple=0
           endelse

       end

    (*pstate).file1browse: begin
        ; Get the file
        infile=dialog_pickfile(/must_exist, path=inpath, /read, get_path=gp)
        ; open up the file in envi
        envi_open_file, infile, r_fid=fid, /no_realize
        ; query intensity image for dimensions
        envi_file_query, fid, ns=nspec, file_type=filetype, spec_names=specnames, sname=sname, fname=fname
        if fid ne -1 then begin
            if filetype eq 4 then begin
                (*pstate).inpath=gp
                (*pstate).fid1=fid
                ; set the text widget to display the selected file
                widget_control, (*pstate).file1text, set_value=sname
                (*pstate).file1fulltext=fname
                ; get and set specnames
                (*pstate).specnames1=ptr_new(specnames)
                (*pstate).selected1=ptr_new(bytarr(n_elements(specnames))+1)
            endif else begin
                result=dialog_message('Not an ENVI spectral library, please enter new file')
                widget_control, (*pstate).file1text, set_value='Required'
                (*pstate).file1fulltext=''
                (*pstate).fid1=-1
            endelse
        endif else begin
           result=dialog_message('No valid ENVI file selected, please select file')
           widget_control, (*pstate).file1text, set_value='None'
           (*pstate).file1fulltext=''
           (*pstate).fid1=-1
           (*pstate).specnames1=ptr_new()
           (*pstate).selected1=ptr_new()
        endelse
    end
    (*pstate).file2browse: begin
        ; Get the file
        infile=dialog_pickfile(/must_exist, path=inpath, /read)
        ; open up the file in envi
        envi_open_file, infile, r_fid=fid, /no_realize
        ; query intensity image for dimensions
        envi_file_query, fid, ns=nspec, file_type=filetype, spec_names=specnames, sname=sname, fname=fname
        if fid ne -1 then begin
            if filetype eq 4 then begin
                (*pstate).fid2=fid
                ; set the text widget to display the selected file
                widget_control, (*pstate).file2text, set_value=sname
                (*pstate).file2fulltext=fname
                ; get and set specnames
                (*pstate).specnames2=ptr_new(specnames)
                (*pstate).selected2=ptr_new(bytarr(n_elements(specnames))+1)
            endif else begin
                result=dialog_message('Not an ENVI spectral library, please enter new file')
                widget_control, (*pstate).file2text, set_value='None'
                (*pstate).file2fulltext=''
                (*pstate).fid2=-1
            endelse
        endif else begin
            result=dialog_message('No valid ENVI file selected, please select file')
            widget_control, (*pstate).file2text, set_value='None'
         (*pstate).file2fulltext=''
            (*pstate).fid2=-1
            (*pstate).specnames2=ptr_new()
            (*pstate).selected2=ptr_new()
        endelse
    end
    (*pstate).file3browse: begin
        ; Get the file
        infile=dialog_pickfile(/must_exist, path=inpath, /read)
        ; open up the file in envi
        envi_open_file, infile, r_fid=fid, /no_realize
        ; query intensity image for dimensions
        envi_file_query, fid, ns=nspec, file_type=filetype, spec_names=specnames, sname=sname, fname=fname
        if fid ne -1 then begin
            if filetype eq 4 then begin
                (*pstate).fid3=fid
                ; set the text widget to display the selected file
                widget_control, (*pstate).file3text, set_value=sname
                (*pstate).file3fulltext=fname
                ; get and set specnames
                (*pstate).specnames3=ptr_new(specnames)
                (*pstate).selected3=ptr_new(bytarr(n_elements(specnames))+1)
            endif else begin
                result=dialog_message('Not an ENVI spectral library, please enter new file')
                widget_control, (*pstate).file3text, set_value='None'
                (*pstate).file1fulltext=''
                (*pstate).fid3=-1
            endelse
        endif else begin
            result=dialog_message('No valid ENVI file selected, please select file')
            widget_control, (*pstate).file3text, set_value='None'
            (*pstate).file1fulltext=''
            (*pstate).fid3=-1
            (*pstate).specnames3=ptr_new()
            (*pstate).selected3=ptr_new()
        endelse
    end
    (*pstate).file5browse: begin
        ; Get the file
        infile=dialog_pickfile(/must_exist, path=inpath, /read)
        ; open up the file in envi
        envi_open_file, infile, r_fid=fid, /no_realize
        ; query intensity image for dimensions
        envi_file_query, fid, ns=nspec, file_type=filetype, spec_names=specnames, sname=sname, fname=fname
        if fid ne -1 then begin
            if filetype eq 4 then begin
                (*pstate).fid5=fid
                ; set the text widget to display the selected file
                widget_control, (*pstate).file5text, set_value=sname
                (*pstate).file5fulltext=fname
                ; get and set specnames
                (*pstate).specnames5=ptr_new(specnames)
                ;(*pstate).selected5=ptr_new(bytarr(n_elements(specnames))+1)
                widget_control, (*pstate).file5subset, set_value=['Select spectrum', specnames]
            endif else begin
                result=dialog_message('Not an ENVI spectral library, please enter new file')
                widget_control, (*pstate).file5text, set_value='None'
                (*pstate).file5fulltext=''
                (*pstate).fid5=-1
            endelse
        endif else begin
            result=dialog_message('No valid ENVI file selected, please select file')
            widget_control, (*pstate).file5text, set_value='None'
            (*pstate).file5fulltext=''
            (*pstate).fid5=-1
            (*pstate).specnames1=ptr_new()
            (*pstate).selected1=ptr_new()
        endelse
    end
    (*pstate).file4browse: begin
        ; Get the filename for the output fraction image
        outfile=dialog_pickfile(/overwrite_prompt, path=inpath, /write)
        (*pstate).outfile=ptr_new(outfile)
        ; Display filename in widget
        widget_control, (*pstate).file4text, set_value=outfile
    end
endcase
filelist=[(*pstate).fid1, (*pstate).fid2, (*pstate).fid3]
sellist=lonarr(3)
if ptr_valid((*pstate).selected1) then sellist[0]=total(*(*pstate).selected1)
if ptr_valid((*pstate).selected2) then sellist[1]=total(*(*pstate).selected2)
if ptr_valid((*pstate).selected3) then sellist[2]=total(*(*pstate).selected3)
;sellist=[total(*(*pstate).selected1), total(*(*pstate).selected2), total(*(*pstate).selected3)]
index=where(filelist gt -1 and sellist gt 0)
if index[0] ne -1 then begin
    sellist=sellist[index]
    nmods=cmapply('*',sellist)
endif else begin
    nmods=0
endelse
widget_control, (*pstate).total_txt, set_value=string(long(nmods))

end


; Subset routine that allows the user to subset the spectral libraries and
; calculates the total number of models that will be run in the.  Uses the
; custom sel_from_list GUI widget function found in sel_from_list.pro
pro viper_mesma_subset, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; Retrieve the pointer to the state structure
widget_control, event.top, get_uvalue = pstate

case event.id of
    (*pstate).file1subset: begin
       ; get specnames ptr from pstate
       specnames=(*pstate).specnames1
       ; make sure it is valid pointer
       valid=ptr_valid(specnames)
       if valid eq 0 then begin
         result=dialog_message('Must select spectral library first')
         return
       endif
       specnames=*specnames
       selected1=(*pstate).selected1
       valid=ptr_valid(selected1)
       if valid eq 0 then begin
         result=dialog_message('Must select spectral library first')
         return
       endif
       selected1=*selected1
        if total(selected1) gt 0 then begin
               presel=where(selected1 eq 1)
               sel=sel_from_list(specnames, title='Subset library 1', presel=presel)
            endif else begin
               sel=sel_from_list(specnames, title='Subset library 1')
            endelse
       (*pstate).selected1=ptr_new(sel)
       ; update if window wasn't canceled
       ;if result[0] ne -1 then (*pstate).selected1=result eq 1
       end
        (*pstate).file2subset: begin
            ; get specnames ptr from pstate
             specnames=(*pstate).specnames2
            ; make sure it is valid pointer
            valid=ptr_valid(specnames)
            if valid eq 0 then begin
              result=dialog_message('Must select spectral library first1')
              return
            endif
            specnames=*specnames
            selected2=(*pstate).selected2
            valid=ptr_valid(selected2)
            if valid eq 0 then begin
              result=dialog_message('Must select spectral library first2')
              return
            endif
            selected2=*selected2
            if total(selected2) gt 0 then begin
               presel=where(selected2 eq 1)
               sel=sel_from_list(specnames, title='Subset library 2', presel=presel)
            endif else begin
               sel=sel_from_list(specnames, title='Subset library 2')
            endelse

            (*pstate).selected2=ptr_new(sel)
            ; update if window wasn't canceled
            ;if result[0] ne -1 then (*pstate).selected1=result eq 1
        end
        (*pstate).file3subset: begin
            ; get specnames ptr from pstate
            specnames=(*pstate).specnames3
            ; make sure it is valid pointer
            valid=ptr_valid(specnames)
            if valid eq 0 then begin
              result=dialog_message('Must select spectral library first1')
              return
            endif
            specnames=*specnames
            selected3=(*pstate).selected3
            valid=ptr_valid(selected3)
            if valid eq 0 then begin
              result=dialog_message('Must select spectral library first2')
              return
            endif
            selected3=*selected3
            if total(selected3) gt 0 then begin
                presel=where(selected3 eq 1)
               sel=sel_from_list(specnames, title='Subset library 3', presel=presel)
            endif else begin
               sel=sel_from_list(specnames, title='Subset library 3')
            endelse
       (*pstate).selected3=ptr_new(sel)
      ; print,*(*pstate).selected3
       ; update if window wasn't canceled
       ;if result[0] ne -1 then (*pstate).selected1=result eq 1
       end
      (*pstate).file5subset: begin
        ; query non-photometric shade spectral library
;        envi_file_query, (*pstate).fid5, spec_names=specnames
;        if total((*pstate).selected5) gt 0 then begin
;                    presel=where(selected1 eq 1)
;               sel=sel_from_list(specnames, title='Select specra', presel=presel)
;            endif else begin
;               sel=sel_from_list(specnames, title='Select specra')
;            endelse
       sel=widget_info((*pstate).file5subset,/droplist_select)
       (*pstate).selected5=sel
           ; update if window wasn't canceled
;        if result[0] ne -1 then begin
;            if total(result) eq 1 then begin
;                (*pstate).selected5=result eq 1
;            endif else begin
;                 result=dialog_message('Select only 1 non-photometric shade spectrum')
;                 return
;            endelse
;        endif
       end
endcase

filelist=[(*pstate).fid1, (*pstate).fid2, (*pstate).fid3]
sellist=lonarr(3)
if ptr_valid((*pstate).selected1) then sellist[0]=total(*(*pstate).selected1)
if ptr_valid((*pstate).selected2) then sellist[1]=total(*(*pstate).selected2)
if ptr_valid((*pstate).selected3) then sellist[2]=total(*(*pstate).selected3)
;sellist=[total(*(*pstate).selected1), total(*(*pstate).selected2), total(*(*pstate).selected3)]
index=where(filelist gt -1 and sellist gt 0)
if index[0] ne -1 then begin
    sellist=sellist[index]
    nmods=cmapply('*',sellist)
endif else begin
    nmods=0L
endelse
widget_control, (*pstate).total_txt, set_value=string(long(nmods))
end


;This proceedure sets the reflectance scale factor
pro viper_mesma_scale, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
list=[0,1,1000,10000]
if event.index eq 0 then begin
    result=dialog_message('Select a reflectance scale factor')
    return
endif
; set scale factor
(*pstate).scale=list[event.index]
print,'scale factor: ',(*pstate).scale
end


; This routine handles the toggle event selecting between photometric shade and non-
; photometric shade
pro viper_mesma_shade, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
if event.select eq 1 then begin
    case event.id of
        (*pstate).photo_shade_but:begin
        ; set shade to 1 and make widgets unsensitive
        (*pstate).photo_shade=1
        widget_control, (*pstate).file5text, sensitive=0
        widget_control, (*pstate).sel5lab, sensitive=0
        widget_control, (*pstate).file5browse, sensitive=0
        widget_control, (*pstate).file5subset, sensitive=0
        (*pstate).fid5=-1
        end
        (*pstate).nonphoto_shade_but:begin
        (*pstate).photo_shade=0
        widget_control, (*pstate).file5text, sensitive=1
        widget_control, (*pstate).sel5lab, sensitive=1
        widget_control, (*pstate).file5browse, sensitive=1
        widget_control, (*pstate).file5subset, sensitive=1
        end
    endcase
endif
print,(*pstate).photo_shade
end


; This proceedure gets the change of status of slider checkboxes, updates the
; 'checks' vector and sets sensitivity of sliders
pro viper_mesma_check, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
; get widget id of slider
case event.id of
    (*pstate).check1:begin
        (*pstate).checks[0]=event.select
        widget_control, (*pstate).p1_slide, sensitive=event.select
        end
    (*pstate).check2:begin
        (*pstate).checks[1]=event.select
        widget_control, (*pstate).p2_slide, sensitive=event.select
        end
    (*pstate).check3:begin
        (*pstate).checks[2]=event.select
        widget_control, (*pstate).p3_slide, sensitive=event.select
        end
    (*pstate).check4:begin
        (*pstate).checks[3]=event.select
        widget_control, (*pstate).maxrmse_text, sensitive=event.select
        widget_control, (*pstate).maxrmse_lab, sensitive=event.select
        end
    (*pstate).check5:begin
        (*pstate).checks[4]=event.select
        widget_control, (*pstate).residv_text, sensitive=event.select
        widget_control, (*pstate).residv_lab, sensitive=event.select
        widget_control, (*pstate).residb_slider, sensitive=event.select
        end
endcase
end


; This proceedure gets the change of status of output checkboxes, updates the
; 'output' vector
pro viper_mesma_out, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
; get widget id of slider
case event.id of
    (*pstate).out_frac:(*pstate).outputs[0]=event.select
    (*pstate).out_class:(*pstate).outputs[1]=event.select
    (*pstate).out_resid:(*pstate).outputs[2]=event.select
endcase
end


; This proceedure that gets and sets all sliders for MESMA constraints
pro viper_mesma_slider, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
; widget id of slider that was moved is in event.id
; compare to slider id's to determine which was moved
; and then update approrpriate info
widget_control, event.id, GET_VALUE = value
case event.id of
    (*pstate).min_slider:(*pstate).constraints[0]=float(round(value*100))/100
    (*pstate).max_slider:(*pstate).constraints[1]=float(round(value*100))/100
    (*pstate).shade_slider:(*pstate).constraints[2]=float(round(value*100))/100
;    (*pstate).maxrmse_slider:begin
;        value=round(value*10)/10.0
;        (*pstate).constraints[3]=float(value)/100.0
;        value=round(value*10)/10.0
;        widget_control, event.id, SET_VALUE = value
;        end
;    (*pstate).residv_slider:begin
;        value=round(value*10)/10.0
;        (*pstate).constraints[4]=float(value)/100.0
;        widget_control, event.id, SET_VALUE = value
;        end
    (*pstate).residb_slider:(*pstate).constraints[5]=value
endcase
print,'Constraints: ',(*pstate).constraints
end


; This proceedure that gets and sets the RMSE and Residual MESMA constraints
pro viper_mesma_text, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate
; widget id of text box that was updated in event.id

widget_control, event.id, GET_VALUE = value
case event.id of
    (*pstate).maxrmse_text:begin
        (*pstate).constraints[3]=float(value)
         widget_control, event.id, SET_VALUE = string(float(value))
        end
    (*pstate).residv_text:begin
         (*pstate).constraints[4]=float(value)
         widget_control, event.id, SET_VALUE = string(float(value))
        end
endcase
print,'Constraints: ',(*pstate).constraints
end


; This routine gets and sets the 'open files' toggle that causes the routine to
; or not open files for viewing or not
pro viper_mesma_open, event
COMPILE_OPT idl2, HIDDEN
on_error, 2

; Retrieve the pointer to the state structure
widget_control, event.top, get_uvalue = pstate

; flip sign on open_files variable
(*pstate).open_files = -(*pstate).open_files

end


;+
; This function is called by the MESMA code and processes on a single image line
; It gets passed the endmember array that is n x m where n is the number of
; endmembers and m is the number of bands per spectrum.  The resids keyword can
; be set to a named variable that will contain the riduals for all bands for all
; spectra in the input image_line.  viper_mesma_fraccalc makes several improvements over the
; previously used fracCalc including much reduced looping.  Returns a structure
; with the following tags: frac, shade, rmse.  Frac is a n X m array where n is the
; number of non-shade endmembers, and m is the number of pixels in the input image
; line, and the values are the fractional abundance of each endmember for each pixel,
; in reflectance units.  Shade and RMSE are m element vectors where m is the number
; of pixels in the input image line, and the values are the fractional abundance of
; shade and the RMSE in reflectance units.  See main viper_mesma routine for an
; example of using viper_mesma_fraccalc.
;-
function viper_mesma_fraccalc, em_array, image_line, resids=resids
compile_opt idl2
on_error, 2

; get the data type of the input endmember array
data_type = size(em_array, /type)

; get the number of endmembers and number of bands
dims = size(em_array, /dimensions)
n_ems = dims[0]
nb = dims[1]

; run a singular value decomposition on the endmember array where:
; w = n element vector containing the eigenvalues
; u = n column by m row orthogonal array used in decomposition
; v = n column by n row orthogonal array used in decomposition
svdc,em_array,w,u,v

; fill w_inv as diagonal matrix where values are 1 over the eigenvalues
w_inv=(1/w)##(bytarr(n_ems)+1)*identity(n_ems)

; calculate inverse matrix
em_inv = v##w_inv##transpose(u)

; do matrix multiplication of em_mat_inv matrix by the image spectra
frac = em_inv##image_line

; calculate the shade fraction, while constraining the fractions to sum to one
; by subtracting the sum of the fractional abundances of all other EMs from one.
if n_ems eq 1 then shade=1-frac else shade=1-total(frac,2)
model = em_array##frac ; calculate the modeled spectra
resids = image_line - model

; calculate root mean square error from residuals
rmse = sqrt(total(resids^2,2)/nb)

; place fractions (ems and shade) and rmse in output structure
sma={frac:frac, $
    shade:shade, $
    rmse:rmse}

return, sma
end


;+
; This routine screens a residuals vector for concecutive bands that exceed a given
; threshold.  Inputs are a residual vector, a threshold value and the number of
; consecutive bands that should be screened for.  Output is either a 0, if value
; was exceeded over sufficient consecutive bands and 1 if not.  This routine can be called by
; cmapply to run in 2nd dimension to allow for all pixels in an image line to be
; processed without an IDL loop.  See viper_mesma for an example
;-
function viper_mesma_thresh_resids, resids, value=value, nbands=nbands
compile_opt idl2
on_error, 2

if total(erode(abs(resids) gt value,bytarr(nbands)+1)) ge 1 then result=0 else result=1
return, result
end


; The _run routine for the viper_mesma
; See main header above for details
pro viper_mesma_run, event
compile_opt idl2
forward_function widget_auto_base, widget_menu, auto_wid_mng
;on_error, 2

; Establish error handler
catch, theError
if theError ne 0 then begin
    catch, /cancel
    help, /last_message, output=errText
    errMsg = dialog_message(errText, /error, title='Error processing request')
    return
endif

; get the pointer to the state structure
widget_control, event.top, get_uvalue = pstate

; get the output image filename and check that at least one output file exists
outfile=(*pstate).outfile
; dereference pointer
if ptr_valid(outfile) then begin
    outfile = *outfile
    ; if in single file mode, then get the output filename from the widget rather than pstate
    if (*pstate).multiple eq 0 then begin
       widget_control, (*pstate).file4text, get_value=outfile
       if strlen(outfile[0]) eq 0 then begin
            result=dialog_message('No output file selected')
            return
        endif
        ; make sure output path is valid
        tmp=file_dirname(outfile[0])
        result=file_test(tmp)
        if result eq 0 then begin
            result=dialog_message('Output directory does not exist, create it?',/cancel)
            if result eq 'OK' then begin
                file_mkdir,tmp
           endif else begin
             return
           endelse
        endif
    endif
    if strlen(outfile[0]) eq 0 then begin
        result=dialog_message('No output file selected')
        return
    endif
    ; find files already open in ENVI and make sure outfile is not in list
    fids = envi_get_file_ids()
    nfids = n_elements(fids)
    if (fids[0] ne -1) then begin
        for i = 0L, nfids - 1 do begin
            envi_file_query, fids[i], fname = fname
            if n_elements(fnamelist) eq 0 then fnamelist=fname else fnamelist=[fnamelist,fname]
        endfor
        tmp=strarr(nfids)
        for i=0L,nfids-1 do tmp[i]= where(outfile eq fnamelist[i])
        if total(tmp ge 0) gt 0 then begin
            result=dialog_message('Output file(s) already open in ENVI.  Select new output filename or close files in ENVI and then run, file(s) will be overwritten.')
            return
        endif
    endif
endif else begin
    result=dialog_message('No output file selected')
    return
endelse

; get the selected (subset) spectra for each of the 3 speclibs
if ptr_valid((*pstate).selected1) then selected1=where(*(*pstate).selected1 eq 1) else selected1=-1
if selected1[0] ne -1 then nsel1=n_elements(selected1) else nsel1=0
if ptr_valid((*pstate).selected2) then selected2=where(*(*pstate).selected2 eq 1) else selected2=-1
if selected2[0] ne -1 then nsel2=n_elements(selected2) else nsel2=0
if ptr_valid((*pstate).selected3) then selected3=where(*(*pstate).selected3 eq 1) else selected3=-1
if selected3[0] ne -1 then nsel3=n_elements(selected3) else nsel3=0
if (*pstate).selected5 eq -1 then nsel4=0 else begin
    nsel4=1
    selected4=(*pstate).selected5-1 ; subtract 1 for the 'Select spectrum' string at start of list
endelse

if nsel1 eq 0 then begin
    result=dialog_message('Library 1 must have a minimum of 1 spectrum selected')
    return
endif
if ((nsel3 gt 0) and (nsel2 eq 0)) then begin
    result=dialog_message('Library 2 must have a minimum of 1 spectrum selected')
    return
endif
;if nsel5 eq 0 then result=dialog_message('Shade library must have a spectrum selected')
;if nsel5 gt 1 then begin
;    result=dialog_message('Shade library must have only one spectrum selected')
;    return
;endif

; get total number of models to run
widget_control, (*pstate).total_txt, get_value=nmods
nmods=long(nmods[0])

; get contstraints from pstate
checks=(*pstate).checks
; update text input constraints in case user did not hit return after entering value(s)
widget_control, (*pstate).maxrmse_text, GET_VALUE = value
widget_control, (*pstate).maxrmse_text, SET_VALUE = string(float(value))
(*pstate).constraints[3]=float(value)
widget_control, (*pstate).residv_text, GET_VALUE = value
widget_control, (*pstate).residv_text, SET_VALUE = string(float(value))
(*pstate).constraints[4]=float(value)
; transfer constraints to new variables
min_frac = (*pstate).constraints[0]
max_frac = (*pstate).constraints[1]
max_shade = (*pstate).constraints[2]
max_rms = (*pstate).constraints[3]
max_resid = (*pstate).constraints[4]
resid_count = (*pstate).constraints[5]

; find out if files should be opened in ENVI
open_files=(*pstate).open_files

; get file IDs and scale factor from pstate
scale = (*pstate).scale
if scale eq 0 then begin
    result=dialog_message('Select a reflectance scale factor')
    return
endif
sl_fid1 = (*pstate).fid1
sl_fid2 = (*pstate).fid2
sl_fid3 = (*pstate).fid3

; get the selected outputs
outputs = (*pstate).outputs

;
if outputs[1] eq 1 and nmods ge 2.0^16/2.0 then begin
    result=dialog_message('Cannot produce output classification when number of models exceeds 32,767. Uncheck ''Classification Image'' output and re-run.')
    return
endif

; get number of input spectral libraries (number of non-shade endmembers in each model)
n_ems = total([sl_fid1, sl_fid2, sl_fid3] ne -1)
out_nb = n_ems+3

; get the list of input files
infile=(*pstate).infile
; dereference pointer
infile=*infile
; count files
nfiles=n_elements(infile)

; make sure that an input image and at least one speclib have been selected
envi_open_file, infile[0], r_fid=image_fid, /no_realize
if image_fid eq -1 then begin
    result=dialog_message('No image selected')
    return
endif
if sl_fid1 eq -1 then begin
    result=dialog_message('No spectral libraries selected')
    return
endif


;Build progress bar
ss=get_screen_size()
base=widget_base(/col, xoffset=ss[0]*.5-150, yoffset=ss[1]*.5, title='SMA/MESMA Progress:', /align_center)
label1=widget_label(base,value='Gathering information     ')
label2=widget_label(base,value='Percent completed:')
slider=widget_slider(base, xsize=300, minimum=0, maximum=100)
widget_control, base, /realize

; query 1st spectral library
envi_file_query, sl_fid1, spec_names=spec_names1, nl=n_spec1, sname=speclib_sname1, ns=sl_bands1, bbl=sl_bbl1
speclib1=envi_get_data(fid=sl_fid1, pos=0, dims=[-1,0,sl_bands1-1,0,n_spec1-1])
; if no bad bands list then make it
if sl_bbl1[0] eq -1 then sl_bbl1=bytarr(sl_bands1)+1

; query 2nd spectral library
if sl_fid2 ne -1 then begin
    envi_file_query, sl_fid2, spec_names=spec_names2, nl=n_spec2, sname=speclib_sname2, ns=sl_bands2, bbl=sl_bbl2
    speclib2=envi_get_data(fid=sl_fid2, pos=0, dims=[-1,0,sl_bands2-1,0,n_spec2-1])
    ; if no bad bands list then make it
    if sl_bbl2[0] eq -1 then sl_bbl2=bytarr(sl_bands2)+1

endif

; query 3rd spectral library
if sl_fid3 ne -1 then begin
    envi_file_query, sl_fid3, spec_names=spec_names3, nl=n_spec3, sname=speclib_sname3, ns=sl_bands3, bbl=sl_bbl3
    speclib3=envi_get_data(fid=sl_fid3, pos=0, dims=[-1,0,sl_bands3-1,0,n_spec3-1])
    ; if no bad bands list then make it
    if sl_bbl3[0] eq -1 then sl_bbl3=bytarr(sl_bands3)+1
endif

; if photo_shade is 0 then need to get non-photometric shade spectrum (npss)
if (*pstate).photo_shade eq 0 then begin
    sl_fid4=(*pstate).fid5
    envi_file_query, sl_fid4, spec_names=spec_names4, nl=n_spec4, sname=speclib_sname4, ns=sl_bands4, bbl=sl_bbl4
    if nsel4 gt 1 then begin
       result=dialog_message('Select only 1 non-zero shade spectrum')
        return
    endif
    speclib4=envi_get_slice(fid=(*pstate).fid5, line=selected4, pos=0, xs=0, xe=sl_bands4-1)
    ; if no bad bands list then make one
    if sl_bbl4[0] eq -1 then sl_bbl4=bytarr(sl_bands4)+1
    ; get shade name
    npss_name=spec_names4[selected4]
endif

;open first image file for testing number of bands and bad bands list
envi_open_file, infile[0], r_fid=image_fid, /no_realize
; query image file
envi_file_query, image_fid, nb=nb

; build lookup table that will store the spectra numbers for each model  each model up to 3 endmembers
lut=intarr(nmods,n_ems)
count=0
case n_ems of
    1: begin
        ; make sure same number of bands in file
        if nb ne sl_bands1 then begin
            result=dialog_message('Image file and spectral library have different number of bands')
            return
        end
        ; build common bad bands list
        bbl_all=sl_bbl1
        ; make lut
        lut=selected1
    end
    2: begin
    ; make sure same number of bands
        if (nb ne sl_bands1) or (nb ne sl_bands2) then begin
            result=dialog_message('Image file and 1 or more spectral libraries have different number of bands')
            return
        end
        ; build common bad bands list
        bbl_all=sl_bbl1*sl_bbl2
        ; make lut
        for i=0.0,nsel1-1 do begin
         for j=0.0, nsel2-1 do begin
             lut[count,*]=[selected1[i],selected2[j]]
             count=count+1
         endfor
       endfor
    end
    3: begin
        if (nb ne sl_bands1) or (nb ne sl_bands2) or (nb ne sl_bands3) then begin
           result=dialog_message('Image file and 1 or more spectral libraries have different number of bands')
            return
       end
       ; build common bad bands list
        bbl_all=sl_bbl1*sl_bbl2*sl_bbl3
        ; make lut
        for i=0.0,nsel1-1 do begin
         for j=0.0, nsel2-1 do begin
             for k=0.0, nsel3-1 do begin
                lut[count,*]=[selected1[i],selected2[j],selected3[k]]
                count=count+1
             endfor
         endfor
       endfor
    end
endcase

; begin loop for each input/output file **
for k=0L,nfiles-1 do begin

    ;get starting time for benchmarking
    starttime=systime(/seconds)

    widget_control, label1, set_value=strcompress('Processing image '+string(k+1)+' of '+string(nfiles))

    ; open input file
    envi_open_file, infile[k], r_fid=image_fid, /no_realize
    ; get output fraction image name
    out_image=outfile[k]
    ; set output class file name
    out_class_image=out_image+'_class'

    ; query image file
    envi_file_query, image_fid, bnames=bnames, byte_swap=byte_swap, fname=fname, nb=nb, interleave=interleave, $
        nl=nl, ns=ns, data_type=data_type, sname=image_sname, bbl=bbl, wl=wl
    image_dims=[-1,0,ns-1,0,nl-1]
    inherit=envi_set_inheritance(image_fid,image_dims,pos,/spatial)
    ;map_info = envi_get_map_info(fid=image_fid)

    ; deal with bad bands list
    if bbl[0] eq -1 then bbl=bytarr(nb)+1

    ; get the bands that are good in image and all spectral libraries
    pos=where(bbl*bbl_all eq 1)
    ngb=n_elements(pos)

    ; if residual constraint is being used, make sure number of goods bands exceeds the residual count
    if checks[4] ne 0 then begin
        if resid_count ge ngb then begin
           result=dialog_message('Residual count ('+strtrim(fix(resid_count),2)+') must be less than the number of good bands ('+strtrim(fix(ngb),2)+'), please fix and rerun')
           return
        endif
    endif

    ; subset to good bands
    speclib1=speclib1[pos,*]
    if n_ems gt 1 then speclib2=speclib2[pos,*]
    if n_ems gt 2 then speclib3=speclib3[pos,*]
    if (*pstate).photo_shade eq 0 then npss=speclib4[pos,*] ; check this in CRES

    ; warn user if image and speclibs have different bad bands selected
    t=array_equal(bbl,bbl_all)
    if t eq 0 then begin
       ;use compound widget warn user about bad band problems
       abase = widget_auto_base(title='Warning',/ybig)
       label = widget_label(abase, value='Image and at least 1 spectral library have different bad bands')
       label = widget_label(abase, value='selected, processing with good bands common to all files')
       label = widget_label(abase, value='')
       list = ['Yes', 'No']
       wm = widget_menu(abase, list=list, uvalue='menu', prompt='Output an ENVI compatible ASCII file with bad bands list?', /excl, /auto, default_ptr=1)
       result = auto_wid_mng(abase)
       ; if output bbl requested then write it
       if result.accept eq 1 and result.menu eq 0 then begin
         openw, u5, out_image+'_bbl.txt', /get_lun
         bbl_out=bbl*bbl_all
         for m=0,nb-1 do printf,u5,bbl_out[m]
         close,u5
         free_lun, u5
       endif
    endif

    ;open a metadata file that has the same name as the output image with a "_meta.txt" at the end and start filling it
    openw,u1,out_image+'_meta.csv', /GET_LUN
    printf,u1,'Metadata file for viper_mesma.pro run,',systime()
    printf,u1
    printf,u1,'Number of endmembers:,',n_ems+1
    printf,u1,'Input image:,', image_sname
    printf,u1,'Input image pathname:,',fname
    printf,u1,'Number of bands:,',nb
    printf,u1,'Number of good bands:,',ngb
    printf,u1,'Scale factor: ',scale
    if n_elements(speclib_sname1) ne 0 then printf,u1,format='(2(A,:,","))','1st input library:',speclib_sname1
    if ptr_valid((*pstate).selected1) then begin
       tmp=where(*(*pstate).selected1 eq 1)
       printf,u1,format='(1000(A,:,","))','Selected spectra:', tmp
    endif
    if n_elements(speclib_sname2) ne 0 then printf,u1,format='(2(A,:,","))','2nd input library:',speclib_sname2 else printf,u1,'2nd input library:,NA'
    if ptr_valid((*pstate).selected2) then begin
       tmp=where(*(*pstate).selected2 eq 1)
       printf,u1,format='(1000(A,:,","))','Selected spectra:', tmp
    endif else begin
       printf,u1, 'Selected spectra:,NA'
    endelse
    if n_elements(speclib_sname3) ne 0 then printf,u1,format='(2(A,:,","))','3rd input library:',speclib_sname3 else printf,u1,'3rd input library:,NA'
    if ptr_valid((*pstate).selected3) then begin
       tmp=where(*(*pstate).selected3 eq 1)
       printf,u1,format='(1000(A,:,","))','Selected spectra:', tmp
    endif else begin
       printf,u1, 'Selected spectra:,NA'
    endelse
    if (*pstate).photo_shade eq 1 then begin
        printf,u1,'Shade spectrum:, Photometric shade used'
    endif else begin
        printf,u1,'Shade spectrum:, Non-photometric shade spectrum from: ',speclib_sname4
        if ptr_valid((*pstate).selected5) then begin
         tmp=where(*(*pstate).selected4) eq 1
         printf,u1,format='(1000(A,:,","))','Selected spectra:', tmp-1
       endif
    endelse
    printf,u1,'Output SMA image:,       ',out_image
    printf,u1,'Number of lines:,        ', nl
    printf,u1,'Number of samples:,      ', ns
    printf,u1,'Output data type:,       ','BIL'
    printf,u1
    printf,u1,'Constraints used:
    if checks[0] ne 0 then printf,u1,'Minimum fraction (shade + nonshade):,',min_frac else printf,u1,'Minimum fraction (shade + nonshade):, NOT USED'
    if checks[1] ne 0 then printf,u1,'Maximum non-shade fraction:,',max_frac else printf,u1,'Maximum non-shade fraction:, NOT USED'
    if checks[2] ne 0 then printf,u1,'Maximum shade fraction:,',max_shade else printf,u1,'Maximum shade fraction:, NOT USED'

    if checks[2] alborz
jf;74nd$sDfgkjoi/UIFG&^%@R^&#BdNK}JDuC(*#BI(CB`NIE(*jnvf~n!ER#iGh^45J*Tbn_O98>Dg
  , NOT USED'
    if checks[3] ne 0 then printf,u1,'Maximum RMS error:,',max_rms else printf,u1,'Maximum RMS error:, NOT USED'
    if checks[4] ne 0 then begin
        printf,u1,'Maximum residual:,',max_resid
        printf,u1,'Residual count:,', resid_count
        endif else begin
        printf,u1,'Maximum residual:,','NOT USED'
        printf,u1,'Residual count:,','NOT USED'
    endelse
    printf,u1
    ; get the number of models that will be run.
    ;nmods = (nsel1 < 1) * (nsel2 < 1) * (nsel3 < 1)
    print,'number of models:,',nmods
    printf,u1,'Total number of models:,',nmods
    printf,u1
    printf,u1,'Model descriptions (endmembers): '

    ; print column headers
    column_headers=['model#','1st_EM','2nd_EM','3rd_EM','4th_EM']
    printf,u1,FORMAT='(5(A, :, ", "))',column_headers[0:n_ems+1]

    ; open output image file
    openw, u2, out_image, /get_lun
    ; if residual outputs selected, open resids file
    if outputs[2] eq 1 then openw, u3, out_image+'_resids', /get_lun

    ; begin a loop to read in image line by line
    for i=0L,nl-1 do begin
        ; retrieve a line of image data for good bands only
        image_line = envi_get_slice(fid=image_fid, line=i, pos=pos, xs=0, xe=ns-1)
        ; if non photometric shade is used, subtract shade spectrum from each spectra
        if (*pstate).photo_shade eq 0 then begin
           npss_line=transpose(rebin(npss,ngb,ns))
           image_line=image_line-npss_line
        endif
        ; set model number counter to 0
        mod_num=0
        ; create/reset the output floating point vectors/arrays for the new line of data
        case n_ems of
            1:out_frac=fltarr(ns)
            2:out_frac=fltarr(ns,2)
            3:out_frac=fltarr(ns,3)
        endcase
        out_shade=fltarr(ns)
        out_rmse=fltarr(ns)+1e20 ; fill out_rmse with very high initial rmse value
        out_model=fltarr(ns)
        ; Kerry Halligan added line of code out below on 12-21-06
        if outputs[2] eq 1 then out_resids=fltarr(ns,nb)

        ; begin j loop (for each model)
        for j=0L,nmods-1 do begin
            ; calculate model number
            mod_num=mod_num+1
            ; get the endmembers for model from speclib arrays using lut
            case n_ems of
                1:begin
                    em_mat=transpose(speclib1[*,lut[j]])
                    em_names=spec_names1[lut[j]]
                end
                2:begin
                    em_mat=transpose([[speclib1[*,lut[j,0]]],[speclib2[*,lut[j,1]]]])
                    em_names=([spec_names1[lut[j,0]],spec_names2[lut[j,1]]])
                    end
               3:begin
                    em_mat=transpose([[speclib1[*,lut[j,0]]],[speclib2[*,lut[j,1]]],[speclib3[*,lut[j,2]]]])
                    em_names=[spec_names1[lut[j,0]],spec_names2[lut[j,1]],spec_names3[lut[j,2]]]
                    end
            endcase
            ;subset em_mat to good bands from image if different
            if ngb ne n_elements(em_mat[0,*]) then begin
                    if ngb lt n_elements(em_mat[0,*]) then begin
                        em_mat=em_mat[*,pos]
                     result=dialog_message('Spectral libraries have fewer bands than image, subset with BBL from image?',/cancel)
                     if result eq 'Cancel' then return
                  endif else begin
                     result=dialog_message('Spectral libraries have more bands than image!',/error)
                     return
                 endelse
            endif


           ; if non-photometric shade is used, subtract from all endmembers
           if (*pstate).photo_shade eq 0 then em_mat=em_mat-npss_line[0:n_ems-1,*]

            ; on first line of image, print out all models and inputs to metadata file
            if i eq 0 then begin
                if (*pstate).photo_shade eq 0 then begin
                   printf,u1,FORMAT='(5(A, :, ", "))', string(mod_num),em_names,npss_name
             endif else begin
                 printf,u1,FORMAT='(5(A, :, ", "))', string(mod_num),em_names,'shade'
             endelse
            endif


            ;if j eq 2 then begin
            ;    print,'stop'
            ; endif

           ; ** CALCULATE SMA **
            ; call the function viper_mesma_fraccalc which calculates fractions and returns a structure with the following tags:
            ; frac, shade, rmse.  Also returns residuals if resids keyword is set
            sma = viper_mesma_fraccalc(em_mat, image_line, resids=resids)

;            ; for testing only
;             if i eq 10 then begin
;                print,sma.frac[590]
;                print,sma.shade[590]
;                print,sma.rmse[590]
;                print,'stop'
;                window,1
;                   wset,1
;                   plot,image_line[590,*]
;                oplot, em_mat, color=255
;                model=em_mat*sma.frac[590]
;                r2=image_line[590,*]-model
;             endif

           ; scale rmse values
           sma.rmse=sma.rmse/float(scale)

           ; ** COMPARE RMSE VALUES **
            ; first check to see if rmse returned by fraccalc is lower than current best models
            lower_rmse=sma.rmse lt out_rmse
            ; if any pixels have lower rmse then best models, then find pixels within selected constraints
            if total(lower_rmse) gt 0 then begin

                ; ** CALCULATE CONSTRAINTS **
                ; build and fill a contraints array that will store the binary results for each pixel
                ; 1 indicates it met a given constraint, 0 means it did not.  If constraint not set then
                ; set all pixels to 1 for that constraint
                con=bytarr(ns,5)
                ; pixels where non-shade em is greater than minimum fraction (e.g., -5%)
                if checks[0] eq 1 then begin
                    case n_ems of
                        ;2: con[*,0]=cmapply('*',sma.frac gt min_frac,2) ; 3em case
                        2: con[*,0]=(sma.frac[*,0] gt min_frac) * (sma.frac[*,1] gt min_frac)
                        ;3: con[*,0]=cmapply('*',sma.frac gt min_frac,2) ; 4em cases
                        3: con[*,0]=(sma.frac[*,0] gt min_frac) * (sma.frac[*,1] gt min_frac) * (sma.frac[*,2] gt min_frac)
                        1: con[*,0]=sma.frac gt min_frac ; 2em case
                    endcase
                endif else begin
                    con[*,0]=1
                endelse
                ; pixels where non-shade em is less than maximum fraction (e.g., 105%)
                if checks[1] eq 1 then begin
                    case n_ems of
                        ;2: con[*,1]=cmapply('*',sma.frac lt max_frac,2) ; 3em case
                        2: con[*,1]=(sma.frac[*,0] lt max_frac) * (sma.frac[*,1] lt max_frac)
                        ;3: con[*,1]=cmapply('*',sma.frac lt max_frac,2) ; 4em cases
                        3: con[*,1]=(sma.frac[*,0] lt max_frac) * (sma.frac[*,1] lt max_frac) * (sma.frac[*,2] lt max_frac)
                        1: con[*,1]=sma.frac lt max_frac ; 2em case
                    endcase
                endif else begin
                    con[*,1]=1
                endelse
                ; pixels where shade em fraction less than maximum shade fraction (e.g., 80%)
                if checks[2] eq 1 then con[*,2]=sma.shade lt max_shade else con[*,2]=1
                ; pixels where rmse is less than or equal to maximum allowable rms (e.g., )
                if checks[3] eq 1 then con[*,3]=sma.rmse lt max_rms else con[*,3]=1
                ; scale residuals
                resids=resids/float(scale)
                ; if residual thresholds set then find where residuals are exceeded
                if checks[4] eq 1 then begin
                    ex={value:max_resid, nbands:resid_count}
                    con[*,4]=cmapply('user:viper_mesma_thresh_resids', resids, 2, functargs=ex)
                endif else begin
                    con[*,4]=1
                endelse
                ; totals in con for each pixel should be 5 for all valid pixels
                valid=total(con,2) eq 5

                ;** UPDATE OUTPUT VECTORS **
                ; Get indices of all pixels that are both valid (meet criteria) and have lower
                ; RMSE than the current best
                update_samps=where(valid*lower_rmse eq 1, count)
                ;print, 'number of pixels updated: ',n_elements(update_samps)
                ; if there are pixels to update, then do so.
                if update_samps[0] ne -1 then begin
                    out_frac[update_samps,*]=sma.frac[update_samps,*]
                    out_shade[update_samps,*]=sma.shade[update_samps,*]
                    out_rmse[update_samps,*]=sma.rmse[update_samps,*]
                    out_model[update_samps,*]=mod_num
                    ; added line below on 12-20-06, modified on 1-1-07
                    if outputs[2] eq 1 then for l=0,ngb-1 do out_resids[update_samps,pos[l]]=resids[update_samps,l]
                endif

            ; end the lower RMSE if statement
            endif
        ; end j loop (model by model loop)
        endfor

        ; reset RMSE for all unmodeled pixels to 0
        index=where(out_rmse eq 1e20)
        if index[0] ne -1 then out_rmse[where(out_rmse eq 1e20)]=0
        ; write output vectors to output SMA image (floating point format)
        writeu, u2, [[out_frac],[out_shade],[out_rmse],[out_model]]
        ; if residual outputs selected, write resids to file
        if outputs[2] eq 1 then begin
         ; Kerry Halligan commented lines of code out below on 12-20 leave residuals in floating point reflectance/dn/radiance
         ; units regardless of input scaling
                ; rescale the residuals and set to original data type
         ;            resids=resids*float(scale)
         ;            case data_type of
         ;                1: resids=byte(resids)
         ;                2: resids=fix(resids)
         ;                4: resids=float(resids)
         ;                else: resids=float(resids)
         ;            endcase
         ; Kerry Halligan added line of code on 12-21-06 below to output the final residuals not the last residuals
            writeu, u3, out_resids
            ;writeu, u3, resids
        endif
        ; update progress bar
        completed=fix(float(i)/float(nl)*100)
        widget_control, slider, bad_id=bi
        if bi ne 0 then return
        widget_control, slider, set_value=completed

    ; end i loop (image line by line loop)
    endfor

    ;close output file(s)
    close, u2
    free_lun, u2
    if outputs[2] eq 1 then begin
        close, u3
        free_lun, u3
    endif

    ;make the array to hold the band names for minrms image
    bnames=['EM1_fraction','EM2_fraction','EM3_fraction']
    bnames=[bnames[0:n_ems-1],'shade_fraction','RMSE','Model#']

    ;write the ENVI header for output SMA image and open file
    if !version.os_family eq 'Windows' then begin
        envi_setup_head, data_type=4, descrip='MESMA image', file_type=0, fname=out_image+'.hdr', $
            interleave=1, nb=out_nb, nl=nl, ns=ns, bnames=bnames, /write, inherit=inherit;, map_info=map_info,
    endif else begin
        envi_setup_head, data_type=4, descrip='MESMA image', file_type=0, fname=out_image, $
            interleave=1, nb=out_nb, nl=nl, ns=ns, bnames=bnames, /write, inherit=inherit;, map_info=map_info,
    endelse
    envi_open_file, out_image, r_fid=r_fid, /no_realize

    if outputs[2] eq 1 then begin
    ; Kerry Halligan changed data type to 4 from input data type (data_type) on 12-20
        if !version.os_family eq 'Windows' then begin
            ;write the ENVI header for output residual image, if selected, and open file
            envi_setup_head, data_type=4, descrip='MESMA residuals image', file_type=0, $
               fname=out_image+'_resids.hdr', interleave=1, nb=nb, wl=wl, $
                nl=nl, ns=ns, /write, inherit=inherit;, map_info=map_info,
        endif else begin
           ;write the ENVI header for output residual image, if selected, and open file
               envi_setup_head, data_type=4, descrip='MESMA residuals image', file_type=0, $
                      fname=out_image+'_resids', interleave=1, nb=nb, wl=wl, $
                   nl=nl, ns=ns, /write, inherit=inherit;, map_info=map_info,
       endelse
        if open_files eq 1 then envi_open_file, out_image+'_resids', /no_realize
    endif

    ; query image file and read in model number band
    envi_file_query, r_fid, interleave=rinterleave, data_type=data_type, nb=rnb, ns=rns, nl=rnl
    ; read in the final models from the output file
    mods = long(envi_get_data(fid=r_fid, dims=image_dims, pos=rnb-1))

    ; calculate the number of unclassified pixels
    npix=total(mods eq 0)
    percent=float(npix)/float(rnl*rns)*100
    printf,u1,' '
    printf,u1,'Unclassified pixels: '
    printf,u1,FORMAT='(2(A, :, ", "))','#pixels','% of image'
    printf,u1,FORMAT='(I,A,F)',npix,',',percent

    ; remove unclassified models
    suc_index=where(mods gt 0,count)
    if suc_index[0] ne -1 then begin
        suc_mods=mods[suc_index]
        ; calculate the successful models and the number of pixels mapped by each model
        printf,u1,' '
        uniq_mods = suc_mods[UNIQ(suc_mods, SORT(suc_mods))]
        n_uniq_mods = n_elements(uniq_mods)
        printf,u1,' '
        printf,u1,'Number of successful models:,',n_uniq_mods
        printf,u1
        printf,u1,'List of successful models:'
        printf,u1,FORMAT='(1000(I, :, ", "))', uniq_mods
        printf,u1,' '
        printf,u1,'Number of pixels mapped for each model:'
        printf,u1,FORMAT='(3(A, :,", "))','model#','#pixels','% of image'
        for i=0,n_uniq_mods-1 do begin
            npix=total(suc_mods eq uniq_mods[i])
            percent=float(npix)/float(rnl*rns)*100
            printf,u1,FORMAT='(I,A,I,A,F)',long(uniq_mods[i]),',',npix,',',percent
        endfor
        printf,u1,' '
    endif else begin
        printf,u1, 'No models succesful'
    endelse

    ; if classification image selected, create it
    if outputs[1] eq 1 then begin
        ; create an array to hold the classnames for the classified image
        classnames=['Unclassified',strcompress('model_'+string(indgen(nmods)+1),/remove_all)]
        ; create a lookup table for class colors, based on the RAINBOW color tabel, resized
        ; to fit the number of models + 1 for unclassified class
        loadct,13; rainbow
        tvlct, r, g, b, /get
        c_lut=transpose([[r],[g],[b]])
        num_classes=nmods+1
        color_lut=congrid(c_lut,3,num_classes,/minus_one, /interp)
        ; open up third output file and write out model numbers to the classified image
        openw, u4, out_class_image, /get_lun
        writeu, u4, mods
        close, u4
        free_lun, u4
        ;write the ENVI header for minrms classification and open file
        if !version.os_family eq 'Windows' then begin
            envi_setup_head, data_type=3, descrip='MESMA classification', file_type=3, fname=out_class_image+'.hdr', interleave=0, nb=1, $
                nl=nl, ns=ns, bnames='MESMA_class', class_names=classnames, lookup=color_lut, num_classes=num_classes, $
                /write, inherit=inherit;, map_info=map_info,
       endif else begin
            envi_setup_head, data_type=3, descrip='MESMA classification', file_type=3, fname=out_class_image, interleave=0, nb=1, $
                   nl=nl, ns=ns, bnames='MESMA_class', class_names=classnames, lookup=color_lut, num_classes=num_classes, $
                   /write, inherit=inherit;, map_info=map_info,
       endelse
        if open_files eq 1 then envi_open_file, out_class_image, /no_realize
    endif

    print, 'Done with file'+i+1+' of '+nfiles+1
    endtime=systime(/seconds)
    elapsed_time=float(endtime-starttime)/60.0
    print,'Elapsed time: ',elapsed_time,' minutes'
    printf,u1,'Elapsed time: ',elapsed_time,' minutes'
    close,u1
    free_lun,u1
    ; if open files selected, then open up metadata file
    if open_files eq 1 then begin
       if !version.os_family eq 'Windows' then begin
            spawnstring='notepad '+out_image+'_meta.csv'
            spawn, spawnstring, /nowait, /noshell
       endif
    endif

    ; close sma file if requested
    if open_files eq -1 then begin
        envi_file_mng, id=r_fid, /remove
    endif
    ; end the loop for multiple input/output files (k loop)
    endfor
;destroy GUI and build progress bar
widget_control, base, /destroy
print,'done'
end



; routine to launch the help index html file
; NOT CURRENTLY IMPLIMENTED!!
pro viper_mesma_help, event
compile_opt idl2
on_error, 2

; find save_add directory for ENVI distribution
save_add_dir = FILE_SEARCH(!DIR, 'save_add',/test_directory)

; check to make sure index.html exists
filename = FILE_SEARCH(save_add_dir, 'viper_mesma.html')

if strlen(filename) ne 0 then begin
    online_help, book=filename, /full_path
endif else begin
    result=dialog_message('viper_mesma help file not found!')
endelse

; end
end



; This proceedure reads in the pstate variable and gathers info from the GUI
; and writes out all the important data to a control file (*.csv) to allow
; for efficient re-running of viper_mesma
pro viper_mesma_save, event
COMPILE_OPT idl2, HIDDEN
;on_error, 2

; Establish error handler
catch, theError
if theError ne 0 then begin
    catch, /cancel
    help, /last_message, output=errText
    errMsg = dialog_message(errText, /error, title='Error processing request')
    return
endif

; Retrieve the pointer to the state structure
widget_control, event.top, get_uvalue = pstate

; select output filename
;outfile=(*pstate).outfile
;if ptr_valid(outfile) then outfile=*outfile
;outfile=outfile[0]
widget_control, (*pstate).file4text, get_value=outfile
if strlen(outfile) ne 0 then begin
    gp=file_dirname(outfile)
    sname=file_basename(outfile)
    tmp=strsplit(sname,'.',/extract)
    tmp=tmp[0]+'_control.csv'
endif else begin
    gp='' & tmp=''
endelse
outfile = dialog_pickfile(title='Enter output filename (*.csv)', path=gp, file=tmp, filter='*.csv', /write, /overwrite_prompt)
if outfile eq '' then return
;outfile=dialog_pickfile(title='Selected file to write',path=(*pstate).inpath, filter='*.csv',/write, /overwrite_prompt);ctl

if strlen(outfile) eq 0 then return
openw,u3,outfile,/get_lun

printf,u3,'viper_mesma.pro control file created:, ',systime()

printf,u3
printf,u3,'Inputs:'

if ptr_valid((*pstate).infile) then begin
    infile = *(*pstate).infile
    printf,u3,format='(1000(A,:,","))','Input image file(s):',infile
endif else begin
printf,u3,'Input image file(s):,-1'
endelse
inpath=(*pstate).inpath
if inpath ne '' then printf,u3,format='(2(A,:,","))','Inpath image file path:',(*pstate).inpath else printf,u3,'Inpath image file path:,-1'
scale=widget_info((*pstate).scale_drop, /droplist_select)
printf,u3,'Scale factor:(0=not selected; 1=1; 2=1000; 2=10000),',scale

speclib=strtrim((*pstate).file1fulltext,2)
if speclib ne '' then printf,u3,format='(2(A,:,","))','1st input library:',speclib else printf,u3,'1st input library:,-1'
if ptr_valid((*pstate).selected1) then printf,u3,format='(1000(A,:,","))','Selected spectra:',*(*pstate).selected1 else printf,u3,'Selected spectra:,-1'

speclib=strtrim((*pstate).file2fulltext, 2)
if speclib ne '' then printf,u3,format='(2(A,:,","))','2nd input library:',speclib else printf,u3,'2nd input library:,-1'
if ptr_valid((*pstate).selected2) then printf,u3,format='(1000(A,:,","))','Selected spectra:',*(*pstate).selected2 else printf,u3,'Selected spectra:,-1'

speclib=strtrim((*pstate).file3fulltext, 2)
if speclib ne '' then printf,u3,format='(2(A,:,","))','3rd input library:',speclib else printf,u3,'3rd input library:,-1'
if ptr_valid((*pstate).selected3) then printf,u3,format='(1000(A,:,","))','Selected spectra:',*(*pstate).selected3 else printf,u3,'Selected spectra:,-1'

printf,u3
printf,u3,'Contraints used (0 = not used; 1 = used):'
printf,u3,'Minimum allowable non-shade fraction:, ',(*pstate).checks[0]
if (*pstate).checks[0] eq 1 then printf,u3, 'Value:, ',(*pstate).constraints[0] else printf,u3, 'Value:, -1'
printf,u3,'Maximum allowable non-shade fraction:, ',(*pstate).checks[1]
if (*pstate).checks[1] eq 1 then printf,u3, 'Value:, ',(*pstate).constraints[1] else printf,u3, 'Value:, -1'
printf,u3,'Maximum allowable shade fraction:, ',(*pstate).checks[2]
if (*pstate).checks[2] eq 1 then printf,u3, 'Value:, ',(*pstate).constraints[2] else printf,u3, 'Value:, -1'

printf,u3,'Maximum allowable RMSE:, ',(*pstate).checks[3]
if (*pstate).checks[3] eq 1 then begin
    widget_control, (*pstate).maxrmse_text, get_value=value
    printf,u3, 'Value:, ',string(float(value))
endif else begin printf,u3, 'Value:, -1'
endelse

printf,u3,'Residual threshold:, ',(*pstate).checks[4]
if (*pstate).checks[4] eq 1 then begin
    widget_control, (*pstate).residv_text, get_value=value
    printf,u3, 'Value:, ',string(float(value))
endif else begin printf,u3, 'Value:, -1'
endelse

if (*pstate).checks[4] eq 1 then printf,u3, 'Number of bands:, ',(*pstate).constraints[5] else printf,u3, 'Number of bands:, '
; save information on the shade to use
printf,u3,'Photometric shade? (1=photometric; 0=non-photometric shade):,',(*pstate).photo_shade
speclib=strtrim((*pstate).file5fulltext, 2)
if speclib ne '' then printf,u3,format='(2(A,:,","))','Non-photometric shade spectral library:',speclib else printf,u3, 'Non-photometric shade spectral library:, -1'
if speclib ne '' then printf,u3,format='(2(A,:,","))','Selected spectrum:',(*pstate).selected5 else printf,u3,'Selected spectrum:,-1'

; write out the outputs to produce and if they should be opened in ENVI
printf,u3
printf,u3,'Outputs (0=not selected; 1=selected):'
if ptr_valid((*pstate).outfile) then begin
    outfile = *(*pstate).outfile
    ; if in single file mode, then get the output filename from the widget rather than pstate
    if (*pstate).multiple eq 0 then widget_control, (*pstate).file4text, get_value=outfile
    printf,u3,format='(1000(A,:,","))','Output file(s):',outfile
endif else begin
printf,u3,'Output file(s):,-1'
endelse
printf,u3,'Fractions, ',(*pstate).outputs[0]
printf,u3,'Classification, ',(*pstate).outputs[1]
printf,u3,'Residuals,' ,(*pstate).outputs[2]
printf,u3,'Open outputs, ',(*pstate).open_files

; let user know that it worked
result=dialog_message('Control file saved',/information)

; close files and end
close,u3
free_lun, u3
end

GUI up to 3 spectral libraries

; This proceedure reads a control file and updates the GUI and the pstate
; structure to reflect stored settings for efficient re-running of viper_mesma
pro viper_mesma_restore, event
COMPILE_OPT idl2, HIDDEN
;on_error, 2

; Establish error handler
catch, theError
if theError ne 0 then begin
    catch, /cancel
    help, /last_message, output=errText
    errMsg = dialog_message(errText, /error, title='Error processing request')
    return
endif

; Retrieve the pointer to the state structure
widget_control, event.top, get_uvalue = pstate

; get file and retrieve data
infile=dialog_pickfile(title='Selected Control File to Load',path=(*pstate).inpath, filter='*.csv',/read)
if strlen(infile) eq 0 then return
nrows=file_lines(infile)
line=''
ctl=strarr(nrows)
openr,u1,infile,/get_lun
for i=0,nrows-1 do begin
    readf,u1,line
    ctl[i]=line
endfor

; get images
; get input filename(s) by splitting string at commas if present
tmp=strsplit(ctl[3],',',/extract)
if tmp[1] ne '-1' then begin
    n=n_elements(tmp)
    inputfile=strtrim(tmp[1:n-1],2)

    ; verify that first file is valid
    result=file_info(inputfile[0])
    if result.exists eq 1 then begin

        (*pstate).inpath=file_dirname(inputfile[0])
        (*pstate).infile=ptr_new(inputfile)

        ; find out if first image at least is valid ENVI file
        envi_open_file, *(*pstate).infile[0], r_fid=fid, /no_realize
        envi_file_query, fid, ns=ns, file_type=filetype, nl=nl, nb=nb, sname=sname
        (*pstate).fid0=fid
        ; set the text widget to display the selected file
        widget_control, (*pstate).file0text, set_value=sname
        slidemax= 20 < nb
        widget_control, (*pstate).residb_slider, SET_SLIDER_MAX=slidemax

        (*pstate).fid0=fid  ; ****************** this won't work for multiple files
    endif else result=dialog_message('Could not locate image file: '+inputfile[0])
endif

;set input path
tmp=strsplit(ctl[4],',',/extract)
if tmp[1] ne '-1' then (*pstate).inpath=strtrim(ctl[4],2)

; get and set the scale
;tmp=strsplit(ctl[5],',',/extract)
;(*pstate).scale=tmp[1]
;widget_control, (*pstate).scale_drop, set_droplist_select=(*pstate).scale

; get and set the scale
list=[0,1,1000,10000]
tmp=strsplit(ctl[5],',',/extract)
(*pstate).scale=list[tmp[1]]
widget_control, (*pstate).scale_drop, set_droplist_select=tmp[1]

; 1st spectral library
tmp=strtrim(strsplit(ctl[6],',',/extract),2)
if tmp[1] ne '-1' then begin
    ; verify that file is valid
    result=file_info(tmp[1])
    if result.exists eq 1 then begin
        envi_open_file, tmp[1], r_fid=fid, /no_realize
        envi_file_query, fid, spec_names=specnames
        (*pstate).specnames1=ptr_new(specnames)
        (*pstate).fid1=fid
        file=file_basename(tmp[1])
        widget_control, (*pstate).file1text, set_value=file
    endif else result=dialog_message('Could not locate spectral library: '+tmp[1])
endif

; 2nd spectral library
tmp=strtrim(strsplit(ctl[8],',',/extract),2)
if tmp[1] ne '-1' then begin
    ; verify that file is valid
    result=file_info(tmp[1])
    if result.exists eq 1 then begin
        envi_open_file, tmp[1], r_fid=fid, /no_realize
        envi_file_query, fid, spec_names=specnames
        (*pstate).specnames2=ptr_new(specnames)
        (*pstate).fid2=fid
        file=file_basename(tmp[1])
        widget_control, (*pstate).file2text, set_value=file
    endif else result=dialog_message('Could not locate spectral library: '+tmp[1])
endif

; 3rd spectral library
tmp=strtrim(strsplit(ctl[10],',',/extract),2)
if tmp[1] ne '-1' then begin
    ; verify that file is valid
    result=file_info(tmp[1])
    if result.exists eq 1 then begin
        envi_open_file, tmp[1], r_fid=fid, /no_realize
        envi_file_query, fid, spec_names=specnames
        (*pstate).specnames3=ptr_new(specnames)
        (*pstate).fid3=fid
        file=file_basename(tmp[1])
        widget_control, (*pstate).file3text, set_value=file
    endif else result=dialog_message('Could not locate spectral library: '+tmp[1])
endif

; set the selected spectra
tmp=strsplit(ctl[7],',',/extract)
if tmp[1] ne '-1' then begin
    n=n_elements(tmp)
    sel=fix(tmp[1:n-1])
    if n_elements(sel) eq 0 then result=dialog_message('No spectra selected in nonphotometric shade spectral library')
    (*pstate).selected1=ptr_new(sel)
endif

; set the selected spectra
tmp=strsplit(ctl[9],',',/extract)
if tmp[1] ne '-1' then begin
    n=n_elements(tmp)
    sel=fix(tmp[1:n-1])
    if n_elements(sel) eq 0 then result=dialog_message('No spectra selected in nonphotometric shade spectral library')
    (*pstate).selected2=ptr_new(sel)
endif

; set the selected spectra
tmp=strsplit(ctl[11],',',/extract)
if tmp[1] ne '-1' then begin
    n=n_elements(tmp)
    sel=fix(tmp[1:n-1])
    if n_elements(sel) eq 0 then result=dialog_message('No spectra selected in nonphotometric shade spectral library')
    (*pstate).selected3=ptr_new(sel)
endif

info=strarr(2,11)
for i=0,10 do info[*,i]=strsplit(ctl[i+14],',',/extract)
info=float(info[1,*])
list=[0,2,4,6,8]
(*pstate).checks=reform(info[list])
list=[1,3,5,7,9,10]
(*pstate).constraints=reform([info[list]])

; check appropriate checkboxes
widget_control, (*pstate).check1, set_button=(*pstate).checks[0]
widget_control, (*pstate).check2, set_button=(*pstate).checks[1]
widget_control, (*pstate).check3, set_button=(*pstate).checks[2]
widget_control, (*pstate).check4, set_button=(*pstate).checks[3]
widget_control, (*pstate).check5, set_button=(*pstate).checks[4]

; set sensitivity and values for sliders
if (*pstate).checks[0] eq 1 then begin
    widget_control, (*pstate).min_slider, set_value=(*pstate).constraints[0]
    widget_control, (*pstate).p1_slide, sensitive=1
    endif else begin
    widget_control, (*pstate).p1_slide, sensitive=0
    endelse
if (*pstate).checks[1] eq 1 then begin
    widget_control, (*pstate).max_slider, set_value=(*pstate).constraints[1]
    widget_control, (*pstate).p2_slide, sensitive=1
    endif else begin
    widget_control, (*pstate).p2_slide, sensitive=0
    endelse
if (*pstate).checks[2] eq 1 then begin
    widget_control, (*pstate).shade_slider, set_value=(*pstate).constraints[2]
    widget_control, (*pstate).p3_slide, sensitive=1
    endif else begin
    widget_control, (*pstate).p3_slide, sensitive=0
    endelse
if (*pstate).checks[3] eq 1 then begin
    widget_control, (*pstate).maxrmse_text, set_value=string((*pstate).constraints[3]), sensitive=1
    widget_control, (*pstate).maxrmse_lab, sensitive=1
    endif else begin
    widget_control, (*pstate).maxrmse_text, sensitive=0
    widget_control, (*pstate).maxrmse_lab, sensitive=0
    endelse
if (*pstate).checks[4] eq 1 then begin
    widget_control, (*pstate).residv_text, set_value=string((*pstate).constraints[4]), sensitive=1
    widget_control, (*pstate).residv_lab, sensitive=1
    widget_control, (*pstate).residb_slider, set_value=round((*pstate).constraints[5]), sensitive=1
    endif else begin
    widget_control, (*pstate).residv_text, sensitive=0
    widget_control, (*pstate).residv_lab, sensitive=0
    widget_control, (*pstate).residb_slider, sensitive=0
    endelse

filelist=[(*pstate).fid1, (*pstate).fid2, (*pstate).fid3]
sellist=lonarr(3)
if ptr_valid((*pstate).selected1) then sellist[0]=total(*(*pstate).selected1)
if ptr_valid((*pstate).selected2) then sellist[1]=total(*(*pstate).selected2)
if ptr_valid((*pstate).selected3) then sellist[2]=total(*(*pstate).selected3)
;sellist=[total(*(*pstate).selected1), total(*(*pstate).selected2), total(*(*pstate).selected3)]
index=where(filelist gt -1 and sellist gt 0)
if index[0] ne -1 then begin
    sellist=sellist[index]
    nmods=cmapply('*',sellist)
endif else begin
    nmods=0L
endelse
widget_control, (*pstate).total_txt, set_value=string(long(nmods))

; get the type of shade and if non-photometric shade, set the filename and the selected spectra
tmp=strsplit(ctl[26],',',/extract)
if tmp[1] eq 0 then begin
    (*pstate).photo_shade=0
    widget_control, (*pstate).nonphoto_shade_but, set_button=1
    widget_control, (*pstate).file5text, sensitive=1
    widget_control, (*pstate).sel5lab, sensitive=1
    widget_control, (*pstate).file5browse, sensitive=1
    widget_control, (*pstate).file5subset, sensitive=1

    ; find out if image at least is valid ENVI file
    envi_open_file, tmp[1], r_fid=fid, /no_realize
    envi_file_query, fid, spec_names=specnames
    (*pstate).specnames5=ptr_new(specnames)
    (*pstate).fid5=fid
    ; write input filename(s) to text widget
    file=file_basename(tmp[1])
    widget_control, (*pstate).file5text, set_value=file
    ; get the selected spectra
    tmp=strsplit(ctl[27],',',/extract)
    if tmp[1] ne '-1' then begin
        (*pstate).selected5=fix(tmp[1])
        widget_control, (*pstate).file5subset, set_value=['Select spectrum', specnames]
        widget_control, (*pstate).file5subset, set_droplist_select=fix(tmp[1])
    endif else result=dialog_message('No spectra selected in nonphotometric shade spectral library')
    print,'stop'
endif

; get and set the outputs to produce and if they should be opened in ENVI
info=strarr(2,4)
for i=0,3 do info[*,i]=strsplit(ctl[i+31],',',/extract)
info=fix(info[1,*])
widget_control, (*pstate).out_frac, set_button=info[0]
widget_control, (*pstate).out_class, set_button=info[1]
widget_control, (*pstate).out_resid, set_button=info[2]
widget_control, (*pstate).open_files_but, set_button=info[3]

; get output filename(s) by splitting string at commas if present
tmp=strsplit(ctl[30],',',/extract)
; convert to pointer and store
if tmp[1] ne '-1' then begin
    n=n_elements(tmp)
    outfile=tmp[1:n-1]
    (*pstate).outfile=ptr_new(outfile)
    ; update text widget
    widget_control, (*pstate).file4text, set_value=*(*pstate).outfile
endif

; let user know that it worked
result=dialog_message('Control file restored',/information)

; close file and end
close,u1
free_lun, u1
end


; The cleanup routine
pro viper_mesma_cleanup, top
COMPILE_OPT idl2, HIDDEN
on_error, 2
    ; Get the state variable from the top-level base.
    widget_control, top, get_uvalue=pstate
    ; Destroy the state variable.
    ptr_free, pstate
end

; the main GUI routine for viper_mesma (see main header above for details)
pro viper_mesma, event
compile_opt idl2
on_error, 2

; make top-level base
tlb = widget_base( title='VIPER Tools: SMA/MESMA', /col, xoffset=300, yoffset=100)

graphic_base=widget_base(tlb,/align_center)
graphic=widget_draw(graphic_base, xsize=202, ysize=52)

; create a main label
;guilab_base=widget_base(tlb, /align_center)
;guilab = widget_label(guilab_base, value='VIPER Tools: MESMA')

; create base to hold all widgets
control_base = widget_base(tlb, /col)
; create base to hold all data input widgets
in_base = widget_base(control_base, /row, frame=2)
; create base to hold all data input widgets
file_base = widget_base(in_base, /col);, frame=2)

; create a row base to hold input file information
file0base = widget_base(file_base, /row)
; create a label for 1st file
sel0lab = widget_label(file0Base, value='Image(s) to unmix:', xsize=150)
; create a text widget to display file once selected
file0text = widget_text(file0base, value='Required', xsize=50)
; create browse button to launch file selection routine
file0browse = widget_button(file0base, value='Browse', event_pro='viper_mesma_browse', xsize=50)

; create scale factor dropdown list
scale_drop=widget_droplist(file0Base, value=['scale factor','1 / DN / Radiance','1,000','10,000'], $
    event_pro='viper_mesma_scale',xsize=115)

; create a row base to hold 1st file information
file1base = widget_base(file_base, /row)
; create a label for 1st file
sel1lab = widget_label(file1Base, value='1st Library:', xsize=150)
; create a text widget to display file once selected
file1text = widget_text(file1base, value='Required', xsize=50)
; create browse button to launch file selection routine
file1browse = widget_button(file1base, value='Browse', event_pro='viper_mesma_browse', xsize=50)
; create dropdown list for subsetting of spectral library
;file1subbase = widget_base(file1base, /row, /nonexclusive)
file1subset = widget_button(file1base, value='Subset spectra', $
    event_pro='viper_mesma_subset',xsize=100)

; create a row base to hold 2st file information
file2base = widget_base(file_base, /row)
; create a label for 2st file
sel2lab = widget_label(file2Base, value='2nd Library:', xsize=150)
; create a text widget to display file once selected
file2text = widget_text(file2base, value='None', xsize=50)
; create browse button to launch file selection routine
file2browse = widget_button(file2base, value='Browse', event_pro='viper_mesma_browse', xsize=50)
; create dropdown list for subsetting of spectral library
file2subset = widget_button(file2base, value='Subset spectra', event_pro='viper_mesma_subset', xsize=100)

; create a row base to hold 3st file information
file3base = widget_base(file_base, /row)
; create a label for 1st file
sel3lab = widget_label(file3Base, value='3rd Library:', xsize=150)
; create a text widget to display file once selected
file3text = widget_text(file3base, value='None', xsize=50)
; create browse button to launch file selection routine
file3browse = widget_button(file3base, value='Browse', event_pro='viper_mesma_browse', xsize=50)
; create dropdown list for subsetting of spectral library
file3subset = widget_button(file3base, value='Subset spectra', event_pro='viper_mesma_subset', xsize=100)

; create a row base to hold 3st file information
shade_tog_base = widget_base(file_base, /row, /exclusive, /align_center)
photo_shade_but = widget_button(shade_tog_base, value='Photometric shade', event_pro='viper_mesma_shade')
nonphoto_shade_but = widget_button(shade_tog_base, value='Non-zero shade', event_pro='viper_mesma_shade')
widget_control, photo_shade_but, /set_button

; create a row base to hold 3st file information
file5base = widget_base(file_base, /row)
; create a label for 1st file
sel5lab = widget_label(file5Base, value='Non-zero Shade Library:', xsize=150, sensitive=0)
; create a text widget to display file once selected
file5text = widget_text(file5base, value='None', xsize=50, sensitive=0)
; create browse button to launch file selection routine
file5browse = widget_button(file5base, value='Browse', event_pro='viper_mesma_browse', xsize=50, sensitive=0)
; create dropdown list for subsetting of spectral library
file5subset = widget_droplist(file5base, value='Select spectrum', event_pro='viper_mesma_subset', xsize=115, sensitive=0)

; create widget to show number of models to be run
total_base=widget_base(in_base,/col,/align_center);, frame=2)
total_lab=widget_label(total_base, value='Total number')
total_lab=widget_label(total_base, value='of models:')
total_txt=widget_text(total_base, xsize=15, value='')

; create a base to hold all sliders
frame_base=widget_base(control_base, /row, event_pro='viper_mesma_slider')
param_base1=widget_base(frame_base, /col, frame=1)
frac_lab=widget_label(param_base1, value='Fraction constraints')
frame2_base=widget_base(frame_base, /col)
param_base2=widget_base(frame2_base, /col, frame=1)
rmse_lab=widget_label(param_base2, value='RMSE constraint')
param_base3=widget_base(frame2_base, /col, frame=1)
resids_lab=widget_label(param_base3, value='Residual constraints')
p1_base=widget_base(param_base1, /row)
p1_check=widget_base(p1_base, /nonexclusive)
check1=widget_button(p1_check,value='', event_pro='viper_mesma_check')
p1_slide_base=widget_base(p1_base, /col)
lab=widget_label(p1_slide_base, value='Minimum Allowable EM Fraction')
p1_slide=widget_base(p1_slide_base, /row, sensitive=0)
min_slider = cw_fslider(p1_slide, xsize=315, minimum=-.50, maximum=1.50, value=-.05, format='(F12.2)', scroll=0.01)
p2_base=widget_base(param_base1, /row)
p2_check=widget_base(p2_base, /nonexclusive)
check2=widget_button(p2_check,value='', event_pro='viper_mesma_check')
p2_slide_base=widget_base(p2_base,/col)
lab=widget_label(p2_slide_base, value='Maximum Allowable EM Fraction')
p2_slide=widget_base(p2_slide_base, /row, sensitive=0)
max_slider = cw_fslider(p2_slide, xsize=315, minimum=-.50, maximum=1.50, value=1.05, format='(F12.2)', scroll=0.01)
p3_base=widget_base(param_base1, /row)
p3_check=widget_base(p3_base, /nonexclusive)
check3=widget_button(p3_check,value='', event_pro='viper_mesma_check')
p3_slide_base=widget_base(p3_base, /col)
lab=widget_label(p3_slide_base, value='Maximum Allowable Shade Fraction')
p3_slide=widget_base(p3_slide_base, /row, sensitive=0)
shade_slider = cw_fslider(p3_slide, xsize=315, minimum=0, maximum=1.00, value=.80, format='(F12.2)')
p4_base=widget_base(param_base2, /row)
p4_check=widget_base(p4_base, /nonexclusive)
check4=widget_button(p4_check,value='', event_pro='viper_mesma_check')
p4_slide=widget_base(p4_base, /row)
p5_base=widget_base(param_base3, /row)
p5_check=widget_base(p5_base, /nonexclusive)
check5=widget_button(p5_check,value='',event_pro='viper_mesma_check')
p5_slide=widget_base(p5_base, /row)
p6_base=widget_base(param_base3, /row)
spacer=widget_base(p6_base, xsize=37)
p6_slide=widget_base(p6_base, /row)

; text widgets for input of rmse and residual value
maxrmse_lab=widget_label(p4_slide, value='Maximum Allowable RMSE', sensitive=0, xsize=260)
maxrmse_text = widget_text(p4_slide, value='0.025', /editable, $
    xsize=10, event_pro='viper_mesma_text', sensitive=0)
residv_lab=widget_label(p5_slide, value='Residual Threshold', sensitive=0, xsize=260)
residv_text = widget_text(p5_slide, value='0.025', /editable, $
    xsize=10, event_pro='viper_mesma_text', sensitive=0)

residb_slider = widget_slider(p6_slide, title='Number of Contiguous Bands', $
    xsize=315, minimum=0, maximum=20, value=7, sensitive=0)

; create a row base to hold 1st file information
out_frame=widget_base(control_base, /col, /align_center, frame=2)
out_base = widget_base(out_frame, /col, /align_center)
file4base = widget_base(out_base, /row, /align_center)
; create a label for 3st file
sel4lab = widget_label(file4Base, value='Output Fraction Image(s):', xsize=150)
; create a text widget to display file once selected
file4text = widget_text(file4base, value='', xsize=50);, /editable)
; create browse button to launch file selection routine
file4browse = widget_button(file4base, value='Browse', event_pro='viper_mesma_browse', xsize=50)

; create base to hold output checkboxes
output_base=widget_base(out_frame, /row, /align_center, /nonexclusive);, frame=1)
out_frac = widget_button(output_base, value='Fraction/RMSE Image', event_pro='viper_mesma_out')
out_class = widget_button(output_base, value='Classification Image', event_pro='viper_mesma_out')
out_resid = widget_button(output_base, value='Residuals Image', event_pro='viper_mesma_out')
widget_control, out_frac, /set_button, sensitive=0
widget_control, out_class, /set_button

open_files_base=widget_base(out_frame, /row, /align_center, /nonexclusive)
open_files_but = widget_button(open_files_base, value='Open output files in ENVI for viewing', event_pro='viper_mesma_open')
widget_control, open_files_but, /set_button

; create run button
run_base = widget_base(control_base, /row, /align_center)
run_but = widget_button(run_base, value='Run', tooltip='', $
    xsize=120, ysize=30, event_pro='viper_mesma_run')
help_but = widget_button(run_base, value='Help', tooltip='', $
    xsize=120, ysize=30, yoffset=20, event_pro='viper_mesma_help', sensitive=0)
save_but = widget_button(run_base, value='Save control file', tooltip='', $
    xsize=120, ysize=30, yoffset=20, event_pro='viper_mesma_save', sensitive=1)
restore_but = widget_button(run_base, value='Restore control file', tooltip='', $
    xsize=120, ysize=30, yoffset=20, event_pro='viper_mesma_restore', sensitive=1)

; Draw the widget heirarchy to the screen.
widget_control, tlb, /realize

; Get the graphic and draw it
if !version.os_family eq 'Windows' then begin
    logofile=FILE_SEARCH(STRSPLIT (!PATH, PATH_SEP(/SEARCH_PATH), /EXTRACT) + '\viper_logo.sav')
endif else begin
    logofile='/usr/local/rsi/idl_6.1/products/envi_4.1/save_add/viper_logo.sav'
    ;logofile=FILEPATH('viper_logo.sav', subdirectory=['save_add'])
endelse
if file_test(logofile[0]) eq 1 then begin
    restore, logofile[0]
    ; Obtain the window index.
    WIDGET_CONTROL, graphic, GET_VALUE = win1
    ; Set the new widget to be the current graphics window
    WSET, win1
    tv, icon, true=1
endif else begin
    result=dialog_message('viper_logo.sav not found')
    ;return
endelse

; Make a structure to store info about the state of the widget program.
state = { $
    infile:ptr_new(), $
    inpath:'',$
    outfile:ptr_new(),$
    nlibs:0, $
    scale:0, $
    scale_drop:scale_drop, $
    specnames1:ptr_new(), $
    specnames2:ptr_new(), $
    specnames3:ptr_new(), $
    specnames5:ptr_new(), $
    selected1:ptr_new(), $
    selected2:ptr_new(), $
    selected3:ptr_new(), $
    selected5:-1, $
    total_txt:total_txt,$
    check1:check1,$
    check2:check2,$
    check3:check3,$
    check4:check4,$
    check5:check5,$
    checks:bytarr(5),$
    min_slider:min_slider, $
    max_slider:max_slider, $
    shade_slider:shade_slider, $
    residv_text:residv_text, $
    maxrmse_text:maxrmse_text, $
    residv_lab:residv_lab, $
    maxrmse_lab:maxrmse_lab, $
    ;maxrmse_slider:maxrmse_slider, $
    ;residv_slider:residv_slider, $
    residb_slider:residb_slider,$
    file0text:file0text, $
    file0browse:file0browse,$
    fid0:-1,$
    file1text:file1text, $
    file1fulltext:'', $
    file1browse:file1browse,$
    fid1:-1,$
    file1subset:file1subset, $
    file2text:file2text, $
    file2fulltext:'', $
    file2browse:file2browse,$
    file2subset:file2subset, $
    fid2:-1,$
    file3text:file3text, $
    file3fulltext:'', $
    file3browse:file3browse,$
    file3subset:file3subset, $
    fid3:-1,$
    sel5lab:sel5lab, $
    file5text:file5text, $
    file5fulltext:'',$
    file5browse:file5browse,$
    file5subset:file5subset, $
    fid5:-1,$
    sel4lab:sel4lab, $
    file4text:file4text, $
    file4browse:file4browse,$
    photo_shade_but:photo_shade_but, $
    nonphoto_shade_but:nonphoto_shade_but, $
    photo_shade:1,$
    p1_slide:p1_slide,$
    p2_slide:p2_slide,$
    p3_slide:p3_slide,$
    constraints:[-.05,1.05,0.80,0.025,0.025,7],$
    out_frac:out_frac,$
    out_class:out_class,$
    out_resid:out_resid,$
    outputs:[1,1,0], $
    open_files_but:open_files_but, $
    open_files:1, $
    multiple:0}

; Make a pointer for the state structure & store in the top-level base's user value.
pstate = ptr_new(state, /no_copy)
widget_control, tlb, set_uvalue=pstate

; Call XMANAGER to register the widget program & start event handling. Do not use cleanup
; routine, instead will get info from pstate, free pointer and return structure below
xmanager, 'viper_mesma', tlb, /no_block, cleanup='viper_mesma_cleanup'

end