;==========================================================================
;PCA TOOL
;Written by David Waddill, April-May 1993
;Modified by Paul Johnson, 1994
;Modified by Paul Johnson for Win32, 1996
;Modified by David Klassen, for vector file list input, command line
;	arguements, presentation GIF output images: 1996
;To compile: .run pcatool
;To run:  pcatool[, filevect][, /fullimage][, /skipzap][, /norm]
;	[, datafile=input.data][, ingStoN=signal-to-noise-ratio]
;	[, incl=included-pixels-list|excl=excluded-pixels-list]
;	[, exclimgf=excluded-pixels-fits-image]
;==========================================================================
; Compiling the rest of the program
@pcalib3
@pcalib1
@pcalib2
@pcalib4
@strno
@gifbar
;
;==========================================================================
PRO pca_init
;
;This procedure runs the images stored in PCArr through PCAINIT.  Windows
;are opened in order to display PCA plots of 1st vs. 0th, 2nd vs. 0th
;and 1st vs. 2nd.
;
;     Image   =  array containing displayed image in window 3 (original).
;     Im_Num  =  image number currently displayed in window 3.
;     Elem1   =  1st eigenvector displayed on PCA Plot 0.
;     Elem2   =  2nd eigenvector displayed on PCA Plot 0.
;     Elem3   =  3rd eigenvector displayed on PCA Plot 1.
;     Elem4   =  4th eigenvector displayed on PCA Plot 1.
;     Elem5   =  5th eigenvector displayed on PCA Plot 2.
;     Elem6   =  6th eigenvector displayed on PCA Plot 2.
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
             X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
;Forming data into an array # of pixels by number of images in size.
NElem = Long(XSize) * Long(YSize)
PCArr = fltarr(NElem,NImages)
temp = fltarr(NElem)
for i = 0,NImages-1 do begin
	temp = reform(Data1(i,XLow:XHigh,YLow:YHigh),NElem)
	PCArr(*,i) = temp
endfor ;i
;
;All the data is now stored in the array PCArr.  Now going to run the data
;through the PCA routine PCAINIT.
pcainit,PCArr,object=good_pix,tek=0 ;in pcalib3.pro
;
N_iter = 1
Elem1 = 0
Elem2 = 1
Elem3 = 0
Elem4 = 2
Elem5 = 2
Elem6 = 1
;xrange=[[min(new_vecs(*,Elem1)),max(new_vecs(*,Elem1))], $
;	[min(new_vecs(*,Elem3)),max(new_vecs(*,Elem3))], $
;	[min(new_vecs(*,Elem5)),max(new_vecs(*,Elem5))]]
;yrange=[[min(new_vecs(*,Elem2)),max(new_vecs(*,Elem2))], $
;	[min(new_vecs(*,Elem4)),max(new_vecs(*,Elem4))], $
;	[min(new_vecs(*,Elem6)),max(new_vecs(*,Elem6))]]
xrange=[[min(new_vecs(good_pix,Elem1)),max(new_vecs(good_pix,Elem1))], $
	[min(new_vecs(good_pix,Elem3)),max(new_vecs(good_pix,Elem3))], $
	[min(new_vecs(good_pix,Elem5)),max(new_vecs(good_pix,Elem5))]]
yrange=[[min(new_vecs(good_pix,Elem2)),max(new_vecs(good_pix,Elem2))], $
	[min(new_vecs(good_pix,Elem4)),max(new_vecs(good_pix,Elem4))], $
	[min(new_vecs(good_pix,Elem6)),max(new_vecs(good_pix,Elem6))]]
;
;Opening up window in order to display PCAPLOT.
window,2,xs=400,ys=400,title='PCA Plot 2',retain=2
window,1,xs=400,ys=400,title='PCA Plot 1',retain=2
window,0,xs=400,ys=400,title='PCA Plot 0',retain=2
;window,2,xs=600,ys=400,title='PCA Plot 2',retain=2
;window,1,xs=600,ys=400,title='PCA Plot 1',retain=2
;window,0,xs=600,ys=400,title='PCA Plot 0',retain=2
;
;Finding pca coordinate for the sky pixel.
;Sky_pca = newcoo(Sky)
;
;Displaying PCA plots
wset,2
;pcaplot,Elem5,Elem6
pcaplot,Elem5,Elem6,0,ps=1,iter=N_iter,outfile=2,$
	xrange=xrange(*,2),yrange=yrange(*,2) ; in pcalib3.pro
w2x = !x
w2y = !y
;!x = w2x
;!y = w2y
;oplot,[Sky_pca(Elem5)],[Sky_pca(Elem6)],psym=1,color=1
wset,1
;pcaplot,Elem3,Elem4
pcaplot,Elem3,Elem4,1,ps=1,iter=N_iter,outfile=1,$
	xrange=xrange(*,1),yrange=yrange(*,1) ;in pcalib3.pro
w1x = !x
w1y = !y
;!x = w1x
;!y = w1y
;oplot,[Sky_pca(Elem3)],[Sky_pca(Elem4)],psym=1,color=1
wset,0
;pcaplot,Elem1,Elem2 ;in pcalib3.pro
pcaplot,Elem1,Elem2,2,ps=1,iter=N_iter,outfile=0,$
	xrange=xrange(*,0),yrange=yrange(*,0) ;in pcalib3.pro
w0x = !x
w0y = !y
;!x = w0x
;!y = w0y
;oplot,[Sky_pca(Elem1)],[Sky_pca(Elem2)],psym=1,color=1
pcaplot,Elem1,Elem2,2,ps=1,iter=N_iter,outfile=3,$
	xrange=xrange(*,0),yrange=yrange(*,0),/thrd ;in pcalib3.pro
clean = intarr(4)
clean(*) = 0
;
return
end

;==========================================================================
PRO specextract,spec,specerr,color
;
; A procedure to extract spectra from the data set
;
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
              X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON COLOR, c_pix, Num_pix
;
XXspec = c_pix1 mod InX
YYspec = c_pix1 / InX
tempspec = fltarr(Num_pix,WB)
;
for ii=0l,Num_pix-1 do begin
	tempspec(ii,*) = Data1(*,XXspec(ii),YYspec(ii))
;	norm=1000.*WB/total(tempspec(ii,*))
;	tempspec(ii,*)=tempspec(ii,*)*norm
endfor ;ii
;
for i = 0,WB-1 do begin
	spec(color-1,i) = Total(tempspec(*,i))/Num_pix
endfor ;i
;
for i = 0,WB-1 do begin
	specerr(color-1,i) = 0.
	for ii=0l,Num_pix-1 do begin
		specerr(color-1,i) = specerr(color-1,i)+ $
		   (tempspec(ii,i)-spec(color-1,i))^2
	endfor ;ii
	specerr(color-1,i) = sqrt(specerr(color-1,i))/Num_pix
endfor ;i
return
end

;==========================================================================
PRO specwin,spec,specerr,color
;
; A procedure to plot up extracted spectra in color
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON PCA1, wavelength, eigvecs, wave_title
nspec = color-2
tempx = !X.S
tempy = !Y.S
window,5,xs=400,ys=400,title='PCA pixel spectra',retain=2
;window,5,xs=600,ys=400,title='PCA pixel spectra',retain=2
!Y.TITLE='Image Pixel Value'
!X.TITLE=wave_title
wset,5
;band = indgen(WB)
plotmin = min(spec(0:nspec,*))
plotmax = max(spec(0:nspec,*))
plot,wavelength,[plotmin,plotmax],/nodata,yrange=[plotmin,plotmax]
;tempspec = spec(0,*)
for i = 0,nspec do begin
;	spec(i,*) = spec(i,*)*1000./tempspec
	oplot,wavelength,spec(i,*),color = i+1
	errplot,wavelength,spec(i,*)-specerr(i,*),spec(i,*)+specerr(i,*)
endfor ;i
openw,6,'spec.tmp',width=300 ;Write out the spectra to a text file
;for ii = 0,nspec do begin
	printf,6,'### ',nspec+1,WB
;	for iii = 0,WB-1 do begin
;		printf,6,wavelength(iii),spec(ii,iii),specerr(ii,iii)
;		printf,6,wavelength(iii),spec(0:nspec,iii),specerr(0:nspec,iii)
;	endfor ;iii
;endfor ;ii
	printf,6,[transpose(wavelength),(spec[0:nspec,*]),$
		(specerr[0:nspec,*])]
close,6
!X.S = tempx
!Y.S = tempy
return
end

;==========================================================================
PRO analyze
;
;Main menu calling procedure for the PCA analysis tools.
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
             X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
;Top-level menu
Choice = 0
while (Choice NE 1) do begin
	menu = 0
	xmenu, ['RETURN to Main Menu',$
		'Plot EIGENVALUES',$
		'Plot EIGENVECTOR Spectra',$
		'Display Eigen-image',$
		'SELECT Endmembers',$
		'MEASURE Fractional Abundances',$
		'Plot ENDMEMBER Spectra',$
		'CHANGE Component Display Cutoffs',$
		'Display RESIDUAL Spectral Image',$
		'Display ALL RESIDUAL Spectral Images',$
	   	'Change COLOR PALETTE',$
		'Analyze Pixel SUBSET Chosen From Fracmaps',$
		'Analyze RESIDUAL Data',$
		'NEGATE an eigenvector'], BASE=base, BUTTONS=B, $
	title='Analysis Tools Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = menu, event.id
	WIDGET_CONTROL, event.top, /DESTROY
	case menu of
	; Return to Main Menu
	0: begin
		Choice = 1
		refr_pca,/color ;in pcalib1.pro
	end ;0
	;Eigenvalue Plot
	1: begin
		lampltsetup ;in pcalib4
		print,'*****       Eigenvalue Plot Complete        *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;1
	;Eigenvector Spectral Plots
	2: begin
		eigpltsetup ;in pcalib4
		print,'*****       Eigenvector Spectral Plots Complete        *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;2
	;Eigen-image Display
	3: begin
		pamtvsetup ;in pcalib4
		print,'*****           Eigen-image Display Complete           *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;3
	; Select Endmembers
	4: begin
		refr_pca,/color ;in pcalib1.pro
		s_endm ;in pcalib2.pro
		print,'*****          Endmembers Selected           *****'
		print,'***** Ready to Measure Fractional Abundances *****'
	end ;4
	; Measure Fractional Abundances
	5: begin
		frac_ab ;in pcalib2.pro
		d_frac ;in pcalib2.pro
		d_resid ;in pcalib2.pro
		print,'*****  Fractional Abundances and Residuals Calculated  *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;5
	;Endmember Spectral Plots
	6: begin
		endspcsetup ;in pcalib4.pro
		print,'*****         Endmember Spectral Plots Complete        *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;6
;	;Particular Pixel Spectral Plots
;	5: begin
;		spcpltsetup ;in pcalib4.pro
;	end ;5
	; Change Component Display Cutoffs
	7: begin
		rscl_comp ;in pcalib1.pro
		print,'*****              Image Contrast Adjusted             *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;7
	; Display Residual Spectral Image
	8: begin
		dispresidsetup ;in pcalib4.pro
		print,'*****         Residual Image Display Complete          *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;8
	; Display All Residual Spectral Images
	9: begin
		mrestv,reswin=21,/verbose ;in pcalib4.pro
		print,'*****          Residual Images Display Complete        *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;9
	;Change color palette
	10: begin
		read,colorvalue,prompt='Enter color palette number (0-41): '
		if colorvalue le 40 then loadct,colorvalue $
		else if colorvalue eq 41 then lct_dmw $
		else loadct,0
		print,'*****               Color Palette Changed              *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;10
	; Analyze Pixel Subset Chosen from Fracmaps
	11: begin
		an_subset ;in pcalib4.pro
	end ;11
	; Analyze Residual Data
	12: begin
		print,'*****    Analyzing Residuals     *****'
		print,'***** Ready to Select Endmembers *****'
		an_resid ;in pcatool.pro
	end ;12
	; Negate an eigenvector and eigen image
	13: begin
		read,negnum,prompt='Enter the number of the eigenvector to negate: '
		eigvecs[*,negnum]=-1.0*eigvecs[*,negnum]
		new_vecs[*,negnum]=-1.0*new_vecs[*,negnum]
		xrange=[[min(new_vecs(good_pix,Elem1)),max(new_vecs(good_pix,Elem1))], $
			[min(new_vecs(good_pix,Elem3)),max(new_vecs(good_pix,Elem3))], $
			[min(new_vecs(good_pix,Elem5)),max(new_vecs(good_pix,Elem5))]]
		yrange=[[min(new_vecs(good_pix,Elem2)),max(new_vecs(good_pix,Elem2))], $
			[min(new_vecs(good_pix,Elem4)),max(new_vecs(good_pix,Elem4))], $
			[min(new_vecs(good_pix,Elem6)),max(new_vecs(good_pix,Elem6))]]
		clean[*]=1
		refr_pca,/color,/ps ;in pcalib1.pro
		print,'*****         Eigenvector Negation Done                *****'
		print,'***** Ready for Display Functions or Residual Analysis *****'
	end ;13
	endcase
endwhile
return
end

;==========================================================================
PRO an_resid
;
;Procedure for analyzing the residual information found in d_resid.
;Modified to remove extra PCAplots

COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
             X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
n_r_pix = n_elements(res_pix)
print, 'Number of Residual pixels',n_r_pix
printf,5,'Number of Residual pixels',n_r_pix
;
;Saving old PCA coordinates.
;old_vecs = new_vecs
;
;Finding subset of good_pix and res_pix for PCAINIT.
good_pix2 = intarr(XSize*YSize)
good_pix2 = long(good_pix2)
good_pix3 = intarr(XSize*YSize)
good_pix3 = long(good_pix3)
good_pix4 = intarr(XSize*YSize)
good_pix4 = long(good_pix4)
good_pix2(res_pix) = 1
good_pix3(good_pix) = 1
good_pix4 = good_pix3 + good_pix2
good_pix = where(good_pix4 eq 2)
print,'# good elements = ', n_elements(good_pix)
;
;Running high residual pixels through PCA again.
pcainit,PCArr,object=good_pix ,tek=0 ;in pcalib3.pro
;
;Finding pca coordinate for the sky pixel.
;Sky_pca2 = newcoo(Sky)
;
;Finding pca coordinates for old end-members.
EM_npca = fltarr(n_oem_pix,NImages)
for i = 0, n_oem_pix-1 do begin
  EM_npca(i,*) = newcoo(PCArr(old_em_pix(i),*)) ;in pcalib3.pro
endfor ;i
;
;Opening up windows in order to display new PCAPLOTs.
window,2,xs=400,ys=400,title='New PCA Plot 2',retain=2
window,1,xs=400,ys=400,title='New PCA Plot 1',retain=2
window,0,xs=400,ys=400,title='New PCA Plot 0',retain=2
;window,2,xs=600,ys=400,title='New PCA Plot 2',retain=2
;window,1,xs=600,ys=400,title='New PCA Plot 1',retain=2
;window,0,xs=600,ys=400,title='New PCA Plot 0',retain=2
;
;Displaying PCAPLOTs.
Elem1 = 0
Elem2 = 1
Elem3 = 0
Elem4 = 2
Elem5 = 2
Elem6 = 1
xrange=[[min(new_vecs(*,Elem1)),max(new_vecs(*,Elem1))], $
	[min(new_vecs(*,Elem3)),max(new_vecs(*,Elem3))], $
	[min(new_vecs(*,Elem5)),max(new_vecs(*,Elem5))]]
yrange=[[min(new_vecs(*,Elem2)),max(new_vecs(*,Elem2))], $
	[min(new_vecs(*,Elem4)),max(new_vecs(*,Elem4))], $
	[min(new_vecs(*,Elem6)),max(new_vecs(*,Elem6))]]
wset,2
pcaplot,Elem5,Elem6,0,ps=1,iter=N_iter+1,outfile=2,$
	xrange=xrange(*,2),yrange=yrange(*,2) ;in pcalib3.pro
;oplot,[Sky_pca2(Elem5)],[Sky_pca2(Elem6)],psym=1,color=1
for jj = 0, n_oem_pix-1 do $
  oplot,[EM_npca(jj,Elem5)],[EM_npca(jj,Elem6)],psym=1,color=4
w2x = !x
w2y = !y
wset,1
pcaplot,Elem3,Elem4,1,ps=1,iter=N_iter+1,outfile=1,$
	xrange=xrange(*,1),yrange=yrange(*,1) ;in pcalib3.pro
;oplot,[Sky_pca2(Elem3)],[Sky_pca2(Elem4)],psym=1,color=1
for jj = 0, n_oem_pix-1 do $
  oplot,[EM_npca(jj,Elem3)],[EM_npca(jj,Elem4)],psym=1,color=4
w1x = !x
w1y = !y
wset,0
pcaplot,Elem1,Elem2,2,ps=1,iter=N_iter+1,outfile=0,$
	xrange=xrange(*,0),yrange=yrange(*,0) ;in pcalib3.pro
;oplot,[Sky_pca2(Elem1)],[Sky_pca2(Elem2)],psym=1,color=1
for jj = 0, n_oem_pix-1 do $
	oplot,[EM_npca(jj,Elem1)],[EM_npca(jj,Elem2)],psym=1,color=4
w0x = !x
w0y = !y
pcaplot,Elem1,Elem2,2,ps=1,iter=N_iter+1,outfile=3,$
	xrange=xrange(*,0),yrange=yrange(*,0),/thrd ;in pcalib3.pro
clean = intarr(4)
clean(*) = 0
;
;Select new set of end_members for residual pixels
Choose = 0
while (Choose ne 1) do begin
	xmenu, ['DISPLAY PCA pixels','New PCA Plot','ZOOM in on PCA Plot', $
		'REFRESH Image and PCA plots',$
		'SELECT Endmembers', '     MEASURE Fractional Abundances     ',$
		'ITERATE Residuals Again','Change COLOR PALETTE', 'STOP'], $
		BASE=base, BUTTONS=B, title='Residual Analysis Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = menu, event.id
	WIDGET_CONTROL,event.top,/DESTROY
;
	case menu of
	;Select pca pixels to be displayed
	0: begin
		Ans = 'y'
		color=1
		refr_pca ;in pcalib1.pro
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		c = reform(b,Long(XSize)*Long(YSize))
		c_pix2 = bytscl(c) > 8b
		spec = fltarr(10,WB)
		specerr = fltarr(10,WB)
		window,4,xs=XSize*scl,ys=YSize*scl,title='PCA pixels',retain=2,$
			yp=0,xp=300
		while (Ans ne 'n') do begin
			s_pcaplot ;in pcalib1.pro
			get_box,color ;in pcalib2.pro
			;Display image with pixels highlighted in color from color_im
			pca_pix ;in pcalib1.pro
			if (Num_pix le 0) then print, 'No points inside boxed area.'
			;Display selected points, in color, in window 4.
			if (Num_pix gt 0) then begin
				c_pix2(c_pix1) = color
				c_pix3 = reform(c_pix2,XSize,YSize)
				c_pix4 = rebin(c_pix3,XSize*scl,YSize*scl,/sa)
				specextract,spec,specerr,color ;in pcatool.pro
				wset,4
				tv, c_pix4
				c_pix=c_pix1
				d_c_pca,color
				color = color+1
			endif
			xmenu,['    No    ','    Yes    '],base=base, $
				buttons=b,title='Again?'
			WIDGET_CONTROL,base,/REALIZE
			event = WIDGET_EVENT(base)
			WIDGET_CONTROL, GET_UVALUE = dispmore, event.id
			WIDGET_CONTROL, event.top, /DESTROY
			case dispmore of
			0: Ans='n'
			1: Ans='y'
			endcase
		endwhile
		specwin,spec,specerr,color ;in pcatool.pro
		clean(*)=1
		print,'*****   Display Routine Complete    *****'
		print,'***** Ready for Endmember Selection *****'
	end ;0
	; New PCA Plot
	1: begin
		N_iter=N_iter+1
		d_pcaplot ;in pcalib1.pro
		N_iter=N_iter-1
		print,'*****    New PCA Plot Displayed     *****'
		print,'***** Ready for Endmember Selection *****'
	end ;1
	; Zoom in on PCA Plot
	2: begin
		ans=''
		for r=0,2 do begin
			print,'Change x-range for plot'+strcompress(string(r))+'? '
			print,'(Current range ='+strcompress(string(xrange(0,r)))+$
				','+strcompress(string(xrange(1,r)))+')'
			read,ans
			if strmid(ans,0,1) eq 'y' then begin
				print,'Enter xlow,xhigh'
				read,x1,x2
				xrange(0,r)=x1
				xrange(1,r)=x2
				clean(r)=1
			endif
            print,'Change y-range for plot'+strcompress(string(r))+'? ' 
			print,'(Current range ='+strcompress(string(yrange(0,r)))+$
				','+strcompress(string(yrange(1,r)))+')'
            read,ans 
            if strmid(ans,0,1) eq 'y' then begin 
                print,'Enter ylow,yhigh' 
                read,y1,y2
				yrange(0,r)=y1
				yrange(1,r)=y2
				clean(r)=1
            endif
		endfor ;r
		refr_pca,/ps,/zoom,/iter1,/color ;in pcalib1.pro
		print,'*****    New PCA Plot Displayed     *****'
		print,'***** Ready for Endmember Selection *****'
	end ;2
	; Refreshing Image and PCA plots
	3: begin
		print,'Refreshing Image and PCA plots...'
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		print,'***** Image and PCA Plot Refreshed  *****'
		print,'***** Ready for Endmember Selection *****'
	end ;3
	;Choose endmembers
	4: begin
;		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		s_endm ;in pcalib2.pro
		print,'*****     Endmembers Selection Complete      *****'
		print,'***** Ready to Calcualte Frational Abuncaces *****'
	end ;4
	;Measure endmember abundances
	5: begin
		N_iter = N_iter+1
		frac_ab ;in pcalib2.pro
		d_frac ;in pcalib2.pro
		d_resid ;in pcalib2.pro
		print,'*****   Frational Abundances and Residuals Calculated    *****'
		print,'*****  For Endmember and Eigenvector Display Functions   *****'
		print,'***** Choose ITERATE RESIDUALS AGAIN and make selections *****'
		print,'*****                 From Analysis Menu                 *****'
	end ;5
	;Repeat iteration
	6: begin
		choose=1
		print,'=============================================='
		print,'CHOOSE ANALYZE RESIDUALS FROM MENU'
		print,'OR CHOOSE ONE OF THE DISPLAY FUNCTIONS'
		print,'=============================================='
		print
	end ;6
	;Change color palette
	7: begin
		read,colorvalue,prompt='Enter color palette number (0-41): '
		if colorvalue le 40 then loadct,colorvalue $
		else if colorvalue eq 41 then lct_dmw $
		else loadct,0
		print,'***** Color Palette Changed *****'
		print,'***** Continue with analysis *****'
	end ;7
	;Save and Stop
	8: begin
		loadct,0
		choose =1
	end ;8
	endcase
endwhile
;
for i = 1,N_EndM-1 do old_em_pix(i-1+n_oem_pix) = EM_pix(i)
n_oem_pix = n_oem_pix + N_EndM - 1
N_EM(N_iter) = N_EndM-1
;
return
end
;==========================================================================
PRO pcatool,filevect,fullimg=fullimg,skipzap=skipzap, $
	datafile=df,imgStoN=imgStoN,norm=norm,scale=scale,$
	incl=incl,excl=excl,eximgf=eximgf
;
;Main menu calling procedure for the PCA tool.
;
;	filevect = input variable containing an array of strings, each
;		element being the name of an image in FITS format
;	fullimg = switch to let program know you want to use the full
;		full image - skips the 'choose sub-image' routine
;	skipzap = switch to let program know you want to use all
;		the image pixels - skips 'zap bad pixels' routine
;	datafile = string input variable containing the name of the
;		input data file - will be asked for if this is not set
;	imgStoN = image signal-to-noise level; will set to a defaul value
;		of 20.0 if not set on command line
;	norm = keyword to normalize each input image by its mean
;	scale = set for display scale to override the autoscaling
;	incl = a vector list of pixels to use in PCA; if set, the
;		skipzap switch will automatically be set
;	excl = a vector list of pixels exclude from  PCA; if set, the
;		skipzap switch will automatically be set
;		NOTE: do not set *both* include and exclude; if you do,
;		only the include will be used!
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
			X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON VECT, vector
COMMON RES2, v_names, StoN
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
if n_params() lt 1 then begin
	print,'USAGE: pcatool[,filevect][,/FULLIMG][,/SKIPZAP]'
	print,'     [,DATAFILE=data_input_file][,IMGSTON=imgStoN][,/NORM]'
	print,'     [,SCALE=scale][,INCL=incl][,EXCL=excl][,EXIMGF=eximgf]'
	print,''
	print,'filevect = input variable containing an array of strings, each'
	print,'           element being the name of an image in FITS format'
	print,'/fullimg = switch to let program know you want to use the full'
	print,'          image - skips the CHOOSE SUB-IMAGE routine'
	print,'/skipzap = switch to let program know you want to use all'
	print,'          the image pixels - skips ZAP BAD PIXELS routine'
	print,'datafile = string input variable containing the name of the'
	print,'           input data file - will be asked for if this is not set'
	print,'imgStoN = image signal-to-noise level; will set to a defaul value'
	print,'          of 20.0 if not set on command line'
	print,'/norm = keyword to normalize each input image by its mean'
	print,'scale = set for display scale to override the autoscaling of display'
	print,'incl = a vector list of pixels to use in PCA'
	print,'       if set, /SKIPZAP will automatically be set'
	print,'excl = a vector list of pixels exclude from  PCA'
	print,'       if set, /SKIPZAP will automatically be set'
	print,'   NOTE: do not set *both* include and exclude; if you do,'
	print,'   only the include will be used!'
   print,'eximgf = the filename of an image mask of pixels to exclude'
   print,'         The image is an INTARR with a value of 1 for pixels'
   print,'         to be excluded and zero for those to be included'
   print,'         This should be used in place of either incl or excl'
endif

if n_elements(filevect) ne 0 then vector=filevect
;Top-level menu
Choice = 1
while (Choice NE 0) do begin
	xmenu,['EXIT PCA Tool','Initialize','RECHOOSE Image', $
		'NEW PCAPLOT','ZOOM in on a PCA Plot', $
		'DISPLAY New Frame','CHANGE Image Display Cutoffs', $
		'DISPLAY PCA Pixels on Image','DISPLAY IMAGE Pixels on PCA Plot', $
		'Change COLOR PALETTE','REFRESH Image and PCA plots',$
		'ANALYSIS Tools Menu' ], BASE = base, BUTTONS=B, $
	   title='PCATool Main Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = Choice, event.id
	WIDGET_CONTROL, event.top, /DESTROY
; Must initialize first - set choice to 1 if no images and not exiting
	if (n_elements(NImages) eq 0)and (Choice ne 0) then Choice=1
	case Choice of
	; Exit PCATool
	0: begin
		loadct,0
		Choice=0
		close,5
		close,4
		!p.charsize=1
	end ;0
	; Initialize
	1: begin
		loadct,0
		print
		print, '=================================================='
		print, '             BEGINNING INITIALIZATION'
		print, '=================================================='
		print, 'Just one moment . . .'
		close,5
		close,4
;Read in the data and set up some files
		if not keyword_set(incl) then include = -1
		if n_elements(eximgf) eq 1 then begin
			exclimg=readfits(eximgf)
			exclude=where(exclimg ne 0)
		endif else $
		if not keyword_set(excl) then exclude = -1
		in_data,datafile=df,norm=keyword_set(norm) ;in pcalib2.pro
;Set Signal-to-Noise stuff; stop-gap for something more sophisticated
		if not(keyword_set(imgStoN)) then begin
			StoN=20.0
			print,'No Signal-to-Noise input so setting S/N = 20'
		endif else begin
			StoN=float(imgStoN)
			print,'Setting S/N = ',StoN
		endelse
;Initialize arrays and variables for future use
		clean = intarr(4)
		clean(*) = 0
		xrange=fltarr(2,3)
		yrange=fltarr(2,3)
		openw,5,'pca.log'
		openw,4,'spectra.dat'
		N_EM=intarr(10)
		old_em_pix=intarr(30) & old_em_pix=long(old_em_pix)
		n_oem_pix=0
		XSize = InX
		YSize = InY
		XLow = 0
		XHigh = InX-1
		YLow = 0
		YHigh = InY-1
		res_pix=make_array(XSize*YSize,/LONG,/INDEX)
;Display the image with the brightest median - presumably continuum
		if keyword_set(scale) then begin
			scl=scale
			print,'Image scale factor = ',scl
		endif  else det_scl ;in pcalib1.pro
		find_med ;in pcalib1.pro
		Im_Num=where(High eq max(High))
			Im_Num=Im_Num(0)
		print, 'Image number for display = ',Im_Num
		Hi=High(Im_Num) & Lo=Low(Im_Num)
		print,'Display scale = ',Lo,Hi
		window,3,xs=InX*scl,ys=InY*scl,title='Initial Image',retain=2
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		d = reform(b,InX,InY)
		Image = rebin(d,InX*scl,InY*scl,/sa)
		tv, bytscl(Image,max=Hi,min=Lo)
;Checking to use the full image or not
		if not(keyword_set(fullimg)) then $
			rsz_query ;in pcalib1.pro
;Checking to zap pixels
		if n_elements(exclude) eq 1 or exclude(0) ne -1 or include(0) ne -1 $
			then skipzap=1
		print
		print, 'Displaying selected sub-image for bad pixel extraction.'
		bad_pix1,skipzap=skipzap ;in pcalib1.pro
;Do the PCA
		print
		print, 'Running data through PCAINIT.'
		pca_init ;in pcatool.pro
;Displaying final image
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		d = reform(b,XSize,YSize)
		Image = rebin(d,XSize*scl,YSize*scl,/sa)
		print, '  Displaying Image ',Im_Num
		print, '  Wavelength = ',wavelength(Im_Num)
		window,3,xs=scl*XSize,ys=scl*YSize, $
			title='Selected Image being used in PCA',retain=2
		print, '  Scale Range = ',Lo,Hi
		tv, bytscl(Image,max=Hi,min=Lo)
;		gifbar,Image,Low(Im_Num),High(Im_Num),ct=0,tc=255,fname='original.gif',$
		gifbar,Image,Low(Im_Num),High(Im_Num),ct=0,tc=255,fname='original',$
			dispwin=20,wintitle='Original - Image '+strno(fix(Im_Num))
		BufImage=rebin(d,Xsize*xygrid,Ysize*xygrid)
		Outimage=fltarr(XSize,YSize)
		Outimage(xmargin:Xsize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
		writefits,'original.fit',Outimage
		print,'***** Initialization Complete *****'
		print,'*****   Ready for Analysis    *****'
	end ;1
	; Rechoose image
	2: begin
		window,3,xs=InX*scl,ys=InY*scl,title='Initial Image',retain=2
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		d = reform(b,InX,InY)
		c = rebin(d,InX*scl,InY*scl,/sa)
		tvscl, c > Low(Im_Num) < High(Im_Num)
		Image = c
		rsz_query ;in pcalib1.pro
		print
		print, 'Displaying selected sub-image for bad pixel extraction.'
		bad_pix1 ;in pcalib1.pro
		print
		print, 'Running data through PCAINIT.'
		pca_init ;in pcatool.pro
		print,'***** Initialization Complete *****'
		print,'*****   Ready for Analysis    *****'
	end ;2
	; New PCA Plot
	3: begin
		d_pcaplot ;in pcalib1.pro
		print,'***** New PCA Plot Displalyed *****'
		print,'*****   Ready for Analysis    *****'
	end ;3
	; Zoom in on PCA Plot
	4: begin
		ans=''
		for r=0,2 do begin
			print,'Change x-range for plot'+strcompress(string(r))+'? '
			print,'(Current range ='+strcompress(string(xrange(0,r)))+$
				','+strcompress(string(xrange(1,r)))+')'
			read,ans
			if strmid(ans,0,1) eq 'y' then begin
				print,'Enter xlow,xhigh'
				read,x1,x2
				xrange(0,r)=x1
				xrange(1,r)=x2
				clean(r)=1
			endif
            print,'Change y-range for plot'+strcompress(string(r))+'? ' 
			print,'(Current range ='+strcompress(string(yrange(0,r)))+$
				','+strcompress(string(yrange(1,r)))+')'
            read,ans 
            if strmid(ans,0,1) eq 'y' then begin 
                print,'Enter ylow,yhigh' 
                read,y1,y2
				yrange(0,r)=y1
				yrange(1,r)=y2
				clean(r)=1
            endif
		endfor ;r
		refr_pca,/ps,/zoom,/color ;in pcalib1.pro
		print,'***** New PCA Plot Displalyed *****'
		print,'*****   Ready for Analysis    *****'
	end ;4
	; Display New Frame
	5: begin
		print
		XD = [XLow,XHigh]
		YD = [YLow,YHigh]
		d_image, XD, YD, scl ;in pcalib1.pro
		print,'***** New Image Frame Displayed *****'
		print,'*****    Ready for Analysis     *****'
	end ;5
	; Change Image Display Cutoffs
	6: begin
		print
		print, '=================================================='
		print, '            CHANGING IMAGE CUTOFFS'
		print, '=================================================='
		print, 'Median =', MEDIAN(Data1(Im_Num,XLow:XHigh,YLow:YHigh))
		print, 'Minimum =', MIN(Data1(Im_Num,XLow:XHigh,YLow:YHigh))
		print, 'Low Cutoff =', Low(Im_Num)
		print, 'High Cutoff =', High(Im_Num)
		print, '===================================================='
		print, 'NEW Cutoffs (Low, High)? '
		read, Lowd, Highd
		High(Im_Num) = Highd
		Low(Im_Num) = Lowd
		Hi=High(Im_Num)
		Lo=Low(Im_Num)
		wset, 3
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		d = reform(b,XSize,YSize)
		Image = rebin(d,XSize*scl,YSize*scl,/sa)
		erase
		tvscl, Image < Hi > Lo
		gifbar,Image,Lowd,Highd,ct=0,tc=255,fname='original', $
			dispwin=20,wintitle='Original - Image '+strno(fix(Im_Num))
		BufImage=rebin(d,Xsize*xygrid,YSize*xygrid)
		OutImage=fltarr(XSize,YSize)
		OutImage(xmargin:XSize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
		writefits,'original.fit',OutImage
		clean(3)=0
		print,'***** Image Contrast Adjusted *****'
		print,'*****   Ready for Analysis    *****'
	end ;6
	; Display PCA Pixels on Image
	7: begin
		Ans = 'y'
		clean(*)=1
		refr_pca ;in pcalib1.pro
		color=1
		kcolor=intarr(n_elements(Data1)) & kcolor(*)=255
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		c = reform(b,Long(XSize)*Long(YSize))
;		c_pix2 = bytscl(c) > 8b
		c_pix2 = bytscl(c,max=High(Im_Num),min=Low(Im_Num))
		window,4,xs=XSize*scl,ys=YSize*scl,title='PCA pixels',retain=2,$
			yp=300,xp=400
		spec=fltarr(10,WB)
		specerr =fltarr(10,WB)
		while (Ans ne 'n') do begin
			s_pcaplot ;in pcalib1.pro
			get_box,color ;in pcalib2.pro
			;Display image with pixels highlighted in color from  color_im
			pca_pix ;in pcalib1.pro
			if (Num_pix le 0) then print, 'No points inside boxed area.'
			;Display selected points, in color, in window 4.
			if (Num_pix gt 0) then begin
				bxcen = c_pix1[Num_pix/2]
				print,'    Image coordinates of box center = ',$
					bxcen mod XSize, bxcen/XSize
				c_pix2(c_pix1) = color
				c_pix3 = reform(c_pix2,XSize,YSize)
				c_pix4 = rebin(c_pix3,XSize*scl,YSize*scl,/sa)
				specextract,spec,specerr,color ;in pcatool.pro
				wset,4
				tv, c_pix4
				c_pix=c_pix1
				d_c_pca,color ;in pcalib2.pro
				kcolor(c_pix1)=color
				color = color+1
			endif
			xmenu,['    No    ','    Yes    '],base=base, $
				buttons=b,title='Again?'
			WIDGET_CONTROL,base,/REALIZE
			event = WIDGET_EVENT(base)
			WIDGET_CONTROL, GET_UVALUE = dispmore, event.id
			WIDGET_CONTROL, event.top, /DESTROY
			case dispmore of
			0: Ans='n'
			1: Ans='y'
			endcase
		endwhile
		specwin,spec,specerr,color ;in pcatool.pro
		wset,0
		pcaplot,Elem1,Elem2,kcolor,xrange=xrange(*,0),yrange=yrange(*,0),$
			zoom=zoom,/ps,outfile='0-im',iter=N_iter
		wset,1
		pcaplot,Elem3,Elem4,kcolor,xrange=xrange(*,1),yrange=yrange(*,1),$
			zoom=zoom,/ps,outfile='1-im',iter=N_iter
		wset,2
		pcaplot,Elem5,Elem6,kcolor,xrange=xrange(*,2),yrange=yrange(*,2),$
			zoom=zoom,/ps,outfile='2-im',iter=N_iter
		clean(*) = 1
		wset,4
		tvlct,r,g,b,/get
		write_png,'pca-im_'+strcompress(string(N_iter), /remove_all)+'.png',$
			c_pix4,r,g,b
		set_plot,'ps',/copy
		device,filename='pca-im_'+strcompress(string(N_iter),/remove_all)+'.eps'
		device,/color
		device,bits_per_pixel=24
		device,/encapsul
		device,xsize=XSize*scl/29.,ysize=YSize*scl/29.
		tv,c_pix4
		device,/close_file
		x
		print,'***** Display Routine Finished *****'
		print,'*****    Ready for Analysis    *****'
	end ;7
	; Display Image Pixels on PCA Plot
	8: begin
		clean(*)=1
		refr_pca ;in pcalib1.pro
		color = 1
		Ans = 'y'
		win = 3
		while (Ans ne 'n') do begin
			get_box, color, 30 ;in pcalib2.pro
			print,'Box size=30 with center at',X_im,Y_im
			image_pix ;in pcalib1.pro
			print,'Number of pixels in box:',Num_pix
			d_c_pca,color ;in pcalib2.pro
			color = color + 1
			xmenu,['    No    ','    Yes    '],base=base, $
				buttons=b,title='Again?'
			WIDGET_CONTROL,base,/REALIZE
			event = WIDGET_EVENT(base)
			WIDGET_CONTROL, GET_UVALUE = dispmore, event.id
			WIDGET_CONTROL, event.top, /DESTROY
			case dispmore of
			0: Ans='n'
			1: Ans='y'
			endcase
		endwhile
		wset,3
		tvscl, Image < High(Im_Num) > Low(Im_Num)
		clean(*)=1
		print,'***** Display Routine Finished *****'
		print,'*****    Ready for Analysis    *****'
	end ;8
	;Change color palette
	9: begin
		read,colorvalue,prompt='Enter color palette number (0-41): '
		if colorvalue le 40 then loadct,colorvalue $
		else if colorvalue eq 41 then lct_dmw $
		else loadct,0
		print,'***** Color Palette Changed *****'
		print,'*****  Ready for Analysis   *****'
	end ;9
	; Refreshing Image and PCA plots
	10: begin
		ans=''
		read,ans,prompt='Refresh PostScript copies too [y/N]? '
		print,'Refreshing Image and PCA plots...'
		clean(*)=1
		if ans eq 'y' then refr_pca,/color,/ps else $
			refr_pca,/color ;in pcalib1.pro
		print,'***** Image and Plot  Display Refreshed *****'
		print,'*****        Ready for Analysis         *****'
	end ;10
	; Analysis Tools Menu
	11: begin
		n_oem_pix=0
		print,'*****     Opening Analysis Menu     *****'
		print,'***** Ready for Endmember Selection *****'
		analyze ;in pcatool.pro
	end ;11
	endcase
endwhile
;
END