;+ ; 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 I/O ; 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 GUI ; 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 |