# A script that produces a series of jpeg images which cvan be used as 
# an animation about forest canopy.
# The animation was used e.g. in my SilviLaser2010 presentation 
# http://cs.uef.fi/~lamehtat/documents/Mehtatalo.pdf
# 
# One of the figures was also used in this book chapter: 
# Mehtätalo, L., Nyblom, J., and Virolainen, A. 2014. A model-based approach for the 
# recovery of forest attributes using ALS data. Book chapter in Maltamo, M., Næsset, E. 
# and Vauhkonen, J. (editors). Forestry applications of airborne laser scanning- Concepts 
# and Case Studies. Springer, Managing Forest Ecosystems 27, pp 193-211.
# http://link.springer.com/chapter/10.1007/978-94-017-8663-8_10
# 
# Lauri Mehtätalo (published in the web on 3.3.2017) 
# http://cs.uef.fi/~lamehtat/

library(fCopulae)

# Palauttaa puun latvusalan korkeudella h,
# jonka pituus on H, latvus on ellipsin muotoinen, 
# la puoliakselit xax*H ja yax*H
# Cross-sectional area of a single tre crown 
# at height h. Total height of the tree is H.
# xax*H and yax*H are the half axes of the ellipse. 

lala.ellipsi<-function(h,H,xax=0.1,yax=0.4) {
              a<-xax*H
              b<-yax*H
              y0<-h-(1-yax)*H
              r<-a*sqrt(pmax(0,1-y0^2/b^2))
              r[h<(1-yax)*H]<-a[h<(1-yax)*H]
              A<-pi*r^2
              dA<--a/sqrt(pmax(0,1-y0^2/b^2))*(h-H+b)/b^2
              dA[h<(H-b)|h>=H]<-0
              list(a=a,b=b,y=y0,r=r,A=A,dA=dA)
              }

# gives the total height of a tree with crown radius of R at reference height h,
# crown is assumed to be ellipsoid with half axles xax*H and yax*H
Rinv.ellipsi<-function(h,R,xax=0.1,yax=0.4) {
              l<-max(length(h),length(R))
              p<-rep(xax,l)
              q<-rep(yax,l)
              h<-rep(h,length.out=l)
              R<-rep(R,length.out=l)
              h<-pmax(h,(1-q)*R/p)
 #             a<-p^2*(2*q-1)
 #             b<-2*p^2*h*(1-q)
 #             c<--q^2*R^2-p^2*h^2
 #             (-b+sqrt(b^2-4*a*c))/(2*a)
              Rinv<-(-p*(1-q)*h+q*sqrt(p^2*h^2+2*q*R^2-R^2))/(p*(2*q-1))
              dRinv<-(q-1)/(2*q-1)+p*q*h/((2*q-1)*sqrt(p^2*h^2+2*q*R^2-R^2))
              list(value=Rinv,dRinv=dRinv)
              }
              
plot(seq(0,25,0.1),sapply(seq(0,25,0.1),function(x) Rinv.ellipsi(x,2)$value),type="l")
plot(seq(0,25,0.1),sapply(seq(0,25,0.1),function(x) Rinv.ellipsi(x,2)$dRinv),type="l")


# calculates the probability that a given point x,y is not covered in a forest where 
# lambda trees per square metrer are located in the nodes of a regular square grid,
# distibution of heights is weibull,
# and crowns are ellipses. 
# The origo is set on a tree location.
# Trees from a 10*10 trees square are used.
P.noclosure.w.e.f<-function(h,x,y,lambda,shape,scale,xax=0.1,yax=0.4) {
                   dist<-xax*qweibull(0.99999999999999,shape,scale)*sqrt(lambda)
                   apu<-seq(0-floor(dist),floor(0.5+dist))/sqrt(lambda)
                   xcoords<-rep(apu,length(apu))
                   ycoords<-rep(apu,each=length(apu))
                   dst<-sqrt((xcoords-x)^2+(ycoords-y)^2)
                   prod(pweibull(Rinv.ellipsi(h,dst,xax,yax)$value,shape,scale))
                   }

# Probability of crown surface at a random point to be below h.
# asuming grid locations of trees,
# ellipsoidal crowns,
# Weibull distribution of tree heights, 
# and stand desnity lambda 
plaser.w.e.f<-function(h,lambda,shape,scale,xax=0.1,yax=0.4,crit=1e-6) {
              f<-function(x) P.noclosure.w.e.f(h,x[1],x[2],lambda,shape,scale,xax,yax)
              4*lambda*adapt(ndim=2,lower=c(0,0),upper=rep(1/(2*sqrt(lambda)),2),functn=f,eps=crit)$value
              }


             
plaser.w.e.f2<-function(h,lambda,shape,scale,xax=0.1,yax=0.4) {
               sapply(h, function(x) plaser.w.e.f(x,lambda,shape,scale,xax,yax)) 
               }
               

g.interior.w.e.f<-function(h,x,y,lambda,shape,scale,xax=0.1,yax=0.4) {
                  dist<-xax*qweibull(0.99999999999999,shape,scale)*sqrt(lambda)
                  apu<-seq(0-floor(dist),floor(0.5+dist))/sqrt(lambda)
                  xcoords<-rep(apu,length(apu))
                  ycoords<-rep(apu,each=length(apu))
                  dst<-sqrt((xcoords-x)^2+(ycoords-y)^2)
                  Rinv<-Rinv.ellipsi(h,dst,xax,yax)
                  F<-pweibull(Rinv$value,shape,scale)
                  f<-dweibull(Rinv$value,shape,scale)
                  sum( sapply(1:length(dst), function(x) prod(F[-x])*f[x]*Rinv$dRinv[x]) )
                  }

# distribution of all observations
dlaser.w.e.f<-function(h,lambda,shape,scale,xax=0.1,yax=0.4) {
              if (h>0) { 
                 f<-function(x) g.interior.w.e.f(h,x[1],x[2],lambda,shape,scale,xax,yax)
                 4*lambda*adapt(ndim=2,lower=c(0,0),upper=rep(1/(2*sqrt(lambda)),2),functn=f)$value
                 } else {
                 plaser.w.e.f(0,lambda,shape,scale,xax,yax)
                 }
              }

# distribution of all observations
dlaser.w.e.f<-function(h,lambda,shape,scale,xax=0.1,yax=0.4) {
#                 cat(h,lambda,shape,scale," ")
                 f<-function(x) g.interior.w.e.f(h,x[1],x[2],lambda,shape,scale,xax,yax)
                 if (h>(1-yax)*qweibull(0.01,shape,scale) & h<qweibull(0.99,shape,scale)) {
                    val<-4*lambda*adapt(ndim=2,lower=c(0,0),upper=rep(1/(2*sqrt(lambda)),2),functn=f,eps=0.000001)$value
                    } else if (h>(1-yax)*qweibull(0.001,shape,scale) & h<qweibull(0.9999,shape,scale)) {
                    val<-4*lambda*adapt(ndim=2,lower=c(0,0),upper=rep(1/(2*sqrt(lambda)),2),functn=f,eps=0.01)$value
                    } else {
                    val<-0
                    }
                 val
                 }

# h can be a vector              
dlaser.w.e.f2<-function(h,lambda,shape,scale,xax=0.1,yax=0.4) {
               f<-h
               f[h>0]<-sapply(h[h>0], function(x) dlaser.w.e.f(x,lambda,shape,scale,xax=xax,yax=yax))
               f[h==0]<-plaser.w.e.f(0,lambda,shape,scale,xax=xax,yax=yax)
               f
               }
               
               dlaser.w.e<-function(h,lamda,shape,scale,xax=0.1,yax=0.4) {
	       #            cat(h,lamda,shape,scale," ")
	                   lamda<-lamda/100
	       #            EY<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$A,lb,ub,rel.tol=10e-14)$value
	        
	                   q<-yax
	                   p<-xax
	                   if (h>0) {
	                      ub<-qweibull(1-1e-16,shape,scale)
	                      if (ub==Inf) print("Infinite uper bound")
	       #               lb<-(1-yax)*qweibull(1e-17,shape,scale)
	                      IG1<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$r^2,h,h/(1-yax),rel.tol=10e-14)$value
	                      IG2<-integrate(f=function(H) dweibull(H,shape,scale)*xax^2*H^2,h/(1-yax),ub,rel.tol=10e-14)$value
	                      EY<-pi*(IG1+IG2)
	       #               EY<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$A,lb,ub,rel.tol=10e-14)$value
	       
	                      s.int<--2*p^2/q^2*integrate(f=function(H) dweibull(H,shape,scale)*(h-(1-q)*H),h,h/(1-q),rel.tol=10e-14)$value
	                      -lamda*pi*s.int*exp(-lamda*EY)
	                      } else {
	                      EY0<-p^2*scale^2*gamma(1 + 2/shape)
	                      exp(-pi*lamda*EY0)
	                      }
	                   }
	       # hyväksyy vektorin h:ksi
	       dlaser.w.e2<-function(h,lamda,shape,scale,xax=0.1,yax=0.4) {
	                    sapply(h,function(x) dlaser.w.e(x,lamda,shape,scale,xax,yax))
             }
             
dlaser.w.e<-function(h,lamda,shape,scale,xax=0.1,yax=0.4) {
#            cat(h,lamda,shape,scale," ")
#            lamda<-lamda/100
#            EY<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$A,lb,ub,rel.tol=10e-14)$value
 
            q<-yax
            p<-xax
            if (h>0) {
               ub<-qweibull(1-1e-16,shape,scale)
               if (ub==Inf) print("Infinite uper bound")
#               lb<-(1-yax)*qweibull(1e-17,shape,scale)
               IG1<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$r^2,h,h/(1-yax),rel.tol=10e-14)$value
               IG2<-integrate(f=function(H) dweibull(H,shape,scale)*xax^2*H^2,h/(1-yax),ub,rel.tol=10e-14)$value
               EY<-pi*(IG1+IG2)
#               EY<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$A,lb,ub,rel.tol=10e-14)$value

               s.int<--2*p^2/q^2*integrate(f=function(H) dweibull(H,shape,scale)*(h-(1-q)*H),h,h/(1-q),rel.tol=10e-14)$value
               -lamda*pi*s.int*exp(-lamda*EY)
               } else {
               EY0<-p^2*scale^2*gamma(1 + 2/shape)
               exp(-pi*lamda*EY0)
               }
            }
# accepts vector-valued h
dlaser.w.e2<-function(h,lamda,shape,scale,xax=0.1,yax=0.4) {
             sapply(h,function(x) dlaser.w.e(x,lamda,shape,scale,xax,yax))
             }
             
plaser.w.e<-function(h,lamda,shape,scale,xax=0.1,yax=0.4) {
            EY<-integrate(f=function(H) dweibull(H,shape,scale)*lala.ellipsi(h,H,xax,yax)$A,0,qweibull(1-1e-16,shape,scale))
            p<-exp(-lamda*EY$value)
            list(EY=EY,F=p)
            }

r.ellipsi<-function(h,H,xax=0.1,yax=0.4) {
           a<-xax*H
           b<-yax*H
           y0<-h-(1-yax)*H
           r<-r2<-a*sqrt(pmax(0,1-y0^2/b^2))
           r[h<(1-yax)*H]<-a[h<(1-yax)*H]
           r2[h<(1-2*yax)*H]<-0.05*a[h<(1-2*yax)*H]
           list(r=r,r2=r2)
           }

ympyra<-function(x,y,r,col="gray",border=NA,lty="solid",lwd=1,grayfill=FALSE) {
        xapu<-sin(seq(0,pi,length=50)-pi/2)
        for (i in 1:length(x)) {
            xv1<-x[i]+xapu*r[i]
            xv2<-x[i]-xapu*r[i]
            yv1<-sqrt(pmax(0,r[i]^2-(xv1-x[i])^2))+y[i]
            yv2<--sqrt(pmax(0,r[i]^2-(xv2-x[i])^2))+y[i]
            yv<-c(yv1,yv2)
            xv<-c(xv1,xv2)
            polygon(xv,yv,border=border,col=col,lty=lty,lwd=lwd)
            }
        }
           
# gives the total height of a tree with crown radius of R at reference height h,
# crown is assumed to be ellipsoid with half axles xax*H and yax*H
Rinv.ellipsi<-function(h,R,xax=0.1,yax=0.4) {
              l<-max(length(h),length(R))
              p<-rep(xax,l)
              q<-rep(yax,l)
              h<-rep(h,length.out=l)
              R<-rep(R,length.out=l)
              h<-pmax(h,(1-q)*R/p)
 #             a<-p^2*(2*q-1)
 #             b<-2*p^2*h*(1-q)
 #             c<--q^2*R^2-p^2*h^2
 #             (-b+sqrt(b^2-4*a*c))/(2*a)
              Rinv<-(-p*(1-q)*h+q*sqrt(p^2*h^2+2*q*R^2-R^2))/(p*(2*q-1))
              dRinv<-(q-1)/(2*q-1)+p*q*h/((2*q-1)*sqrt(p^2*h^2+2*q*R^2-R^2))
              list(value=Rinv,dRinv=dRinv)
              }

plot.a.tree<-function(u,h,z,xax=0.1,yax=0.4,radius=TRUE,diam=TRUE,disc=TRUE,shade=TRUE,border=NA,col="gray",col.line="black",v=0,col.disc=NA) {
             a<-xax*h
             b<-yax*h
             h1<-seq(0,h,length=100)
             h2<-seq(h,0,length=100)
             hv<-rep(h,length(h1))
             r1<-r.ellipsi(h1,hv,xax,yax)
             r2<-r.ellipsi(h2,hv,xax,yax)
             if (disc) {
                for (i in 1:length(z)) {
                    a2<-r.ellipsi(z[i],h,xax,yax)$r
                    b2<-sin(angle)*a2
                    x<-seq(u-a2+0.0001,u+a2-0.0001,length=100)
                    y1<-z[i]+sqrt((1-(x-u)^2/a2^2)*b2^2)
                    y2<-z[i]-sqrt((1-(x-u)^2/a2^2)*b2^2)
                    x<-c(x,x[100:1])
                    y<-c(y1,y2[100:1])
                    x<-x+tan(angle)*(y-z[i])
                    polygon(x,y+v,border=NA,col=col.disc[i])
                    }
                }

             polygon(c(u-r1$r2,u+r2$r2),c(h1,h2)+v,col=col,border=border)
             if(shade) lines(c(u-r1$r,u+r2$r),c(h1,h2)+v,lty="solid")

             if (radius) {
                lines(c(u,u+a2),c(z,z)+v,lwd=2,col=col.line)
                }
             if (diam) {
                if(a2>0) lines(c(u-a2,u+a2),c(z,z)+v,lwd=2,col=col.line)
                }
             }

puut<-data.frame(x=runif(63,0,30),y=runif(63,0,30),h=rweibull(63,5,20))
puut<-puut[order(puut$y,decreasing=TRUE),]
angle<-20/90*pi/2

z1<-seq(0,27,1) 
colours<-gray(c(0.4,0.55,0.7,0.85))
F<-c()

for (j in 1:length(z1)) {
    z<-z1[1:j]
    colours<-gray(1-z/27)
    colours<-rgb(0.5-z1/100,z1/27,1-z1/27,0.4)
    jpeg(filename = paste("c:/laurim/temp/3D",j,".jpg",sep=""), width = 3*480, height = 1.5*480,
     units = "px", pointsize = 12, quality = 100, bg = "white",
     res = NA, restoreConsole = TRUE)
     
    par(xaxt="n",yaxt="n",bty="n",mfcol=c(1,2),mai=c(2,2,0.1,0.1),cex=2)
    plot(1,1,type="n",xlim=c(0,30*(1+cos(angle))),ylim=c(0,30*(1+cos(angle))),xlab=NA,ylab=NA)
    polygon(c(0,30,30*(1+cos(angle)),30*cos(angle)),c(0,0,rep(30*sin(angle),2)),border=NA,lwd=2,col=gray(0.3))

    for (i in 1:dim(puut)[1]) {
        plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=z1[1],
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc="white",xax=0.13,yax=0.4)
                   } 

    arrows(0,0,0,30)


    for (i in 1:dim(puut)[1]) {
        plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=z1[1:j],
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc=colours,xax=0.13,yax=0.4)
                   }
                   
    polygon(c(0,30,30*(1+cos(angle)),30*cos(angle)),c(0,0,rep(30*sin(angle),2))+j-1,border=NA,lwd=2,col=gray(0))
    for (i in 1:dim(puut)[1]) {
        plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=z1[j],
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc="white",xax=0.13,yax=0.4)
                   }    
    

    if (j>1) {
       for (i in 1:dim(puut)[1]) {
           plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=z1[-(1:(j-1))],
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc=colours[-(1:(j-1))],xax=0.13,yax=0.4)
                   }
       } else {
       for (i in 1:dim(puut)[1]) {
           plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=z1,
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc=colours,xax=0.13,yax=0.4)
                   }
       }
         

                   
   for (i in 1:dim(puut)[1]) {
       plot.a.tree(puut$x[i]+puut$y[i]*cos(angle),h=puut$h[i],z=c(0,10,15,20),
                   v=puut$y[i]*sin(angle),
                   col=NA,border="black",
                   shade=FALSE,
                   radius=FALSE,diam=FALSE,
                   col.disc=gray(c(0.1,0.3,0.5,0.7)),disc=FALSE,xax=0.13,yax=0.4)
                   }
#    dev.off()
    
    
    z<-z1[1:j]
    
    colours<-rgb(0.5-z/100,z/27,1-z/27)
    
     
    par(xaxt="s",yaxt="s",bty="l")
    F<-c(F,plaser.w.e(z1[j],50/900,5,20,xax=0.13,yax=0.4)$F)
    
    plot(z[1],F[1],col="black",xlab="height, m",ylab="P(Z<z) = F(z)",ylim=c(0,1),xlim=c(0,29),cex.axis=1,cex.lab=1)
    if (j>1) {
       for (i in 2:j) {
            lines(z[(i-1):i],F[(i-1):i],col=colours[i],lwd=3)
            }
       }
    dev.off()


}