LAI_otv_gf_program.txt
;********************************************************************************
;PURPOSE:
;      This software estimates the Leaf Area Index from hemispherical photographs
;      using the following steps, in order:
;      1. Gets the Histogram for each hemispherical Image
;      2. Gets the OTV following Sahoo et al. 1997
;      3. Gets the LAI using the gap fraction method. Norman and Cambell 1989
;
;AUTHOR:
;      Robinson I. Negrón Juárez @ University of Sao Paulo. May 2001
;
;NOTE:
;      suggestions at lcb@model.iag.usp.br
;
;WARNINGS.
;      1. Hemispherical photographs must be gray-level images, such that
;           for a given pixel R=G=B
;      2. You must know the image dimension as well as the effective area
;           of the hemispherical photograph.
;      3. Other warnings will be updated in brief.
;      4. LAI_otv_gf.pro , tvcircle.pro and the hemis. photo must be in the same directory
;
;REVISON HISTORY:
;         RNJ: 20070818. The Software was merged (HISTOGRAMS, OTV, GAP Fraction)
;
;********************************************************************************
;
;################################################################################
;########################### FAIR USE (read carefullly)##########################
;########################### DO NOT Remove this note. ###########################
;################################################################################
;
;This software is available as-is without any warranty. You must not claim authorship of
;the original software. If you use this software or part of its concepts in your
;RESEARCH WORK an acknowledgment and citation (see below) are required.
;For purposes other than research please contact to:
;     Dr. Humberto R. da Rocha.
;     E-mail: humberto@model.iag.usp.br
;     Departamento de Ciências Atmosféricas/IAG/Universidade de São Paulo
;     Rua do Matão, 1226 - Cidade Universitária - São Paulo, SP  - Brasil
;     Cep 05508-090 -  Tel +55(11)30914713,30914705   Fax +55(11)30914714
;
;Acknowledgments:
;     Thanks to the Laboratory of Climate and Biosphere at IAG- Univerity of Sao Paulo
;     (http://www.dca.iag.usp.br/www/material/humberto/index.htm) for making available
;     its LAI-OTV-GF software supported by CAPES and FAPESP (Contract 02/09289-9).
;
;Citation:
;     Negron Juarez, R.I, H.R. da Rocha, M.A. Silva Figueira, M. L. Goulden, and
;     S. Miller, An improved estimate of leaf area index based on the histogram analysis
;     of hemispherical photographs. Agric. Forest Meteorol. (2008), 
;     doi:10.1016/j.agrformet.2008.11.012
;##################################################################################

Close, /All

;#########################################################################
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
;=================
; Controlers
;=================
;inputs and Outputs Files
;------------------------
CD, 'C:\Users\robinson\USPP003_Lai_AgrForMet\LAI_IDL_Programs\LAI_OTV_GF\'
ImgDir='C:\Users\robinson\USPP003_Lai_AgrForMet\LAI_IDL_Programs\LAI_OTV_GF\'
ImgFiles=File_Search(ImgDir,'*.bmp', Count=nimg)

FileOut='LAI_OTV_GF.txt'

;Image Dimensions: You must know it
;-----------------
NCOL=640
NFIL=480

;Effective area (the hemispheric image): You must know it
;----------------------------------------
NCini=100
NCfin=579
NFini=0
NFfin=479

;leaves Clumping factor: You must provide this information
;----------------------
ClFc=0.8

;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
;#########################################################################
;
OpenW,2, FileOut
thisformat='(A70, 1x, 5(A10, 1x))'
PrintF,2,'IMAGE', 'LAI', 'X', 'MSE', 'K', 'Tmean', Format=thisformat

For im=0, nimg-1 Do Begin
;
;=========================
;1.getting the Histograms
;=========================

   HMIN=100
   HMAX=255
   PTSHISTO=HMAX-HMIN+1

   B=IntArr(NFIL, NFIL)

;======================================================================
;	____________		________ _ _row 0   (A) the total size of Image
;   |	640		|       | 480  |			(B) the efective area
;	| 	 (A)   	|480    | (B)  |			    of the hemisph.
;	|			|		|	   |480             Photo. and radius=NFIL/2
;	|___________|		|______|_ _ row 479
;	0			639	  cl:100   cl:579		RNJ 24/04/01
;  Original Size      Area containing the Hemispherical Photograph
;======================================================================
window, xsize=NFIL, ysize=NFIL,TITLE='FOTO HEMISFERICA -- LCB | Imagem: '+ImgFiles[im]

    HISTO=FltArr(PTSHISTO)
	Image= Read_BMP(ImgFiles[im])
	R = Reform(Image(0,*,*))      ; Because in gray-level images for a given pixel: R=G=B
	B=R[NCini:NCfin, NFini:NFfin]           ; TOMO EL SECTOR CUADRADO DE LADO 2r,

	For I=0, NFIL-1 Do Begin
		For J=0, NFIL-1 Do Begin
	 		If (B(I,J) LT 10) Then B(I,J)=10 Else B(I,J)=B(I,J)
	 	EndFor
	EndFor
	BB= B

	radio=(NFIL/2.0)
	Tvcircle, radio , NFIL/2, NFIL/2, /FILL; 1=BLACK, NOTHING= WHITE
	A=TVRD(TRUE=1)
	AA=Reform(A(0,*,*))

	For I=0, NFIL-1 Do Begin
		For J=0, NFIL-1 Do Begin
		  If (AA(I,J) GE 150) Then AA(I,J)=255 ELSE AA(I,J)=0
		EndFor
	EndFor
	ARO=AA

	CC=ARO+BB

	For M=0, NFIL-1 Do Begin
		For N=0, NFIL-1 Do Begin
		 If (CC(M,N) EQ BB(M,N)) Then CC(M,N)=0 Else CC(M,N)=BB(M,N)
		EndFor
    CC=CC
	EndFor

 	DISCO=CC
 	TV, DISCO
 	HISTO(*)=Histogram(DISCO, Min=HMIN, Max=HMAX)

;==================================================
;2.getting the Optimal Threshold value
;(Sahoo and Albert 1997, Opt. Eng. 36(7), 1976-1981
;===================================================

   EPixP=FltArr(PTSHISTO)
   EPixB=FltArr(PTSHISTO)
   dif=FltArr(PTSHISTO)

   For t=0, PTSHISTO-1 Do Begin

      Pi=HISTO/(Total(HISTO, /NaN)) ; Probability of gray-level values.
      PiP=Total( Pi[0:t],/NaN )

      IF ((t+1) EQ 156) Then Begin
       PiB=0
      EndIf Else Begin
       PiB=Total( Pi[(t+1):PTSHISTO-1], /NaN )
      EndElse
	   PiB=PiB

    ; a priori Total entropy
      ET=-1*Total( Pi[0:PTSHISTO-1]*ALog(Pi[0:PTSHISTO-1]), /NaN)

    ; black (preto:P) pixels entropy
      EPixP(t)=-1*Total( (Pi[0:t]/PiP)*ALog( Pi[0:t]/PiP), /NaN )

    ; white (brancos:B) pixels entropy
      IF PiB EQ 0 Then Begin
       EPixB(t)=0
      EndIf Else Begin
       EPixB(t)=-1*Total( (Pi[t+1:PTSHISTO-1]/PiB)*ALog( Pi[t+1:PTSHISTO-1]/PiB), /NaN )
      EndElse
       EPixB(t)=EPixB(t)

       Dif(t)=(EPixP(t)-EPixB(t))^2
    EndFor

	Dif=Dif
	DifMin=Min(Dif, /NaN)
	Threshold=Where(Dif EQ DifMIn, count)
    Threshold=Threshold+100  ; My scale starts at 100 (and runs to 255)

;==============================================================
;3.getting the Leaf Area Index by using the Gap Fraction Method
;Norman, J.M., and Campbell, G.S., 1989. Canopy Structure.
;In: Pearcy R.W., Ehlesinger J., Mooney H.A. Rundel P.W. (Editors),
;Plant physiological Ecology. Field methods and
;instrumentation. Chapman and Hall, London, pp. 301-325
;==============================================================

   ZENITH=[10,20,30,40,50,60,70]
   NAZ=N_Elements(ZENITH)
   GFd= FltArr(NAZ)

   window, xsize=NFIL, ysize=NFIL,TITLE='GAP-FRACTION -  LCB IAG/USP  lcb@model.iag.usb.br'

   For ZZ=0, NAZ-1 Do Begin         ; Changing the Zenith angle
      For I=0, NFIL-1 Do Begin
      For J=0, NFIL-1 Do Begin
        IF (B(I,J) GE Threshold) Then B(I,J)=255 Else B(I,J)=10
      EndFor
      EndFor
       BB=B	                    ; Original Image with threshold

      radio=(NFIL/2.0)/15.0     ; In the CI-110 the zenith range from 0 to 75 dg
								; ASI QUE TOMO EL RADIO =75/15 => (NFIL/2)/15
								; ZENITH (+-5°) van GARDINGEN ET AL. 1999 AFM, PG.247
		Case ZENITH(ZZ) Of
			10:Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*3.0 , NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*1.0 , NFIL/2, NFIL/2,1,/FILL
			EndCase
			20:Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*5.0 , NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*3.0 , NFIL/2, NFIL/2,1,/FILL
			EndCase
			30:Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*7.0 , NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*5.0 , NFIL/2, NFIL/2,1,/FILL
			EndCase
			40:Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*9.0 , NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*7.0 , NFIL/2, NFIL/2,1,/FILL
			EndCase
			50: Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*11.0, NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*9.0 , NFIL/2, NFIL/2,1,/FILL
			EndCase
			60:Begin
			tvcircle, radio*15.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*13.0, NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*11.0, NFIL/2, NFIL/2,1,/FILL
			EndCase
			70:Begin
			tvcircle, radio*16.0, NFIL/2, NFIL/2,1,/FILL
			tvcircle, radio*15.0, NFIL/2, NFIL/2,  /FILL
			tvcircle, radio*13.0, NFIL/2, NFIL/2,1,/FILL
			EndCase
		End

		A=TVRD(TRUE=1)
		AA=Reform(A(0,*,*))

		For K=0, NFIL-1 Do Begin
			For L=0, NFIL-1 Do Begin
			 If (AA(K,L) GE 150) Then AA(K,L)=255 Else AA(K,L)=0
			EndFor
		EndFor
		ARO=AA

		CC=ARO+BB

		; --- Segment of interest
		For M=0, NFIL-1 Do Begin
			For N=0, NFIL-1 Do Begin
			 If (CC(M,N) EQ BB(M,N)) Then CC(M,N)=0 Else CC(M,N)=BB(M,N)
			EndFor
		EndFor
		DISCO=CC

		;--------Getting the Gap Fraction
		BR=Where(DISCO EQ 255, COUNT)
		 W=Double(COUNT)
		PR=Where(DISCO EQ 10, COUNT)
		 P=Double(COUNT)
		GF=W/(P+W)
		If (GF EQ 0.0) Then GF=0.000005 Else GF=GF
		GFd(ZZ)=GF
	  Erase   ; erases the current screen
   EndFor

    ;calculating the LAI

      ZENITH=[10 , 20 , 30 , 40 , 50 , 60, 70]    ; Zenith Angles 10-70
      Z=Tan(Double(ZENITH*!PI/180.0))
      NZ=7      ; number of zenith angles
      TR=DblArr(NZ)
; --------- GAP PRACTION for each AZ and LAI calculation
     TR(0)=GFd(0)
	 TR(1)=GFd(1)  & TR(2)=GFd(2)  & TR(3)=GFd(3)
	 TR(4)=GFd(4)  & TR(5)=GFd(5)  & TR(6)=GFd(6)

	For K=0, NZ-1 Do Begin     ; PARA ZENITH 10-70
	 	IF TR(K) EQ 0.0 Then TR(K)=99999 Else TR(K)=TR(K)
    EndFor
	TR=TR
    T=Double(Alog(TR))
    TMEAN=Mean(T)

;;-------- FINDING X USING THE BISECTION METHOD, Norman and Campbel, 89---------
	A=0.1 & B=10. & X=1. & DX=0.01
	While ABS(A-B) GT DX Do Begin
		S1=0.0    & S2=0.0   & S3=0.0  & S4=0.0
	  	For J=0, NZ-1 Do Begin
		TZ=Z(J)
		KB=SQRT(X^2 + TZ^2)/(X+1.774*(X+1.182)^(-0.733))
		DK=(SQRT((X+DX)^2 + TZ^2)/((X+DX)+1.774*((X+DX)+1.182)^(-0.733))-KB)
   		S1=S1+KB*T(J)
    	S2=S2+KB*KB
    	S3=S3+KB*DK
    	S4=S4+DK*T(J)
	  	EndFor
	  	F=S2*S4-S1*S3
	  	;print, X, F
		If (F LT 0.0) Then A=X Else B=X
		X=(A+B)/2.0
	EndWhile

	IF (S2 NE 0.0) Then Begin
	L=-S1/S2
	EndIf Else Begin
	L=99999
	EndElse

	LAI=L/ClFc  ;ClFc:  clumping factor

;; MSE between o measured K and calculated K
    sum_erro=0.
    Tmedido=DBLARR(NZ) & Tpredito=DBLARR(NZ) & difsqrt=DBLARR(NZ)
    	For Kk=0, NZ-1 Do Begin
    	Tmedido(Kk)=EXP(T(Kk))
    	Tpredito(Kk)=EXP(-(SQRT(X^2 + Z(Kk)^2)/(X+1.774*(X+1.182)^(-0.733)))*LAI)
    	difsqrt(Kk)=(Tpredito(Kk) - Tmedido(Kk))^2
    	sum_erro=sum_erro+difsqrt(Kk)
		EndFor
    MSE=sqrt(sum_erro/7.)
    If (MSE GT 1) Then MSE=1 Else MSE=MSE
;; printing out
     PrintF,2,Format='(A70, 1X, 5(F10.6, 1X))',ImgFiles(im),LAI, X, MSE, KB, TMEAN

   Print, ImgFiles[im], LAI
EndFor
Close, 2

Print, 'Done  :-)'
END