;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Reads the first line of a meta file which should contain the number of years
;;	we will need to analyze
;;
;;	OUT-DATED, WE CAN NOW GET THIS INFO FROM THE ASSOCIATED ENVI HEADER FILES
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function readyears, unit
	 ny=0
	 readf, unit, ny
	 print, ny
	 return, ny
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Reads the name of the "spatial" century file, and additional data from
;;	the necessary associated header file.  
;;	    returns a struct :
;;		    ns:number of samples in the file
;;		    nl:number of lines in the file
;;		    nb:number of bands in the file
;;		    map: an envi map structure (includes utm coords and pixelsize)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getCenturyFileInfo, unit
	 name=''
	 readf, unit, name

;	 envi_open_file, name, r_fid=id
;	 
;	 if id eq -1 then return, {errorstruct, name:name, ns:-1, nl:-1, nb:-1}
;
;	 envi_file_query, id, nb=nb, nl=nl, ns=ns, h_map=maph
;	 handle_value, maph, map
;	 envi_file_mng, id=id, /remove
;	 print, name
;	 return, {centstruct, name:name, ns:ns, nl:nl, nb:nb, map:map}
         return, getFileInfo(name)
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Reads the name of the landsat ndvi file and returns a structure identical to 
;;	the structure returned by getCenturyFileInfo
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getLandsatInfo, unit
	 name=''
	 readf, unit, name

;	 envi_open_file, name, r_fid=id
;	 
;	 if id eq -1 then return, {errorstruct, name:name, ns:-1, nl:-1, nb:-1}
;
;	 envi_file_query, id, nb=nb, nl=nl, ns=ns, h_map=maph
;	 handle_value, maph, map
;	 envi_file_mng, id=id, /remove
;	 print, name
;	 return, {ndvistruct, name:name, ns:ns, nl:nl, nb:nb, map:map}
         return, getFileInfo(name)
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Finds the offsets in number of pixels for the NDVI landsat images
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function calcOffset, ls, cent
	 xoff=(ls.map.mc(2) - cent.map.mc(2))/cent.map.ps(0)
	 yoff=(cent.map.mc(3) - ls.map.mc(3))/cent.map.ps(1)
	 ratio = cent.map.ps(0) / ls.map.ps(0)
	 
	 if xoff lt 0 or yoff lt 0 then return, {offsetstruct, x:-1, y:fix(yoff), factor:ratio}

	 return, {offsetstruct, x:fix(xoff), y:fix(yoff), factor:ratio}
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Main program.  
;;	Gets data from a meta file then reads century data and ndvi data
;;	gets map data from header files (ENVI MUST BE LOADED WHEN THIS IS RUN)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ndvivcent, meta
  
  openr, metaun, /get, meta
  
  readf, metaun, nclasses
  
;;Reads the file names from a single meta file and reads file info
;;actually these two procedures are identical, but are seperated for clarities sake?
;;	   mostly because at one point they weren't identical
  CentInfo=getCenturyFileInfo(metaun)
  LsInfo=getLandsatInfo(metaun)
;	CentInfo=getFileInfo(metaun)
;	LsInfo=getFileInfo(metaun)
  if LsInfo.nb eq -1 or CentInfo.nb eq -1 then begin
     print, 'ERROR OPENING NDVI or CENTURY FILES'
     print, 'NDVI Info = ', LsInfo
     print, 'Century Info = ', CentInfo
     close, metaun	&free_lun, metaun
     return, -1
  endif
  
  classes=' '
  readf, metaun, classes
  close, metaun	&free_lun, metaun
  
  
;;Calculates the offset of one file from the other
;;infact this program assumes the ndvi is a subset within the century, this is true
;;as long as the ndvi image was created with combunequal.pro, and the cent image used
;;one of those ndvi images for it's header info
  offset=calcOffset(LsInfo, CentInfo)
  if offset.x eq -1 then begin
     print, 'ERROR CALCULATEING OFFSETS'
     print, 'NDVI Info = ', LsInfo
     print, 'Century Info = ', CentInfo
     return, -1
  endif
  
  print, LsInfo
  print, '**************************************************'
  print, CentInfo
  print, '**************************************************'
  print, offset

  openr, class_un, /get, classes
  openr, ndvi_un, /get, LsInfo.name
  openr, cent_un, /get, CentInfo.name
  
  centline=fltarr(CentInfo.ns)
  ndviline=make_array(LsInfo.ns, type=LsInfo.type)
  classline=bytarr(LsInfo.ns)
  classcount=lonarr(nclasses)
  
;;count how many points are in each class.  
;;could be done on the fly by concatenating arrays instead?  
  print,  'Getting Class data'
  for i=0, LsInfo.nl-1 do begin
     readu, class_un, classline
     for j=0, nclasses-1 do begin
        index=where(classline eq j)
        if index(0) ne -1 then $
          classcount(j) = classcount(j) + n_elements(index)
     endfor
  endfor
  index=0
  
;;return to the begining of the file
  point_lun, class_un, 0	
;	close, class_un
;	openr, class_un, classes
  
  for i=0, offset.y-1 do begin
     readu, cent_un, centline
  endfor
  
  dataout=string(meta, '.out')
  openw, oun, /get, dataout

;  old_device=!D.name  
;  set_plot, 'ps'
;  outfile=string(meta, '_full', '.ps')
;  device,file= outfile ,xsize=7.5, ysize=10, xoff=.5, yoff=.5, /inches, $
;         set_font='Times', /tt_font, font_size=14
  
  data=fltarr(3,nclasses-1)
  print, 'starting...'
  start=(classcount(0) eq 0) + 1	;;if we screwed up and classified dark area as a class
  ;;then the first class (unclassified) will be empty, and we should
  ;;skip the second class (classcount(1)) too rare occurance, unclear.
  
  for i=start, nclasses-1 do begin	;;note we start at class 1 because
     ;; class 0 is unclassified
     if classcount(i) lt 10 then begin
        print, 'Not Enough points in class ', i
;			close, class_un, ndvi_un, cent_un, oun
;			free_lun, class_un, ndvi_un, cent_un, oun
;			return, -1
        data(*,i-1) = -1
     endif else begin
        ndvidata=make_array(classcount(i)*LsInfo.nb, type=LsInfo.type)
        centdata=fltarr(classcount(i)*CentInfo.nb)
        curcount = 0
        
        point_lun, ndvi_un, 0
        point_lun, cent_un, 0
        
        for j=0, LsInfo.nb-1 do begin
           point_lun, class_un, 0
           for k=0, LsInfo.nl-1 do begin
              
              readu, class_un, classline
              readu, ndvi_un, ndviline
              
;;if this is one in offset.factor lines then read a line from spatCent too
              if k mod offset.factor eq 0 then $
                readu, cent_un, centline
              newcount=0
              index=where(classline eq i)
              if index(0) ne -1 and newcount le n_elements(ndvidata) then begin
                 newcount=n_elements(index)+curcount
                 ndvidata(curcount:newcount-1) = ndviline(index)
                 
                 index=fix((index / offset.factor)+offset.x)
                 centdata(curcount:newcount-1) = centline(index)
                 
                 curcount=newcount
              endif else if newcount gt n_elements(ndvidata) then begin
                 print, 'ERROR, ',newcount,' is greater then the number of elements ',$
                        'in NDVI or Century Data = ', n_elements(ndvidata)
                 print, 'Class =',i,' Current line=',k
              endif
           endfor
        endfor
        
;    plot, centdata, ndvidata, psym=3
        if n_elements(centdata) eq n_elements(ndvidata) then begin
           data(0:1,i-1) = linfit(centdata, ndvidata)
           data(*,i-1) = data(*,i-1) / 255.
;			r = correlate(centdata, ndvidata/255.)
           res = regress(centdata, ndvidata/255., ftest=ftst, correlation=r)
           print, 'Coefficients ', data(0:1,i-1)
           print, 'R^2 = ', r^2
           print, 'F = ', ftst
           dfn = 1              ;degrees of freedom in the numerator of the F-statistic
           dfd = n_elements(centdata) - 2 ; degrees freedom in denominator
           p = f_pdf(ftst, dfn, dfd)
           print, 'p = ', p
           data(2,i-1) = r^2
           writeu, oun, classcount(i), fix(ndvidata), centdata
        endif else print, 'incorrect number of points', $
                          n_elements(centdata), n_elements(ndvidata)
     endelse    
     print, 'Points in this class = ', classcount(i)
     
  endfor ;;loop over n_classes
  
;  device, /close
;  set_plot, old_device
  
  close, class_un, ndvi_un, cent_un, oun
  free_lun, class_un, ndvi_un, cent_un, oun
  
  return, data
end