######################################################################
###				Einführung in spatstat				 ###
######################################################################
###         Stochastische Geometrie II   -  WS 2016/17      	 ###
###       			   Dr. Jürgen Kampf       	    		 ###
######################################################################


# spatstat ist ein R-Paket für Stochastische Geometrie und Räumliche Statistik

# Ggf. muss es zunächst aus den Internet heruntergeladen und installiert werden.
# In jedem Fall muss es geladen werden:
library(spatstat)

# Überblick
demo(spatstat)
help(spatstat)

#Simulation von Poisson-Voronoi-Mosaiken

Xt=rpoispp(100)		# Simulation eines stationären Poissonprozesses mit
plot(Xt, main ="")	# Intensität 100 in Beobachtungsfenster [0,1]x[0,1]

X=dirichlet(Xt)		# In R heißt das Voronoi-Mosaik Dirichlet-Mosaik
plot(X, main="",add=TRUE)

# Bei der Analyse von Poisson-Voronoi-Mosaiken ist doppelte Randbehandlung 
# notwendig. 
# Der Poisson-Prozess muss auch außerhalb des Beobachtungsfensters simuliert 
# werden, damit im Beobachtungsfenster ein korrektes Poisson-Voronoi-Mosaik
# angezeigt wird. 
# Für die statistischen Analysen dürfen nicht einfach die Zellen, die ganz 
# in dem Beobachtungsfenster enthalten sind, genutzt werden (erst recht nicht
# der Schnitt der Zellen mit dem Beobachtungsfenster). Statt dessen muss eine 
# Bedingung nur an die Zentrumsfunktion definiert werden, die sicherstellt, 
# dass die Zelle ganz im Fenster enthalten ist.  

Xt=rpoispp(100,win=owin(c(-1,2),c(-1,2)) ) 	# Wir simulieren den Poisson-						
plot(Xt, main ="")	# Prozess nun auf dem vergrößerten Fenster [-1,2]x[-1,2]

X=intersect.tess(dirichlet(Xt),owin(c(0,1),c(0,1)))	# Da R per default
			# das Voronoi-Mosaik im selben Fenster wie den zu Grunde
			# liegenden Punktprozess simuliert, muss dieses mit dem
plot(X, main="")	# dem Fenster [0,1]x[0,1] geschnitten werden. 

# Nun müssen wir ein Programm schreiben, dass zu einer Zelle ihr Zentrum 
# berechnet 
zentrum <- function(Cell){
	x=edges(Cell)$ends$x0;	# edges liefert die Menge der Kanten einer Zelle
	y=edges(Cell)$ends$y0;	# ends die Endpunkte dieser Kanten
	c(min(x),min(y))
	}
# Alternative: Schwerpunkt centroid.owin()


for(P in tiles(X)){
print(zentrum(P))
	points(zentrum(P)[1], zentrum(P)[2])
	}

# Die Flächen vom Schnitt von Beobachtungsfenster und Zellen
sort(tile.areas(X))
# Wie oben ausgeführt, macht es keinen Sinn diese Stichprobe zu analysieren.
# Statt dessen ermitteln wir diejenigen Zellen, deren Zentrum in 
# (0, 0.8) x (0, 0.8) liegt.
T=tiles(X)
zentral=rep(TRUE, length(T) )
for(i in ( 1: length(T) ) ){
	z= zentrum(T[[i]])
	zentral[i]= inside.owin(z[1], z[2], w=owin(c(10^{-6}, 0.8), c(10^{-6}, 0.8) ) )
	# text(z[1],z[2], round(10000*tile.areas(X)[i]))
	}
sort(tile.areas(X)[zentral])