Nothing
scontour=function(P,tracepoints=FALSE,colpoints="black",tracemed=TRUE,maxdepth=FALSE,xlim=c(0,2*pi),displaymed=FALSE,title="Circular Tukey contours",
ylab="Tukey's circular depth",xlab=expression(theta),colmed=2,colarc="red",sizepoints=3){
if(is.data.frame(P)) P=as.matrix(P)
if(is.list(P)){
m=length(P)
n=length(P[[1]])
y=matrix(0,n,m)
for(i in 1:m){
y[,i]=P[[i]]
if(length(P[[i]])!=n){ stop("When using a list, each element must be a vector of the same length.") }
}
P=y
}
if(is.vector(P)) p=1
if(is.matrix(P)) p=ncol(P)
if(p==1){
if(max(P)>2*pi | min(P)<0) stop("In 2D, the dataset must be a vector of angles")
contourc(P,xlim=xlim,tracepoints=tracepoints,mediane=tracemed,maxprof=maxdepth,
affichemed=displaymed,title=title,ylab=ylab,xlab=xlab,medianecol=colmed,pointscol=colpoints)
}
if(p==2){
if(max(P[,1])>2*pi | min(P[,1])<0 | max(P[,2])>pi | min(P[,2])<0) stop("Your data is interpreted as spherical coordinates of data on the sphere, but the valuies provided exceed the natural ranges for the angles : [0,2pi] and [0,pi].")
npoints=nrow(P)
points=matrix(0,nrow=npoints,ncol=3)
for(i in 1:npoints) {
points[i,1]=cos(P[i,1])*sin(P[i,2])
points[i,2]=sin(P[i,1])*sin(P[i,2])
points[i,3]=cos(P[i,2])
}
P=points
STD(P,tracepoints=tracepoints,colpoints=colpoints,colarc=colarc,sizepoints=sizepoints,tracemed=tracemed)
}
if(p==3){
if(sum(!sapply(apply(P,1,function(x){sum(x^2)}),all.equal,current=1))>0) stop("Points are not on the sphere.")
STD(P,tracepoints=tracepoints,colpoints=colpoints,colarc=colarc,sizepoints=sizepoints,tracemed=tracemed)
}
if(p>3) stop("This function supports data on the cricle or the sphere only.")
}
# Functions used by scontour
contourc<-function(P,xlim=c(0,2*pi),tracepoints=FALSE,mediane=FALSE,maxprof=FALSE,
affichemed=FALSE,title="Contour circulaire de Tukey",
ylab="Profondeur circulaire de Tukey",xlab=expression(theta),medianecol=2,
pointscol=1)
#Cette fonction trace le graphe complet de la profondeur de Tukey sur
#le cercle par rapport a un echantillon de n points en position gene-
#rale.
#Entree obligatoire :
#-Un vecteur P contenant un echantillon de n points en coordonnees
#polaires. le nombre de points doit etre superieur ou egale a 2.
#Entrees optionnelles:
#-tracepoints=T: ajoute les points de l'echantillon sur le graphique a
#la hauteur de leur profondeur respective.
#-mediane=T: ajoute le point median sur le graphique a la hauteur de
#sa profondeur.
#-maxprof=T: affiche la valeur de la profondeur maximale.
#-affichemed=T: affiche la valeur du point median.
#-title="titre": ajoute le titre souhaite au graphique
#-ylab="axe des y": ajoute le titre souhaite pour l'axe des ordonnees.
#-xlab="axe des x": ajoute le titre souhaite pour l'axe des abcisses.
#-medianecol="couleur": appose la couleur souhaitee au point median si
#mediane=T.
#-pointscol="couleur": appose la couleur souhaitee pour les points de
#l'echantillon apparaissant sur le graphique.
#Sortie :
#-Le graphe complet de la profondeur de Tukey
{
points=sort(P)
npoints=length(points)
tukdepthpoints=rep(0,npoints)
bornes=rep(0,npoints)
for(i in 1:npoints)
{
tukdepthpoints[i]=tukdepthc3(points,points[i])[[2]]
}
if(tukdepthpoints[2]<tukdepthpoints[1])
{
bornes[1]=2
}
if(tukdepthpoints[npoints]<tukdepthpoints[1])
{
bornes[1]=1
}
if(tukdepthpoints[2]<tukdepthpoints[1] & tukdepthpoints[npoints]<
tukdepthpoints[1])
{
bornes[1]=3
}
if(tukdepthpoints[1]<tukdepthpoints[npoints])
{
bornes[npoints]=2
}
if(tukdepthpoints[(npoints-1)]<tukdepthpoints[npoints])
{
bornes[npoints]=1
}
if(tukdepthpoints[1]<tukdepthpoints[npoints] &
tukdepthpoints[npoints-1]<tukdepthpoints[npoints])
{
bornes[npoints]=3
}
if(npoints>2)
{
for(j in 2:(npoints-1))
{
if(tukdepthpoints[j+1]<tukdepthpoints[j])
{
bornes[j]=2
}
if(tukdepthpoints[j-1]<tukdepthpoints[j])
{
bornes[j]=1
}
if(tukdepthpoints[j+1]<tukdepthpoints[j] & tukdepthpoints[j-1]<
tukdepthpoints[j])
{
bornes[j]=3
}
}
}
maximumprof=max(tukdepthpoints)
minimumprof=min(tukdepthpoints)
if(minimumprof==maximumprof)
{
ylim=c(0,0.55)
}
else
{
if(range.circular(circular(points))[[1]]<=pi)
{
ylim=c(0,maximumprof+2*(maximumprof-minimumprof)/20)
}
else
{
ylim=c(minimumprof,maximumprof+2*(maximumprof-minimumprof)/20)
}
}
if(tracepoints==FALSE)
{
plot(NULL,NULL,xlim=xlim,ylim=ylim,ylab=ylab,xlab=xlab,pch=16,
col=pointscol)
}
else
{
plot(points,tukdepthpoints,xlim=xlim,ylim=ylim,ylab=ylab,xlab=xlab,
pch=16,col=pointscol)
}
title(title)
maximumprof=max(tukdepthpoints)
minimumprof=min(tukdepthpoints)
v=which(tukdepthpoints==maximumprof)
if(length(v)==2 && (v[1]!=1 || v[2]!=npoints))
{
lines(c(points[min(v)],points[max(v)]),c(tukdepthpoints[min(v)],
tukdepthpoints[min(v)]))
}
if(length(v)==2 && v[1]==1 && v[2]==npoints)
{
lines(c(points[npoints],2*pi),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
lines(c(0,points[1]),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
}
if(bornes[1]==0)
{
lines(c(0,points[1]),c(tukdepthpoints[1],tukdepthpoints[1]))
lines(c(points[npoints],2*pi),c(tukdepthpoints[1],
tukdepthpoints[1]))
lines(c(points[1],points[2]),c(tukdepthpoints[1],tukdepthpoints[1]))
}
if(bornes[1]==1)
{
lines(c(0,points[1]),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
lines(c(points[npoints],2*pi),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
points(points[1],tukdepthpoints[1],pch=16,col=pointscol)
points(points[1],tukdepthpoints[npoints])
}
if(bornes[1]==2)
{
lines(c(points[1],points[2]),c(tukdepthpoints[2],tukdepthpoints[2]))
points(points[1],tukdepthpoints[1],pch=16,col=pointscol)
points(points[1],tukdepthpoints[2])
}
if(bornes[1]==3)
{
lines(c(0,points[1]),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
lines(c(points[1],points[2]),c(tukdepthpoints[2],tukdepthpoints[2]))
lines(c(points[npoints],2*pi),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
points(points[1],tukdepthpoints[1],pch=16,col=pointscol)
points(points[1],tukdepthpoints[2])
}
if(bornes[npoints]==0)
{
lines(c(points[(npoints-1)],points[npoints]),
c(tukdepthpoints[npoints],tukdepthpoints[npoints]))
lines(c(points[npoints],2*pi),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
lines(c(0,points[1]),c(tukdepthpoints[npoints],
tukdepthpoints[npoints]))
}
if(bornes[npoints]==1)
{
lines(c(points[(npoints-1)],points[npoints]),
c(tukdepthpoints[(npoints-1)],tukdepthpoints[(npoints-1)]))
points(points[npoints],tukdepthpoints[npoints],
pch=16,col=pointscol)
points(points[npoints],tukdepthpoints[(npoints-1)])
}
if(bornes[npoints]==2)
{
lines(c(points[npoints],2*pi),c(tukdepthpoints[1],
tukdepthpoints[1]))
lines(c(0,points[1]),c(tukdepthpoints[1],tukdepthpoints[1]))
points(points[npoints],tukdepthpoints[npoints],pch=16,
col=pointscol)
points(points[npoints],tukdepthpoints[1])
}
if(bornes[npoints]==3)
{
lines(c(points[(npoints-1)],points[npoints]),
c(tukdepthpoints[(npoints-1)],tukdepthpoints[(npoints-1)]))
lines(c(points[npoints],2*pi),c(tukdepthpoints[1],
tukdepthpoints[1]))
lines(c(0,points[1]),c(tukdepthpoints[1],tukdepthpoints[1]))
points(points[npoints],tukdepthpoints[npoints],pch=16,
col=pointscol)
points(points[npoints],tukdepthpoints[1])
}
if(npoints>2)
{
for(k in 2:(npoints-1))
{
if(bornes[k]==0)
{
lines(c(points[k-1],points[k]),c(tukdepthpoints[k],
tukdepthpoints[k]))
lines(c(points[k],points[k+1]),c(tukdepthpoints[k],
tukdepthpoints[k]))
}
if(bornes[k]==1)
{
lines(c(points[k-1],points[k]),c(tukdepthpoints[k-1],
tukdepthpoints[k-1]))
points(points[k],tukdepthpoints[k],pch=16,col=pointscol)
points(points[k],tukdepthpoints[k-1])
}
if(bornes[k]==2)
{
lines(c(points[k],points[k+1]),c(tukdepthpoints[k+1],
tukdepthpoints[k+1]))
points(points[k],tukdepthpoints[k],pch=16,col=pointscol)
points(points[k],tukdepthpoints[k+1])
}
if(bornes[k]==3)
{
lines(c(points[k-1],points[k]),c(tukdepthpoints[k-1],
tukdepthpoints[k-1]))
lines(c(points[k],points[k+1]),c(tukdepthpoints[k+1],
tukdepthpoints[k+1]))
points(points[k],tukdepthpoints[k],pch=16,col=pointscol)
points(points[k],tukdepthpoints[k+1])
}
}
}
if(((mediane==TRUE || affichemed==TRUE) & (minimumprof!=maximumprof ||
npoints==2)) || maxprof==TRUE)
{
minimummed=min(which(tukdepthpoints==maximumprof))
maximummed=max(which(tukdepthpoints==maximumprof))
pointsbon=points[round(which(tukdepthpoints==maximumprof),14)]
pointmax=max(pointsbon)
pointmin=min(pointsbon)
if((pointmax-pointmin)<pi)
{
pointmilieu=(pointmax+pointmin)/2
}
else
{
premierpoint=max(pointsbon[which(pointsbon<pi)])
deuxiemepoint=min(pointsbon[which(pointsbon>pi)])
pointmilieu=(premierpoint+deuxiemepoint)/2+pi
if(pointmilieu>2*pi)
{
pointmilieu=pointmilieu-2*pi
}
}
tukmed=pointmilieu
if(mediane==TRUE & (minimumprof!=maximumprof || npoints==2))
{
points(pointmilieu,maximumprof,pch=8,col=medianecol)
}
if(maxprof==TRUE & affichemed==FALSE)
{
if(maximumprof==minimumprof)
{
text((xlim[1]+xlim[2])/2,0.54,"maximum depth:")
text((xlim[1]+xlim[2])/2,0.52,tukdepthpoints[minimummed])
}
else
{
text((xlim[1]+xlim[2])/2,
maximumprof+2*(maximumprof-minimumprof)/20,
"pronfondeur maximale:")
text((xlim[1]+xlim[2])/2,
maximumprof+(maximumprof-minimumprof)/20,
tukdepthpoints[minimummed])
}
}
if(maxprof==FALSE & affichemed==TRUE)
{
if(maximumprof==minimumprof & npoints==2)
{
text((xlim[1]+xlim[2])/2,0.54,"median:")
text((xlim[1]+xlim[2])/2,0.52,pointmilieu)
}
else
{
text((xlim[1]+xlim[2])/2,
maximumprof+2*(maximumprof-minimumprof)/20,"median:")
text((xlim[1]+xlim[2])/2,
maximumprof+(maximumprof-minimumprof)/20,pointmilieu)
}
}
if(maxprof==TRUE & affichemed==TRUE)
{
if(maximumprof==minimumprof)
{
text((xlim[2]-xlim[1])/4+xlim[1],0.54,"maximum depth:")
text((xlim[2]-xlim[1])/4+xlim[1],0.52,
tukdepthpoints[minimummed])
if(npoints==2)
{
text(3*(xlim[2]-xlim[1])/4+xlim[1],0.54,"median:")
text(3*(xlim[2]-xlim[1])/4+xlim[1],0.52,pointmilieu)
}
}
else
{
text((xlim[2]-xlim[1])/4+xlim[1],
maximumprof+2*(maximumprof-minimumprof)/20,
"pronfondeur maximale:")
text((xlim[2]-xlim[1])/4+xlim[1],
maximumprof+(maximumprof-minimumprof)/20,
tukdepthpoints[minimummed])
text(3*(xlim[2]-xlim[1])/4+xlim[1],
maximumprof+2*(maximumprof-minimumprof)/20,"median:")
text(3*(xlim[2]-xlim[1])/4+xlim[1],
maximumprof+(maximumprof-minimumprof)/20,pointmilieu)
}
}
}
if((minimumprof==maximumprof) & (affichemed==TRUE || mediane==TRUE) &
npoints>2)
{
stop("All points have the same depth. Tukey's median cannot be identified.")
}
differences=rep(0,npoints)
ecart=0
for(i in 1:(npoints-1))
{
if(points[i+1]-points[i]>=pi)
{
ecart=i
}
}
if((2*pi-(points[npoints]-points[1]))>=pi)
{
ecart=npoints
}
if(ecart!=0)
{
if(ecart<npoints)
{
lines(c(points[ecart],points[ecart+1]),c(minimumprof,
minimumprof),col="white")
lines(c(points[ecart],points[ecart+1]),c(0,0))
points(c(points[ecart],points[ecart+1]),c(0,0))
points(c(points[ecart],points[ecart+1]),c(minimumprof,
minimumprof),pch=16)
}
else
{
lines(c(points[ecart],2*pi),c(minimumprof,minimumprof),
col="white")
lines(c(0,points[1]),c(minimumprof,minimumprof),col="white")
lines(c(points[ecart],2*pi),c(0,0))
lines(c(0,points[1]),c(0,0))
points(c(points[ecart],points[1]),c(0,0))
points(c(points[ecart],points[1]),c(minimumprof,minimumprof),
pch=16)
}
}
}
STD=function(P,tracepoints=FALSE,colpoints="black",colarc="red",
sizepoints=3,tracemed=TRUE)
#Cette fonction trace tous les contours de Tukey a partir d'un
#echantillon sur la sphere en position generale.
#Entree obligatoire:
#-Une matrice dont chaque ligne represente un point de l'echantillon
#donne en coordonnees cartesiennes.
#Entrees optionelles:
#1:tracepoints=T pour ajouter les points de l'echantillon sur le
#graphique.
#2:colpoints="couleur" controle la couleur des points.
#3:colarc="couleur" controle la couleur des arcs de cercles formant les
#cotes des polygones spheriques.
#4:sizepoints controle la grosseur des points.
#5:tracemed=T permet d'ajouter la mediane au graphique.
#Sorties: Une liste a 2 composantes:
#1:Un vecteur contenant la liste des profondeurs des regions tronquees
#dans l'ordre croissant.
#2:Une liste dont chaque composante est une matrice contenant les
#sommets de chaque region tronquee.
{
npoints=nrow(P)
tukdepths=NA
pointshemi=pointshemi(P)
nombrepoints=pointshemi[[1]]
normale=pointshemi[[3]]
demi=pointshemi[[4]]
for(i in 1:npoints)
{
tukdepths[i]=tukdepthsdemi(P,demi,nombrepoints,normale,P[i,])
}
tuk=sort(round(unique(tukdepths),14))
if(length(tuk)==1)
{
if(round(max(nombrepoints)==npoints,14)||round(min(nombrepoints)==2,
14))
{
rgl.open()
rgl.spheres(0,0,0,1)
B=traceregion(P,tuk,rgl=TRUE,open=TRUE,tracepoints=FALSE,
colpoints=colpoints,colarc=colarc,sizepoints=sizepoints,
demicalcul=TRUE,pointshemi=pointshemi,tukdepthscalcul=TRUE,
tukdepths=tukdepths)
if(tracepoints==TRUE)
{
rgl.points(P[,1],P[,2],P[,3],size=sizepoints,col=colpoints)
}
return(list(tuk,B[[1]]))
}
stop("All points have the same depth. The cloud cannot be split using Tukey's depth.")
return(0)
}
rgl.open()
rgl.spheres(0,0,0,1)
region=NA
comptreg=0
comptreg=1
npointsregions=NA
profondeur=tuk[1]+1/(npoints)
B=traceregion(P,profondeur,open=TRUE,rgl=TRUE,tracepoints=FALSE,
colpoints=colpoints,colarc=colarc,sizepoints=sizepoints,demicalcul=TRUE,
pointshemi=pointshemi,tukdepthscalcul=TRUE,tukdepths=tukdepths,
derniereregion=FALSE)
while(B[[3]]==1)
{
region[comptreg]=profondeur
profondeur=profondeur+1/(npoints)
pointsregion=B[[1]]
if((is.matrix(pointsregion))==0)
{
npointsregion=1
}
else
{
npointsregion=nrow(pointsregion)
}
faces=B[[2]]
norms=B[[4]]
if(comptreg==1)
{
sommets=pointsregion
npointsregions[1]=npointsregion
}
else
{
sommets=rbind(sommets,pointsregion)
npointsregions[comptreg]=npointsregion
}
B=traceregion(P,profondeur,rgl=TRUE,open=TRUE,tracepoints=FALSE,
colpoints=colpoints,colarc=colarc,sizepoints=sizepoints,
demicalcul=TRUE,pointshemi=pointshemi,tukdepthscalcul=TRUE,
tukdepths=tukdepths,derniereregion=TRUE,normsdern=norms)
comptreg=comptreg+1
}
nreg=comptreg-1
if(tracepoints==TRUE)
{
rgl.points(P[,1],P[,2],P[,3],size=sizepoints,col=colpoints)
}
Som=as.list(rep(NA,nreg))
if(npointsregions[1]==1)
{
Som[[1]]=sommets
}
else
{
Som[[1]]=sommets[1:npointsregions[1],]
}
sum=0
if(nreg>=2)
{
for(i in 2:nreg)
{
sum=sum+npointsregions[i-1]
Som[[i]]=sommets[c((sum+1):(sum+npointsregions[i])),]
}
}
if(tracemed==TRUE)
{
if(is.matrix(Som[[nreg]])==0)
{
rgl.points(Som[[nreg]][1],Som[[nreg]][2],Som[[nreg]][3],size=3,
col="green",pch=8)
}
else
{
med=spherpolycentroide(Som[[nreg]],faces)
rgl.points(med[1],med[2],med[3],size=sizepoints,col="green",pch=8)
}
}
return(list(region,Som,med))
}
espacetronque=function(P,alpha,demicalcul=FALSE,pointshemi=NA,
tukdepthscalcul=FALSE,tukdepths=NA,derniereregion=FALSE,normsdern=NA)
#Cette fonction calcule le polygone spherique formant la frontiere d'une
#region tronquee par rapport a un echantillon en position generale sur
#la sphere a l'aide de la notion de profondeur spherique empirique de
#Tukey.
#Entrees obligatoires:
#-Une matrice decrivant l'echantillon. Chaque ligne est un point
#de la sphere en coordonnees cartesiennes. Le nombre de points doit
#etre d'au moins 3.
#-Un nombre "alpha" representant la borne inferieure de la profondeur
#de la egion tronquee que l'on souhaite obtenir. alpha doit etre de la
#forme k/n ou n est la taille de l'echantillon.
#Entrees optionnelles:
#-demicalcul=T indique que pointshemi a deja ete calcule. La valeur F
#est donnee par defaut.
#-Dans le cas ou demicalcul=T, pointshemi doit contenir une liste a
#quatre composantes qui est calculee par pointshemi (voir la fonction
#pointshemi).
#-tukdepthscalcul=T indiquent si la profondeur en chaque points de
#l'echantillon a ete calculee.
#-Dans le cas ou tukdepthscalcul=T, tukdepths doit etre un vecteur ou
#l'entree j est la profondeur de Tukey au j ieme point de
#l'echantillon.
#-derniereregion=T indique que le calcul d'une region tronquee de
#profondeur inferieur a deja ete cacule.
#-Dans le cas ou derniereregion=T, normsdern doit contenir la matrice
#contenant les vecteurs normaux associes aux cotes du polygone
#spherique correspondant a la region tronquee anterieure.
#Les 4 dernieres options sont utilises par defaut par la fonction STD()
#et il n'est pas essentiel pour l'utilisateur de bien comprendre leur
#fonctionnement.
#Sorties:
#-Une liste a 4 composantes:
#1:Une matrice de points de la sphere. Ce sont les sommets du polygone
#spherique.
#2:Une matrice de deux colonnes dont chaque ligne sont les deux indices
#des points de la matrice de la sortie 1 formant un cote du polygone
#spherique.
#3:Un scalaire valant 0,1 ou 2 selon le cas ou l'espace tronque est
#vide, un polygone spherique ou la sphere au complet, respectivement.
#4:La matrice des vecteurs normaux associes a chaque face du polygone
#spherique.
{
npoints=nrow(P)
alpha=alpha*npoints
alpha=ceiling(round(alpha,14))-(1)*(ceiling(round(alpha,14))-
alpha>=1/2)
nombrepointsmin=alpha+1
if(demicalcul==FALSE)
{
pointshemi=pointshemi(P)
}
G=pointshemi[[4]]
nombrepoints1=pointshemi[[1]]
nombrepoints2=pointshemi[[2]]
normale=pointshemi[[3]]
if(tukdepthscalcul==FALSE)
{
tukdepths=NA
for(i in 1:npoints)
{
tukdepths[i]=tukdepthsdemi(P,G,nombrepoints1,normale,P[i,])
}
}
bonddemis1=which(nombrepoints1<=nombrepointsmin)
bonddemis2=which(nombrepoints2<=nombrepointsmin)
nbind=nrow(G)
bonpoints=NA
comptbonp=0
bonpoints2=NA
comptbonp2=0
for(i in 1:npoints)
{
if(abs(tukdepths[i]-alpha/npoints)<=(10^(-12)))
{
comptbonp=comptbonp+1
bonpoints[comptbonp]=i
}
if(tukdepths[i]>=(alpha/npoints-(10^(-12))))
{
comptbonp2=comptbonp2+1
bonpoints2[comptbonp2]=i
}
}
if(comptbonp2==0)
{
return(list(NA,NA,0,NA))
}
if((comptbonp2==npoints && max(c(nombrepoints1,nombrepoints2))!=
npoints)
||(alpha==0))
{
return(list(NA,NA,2,NA))
}
if(comptbonp2==1)
{
return(list(P[bonpoints2,],1,1,NA))
}
normale=(-1)*normale
bonddemis=c(bonddemis1,bonddemis2)
pointsintersection=matrix(NA,nrow=1,ncol=3)
verification=rep(0,nbind)
ndemi=length(bonddemis)
type1=NA
type2=NA
comptint=0
for(k in 1:(ndemi-1))
{
for(l in (k+1):ndemi)
{
inter=0
verif=0
if(G[bonddemis[k],][1]==G[bonddemis[l],][1])
{
inter=1
v=c(1,1)
}
if(G[bonddemis[k],][2]==G[bonddemis[l],][2])
{
inter=1
v=c(2,2)
}
if(G[bonddemis[k],][1]==G[bonddemis[l],][2])
{
inter=1
v=c(1,2)
}
if(G[bonddemis[k],][2]==G[bonddemis[l],][1])
{
inter=1
v=c(2,1)
}
if(inter==1)
{
point1=P[G[bonddemis[k],v[1]],]
point2=-point1
verification[G[bonddemis[k],v[1]]]=verification[G[bonddemis[k],
v[1]]]+1
if(verification[G[bonddemis[k],v[1]]]>=2)
{
verif=1
}
}
else
{
points=intersect(normale[bonddemis[k],],normale[bonddemis[l],])
point1=points[[1]]
point2=points[[2]]
}
bl=1
if(comptint==0)
{
bl=0
comptint=comptint+2
pointsintersection=rbind(point1,point2)
}
if(comptint>=1&&bl==1)
{
if(verif==0)
{
pointsintersection=rbind(pointsintersection,rbind(point1,
point2))
comptint=comptint+2
}
}
}
}
moybonp1=sum(P[bonpoints2,1])/comptbonp2
moybonp2=sum(P[bonpoints2,2])/comptbonp2
moybonp3=sum(P[bonpoints2,3])/comptbonp2
moybonp=c(moybonp1,moybonp2,moybonp3)
moybonp=moybonp/(sqrt(sum(moybonp^2)))
mauvais=NA
comptmauv=0
if(derniereregion==TRUE)
{
for(j in 1:comptint)
{
if(verifextpoly(normsdern,pointsintersection[j,])==1)
{
comptmauv=comptmauv+1
mauvais[comptmauv]=j
}
else
{
dev=pointsintersection[j,]+0.00000001*(moybonp-
pointsintersection[j,])
dev=dev/(sqrt(sum(dev^2)))
if(tukdepthsdemiineg(P,G,nombrepoints1,-normale,dev,
alpha/npoints)==1)
{
comptmauv=comptmauv+1
mauvais[comptmauv]=j
}
}
}
}
else
{
for(j in 1:comptint)
{
dev=pointsintersection[j,]+0.00000001*(moybonp-
pointsintersection[j,])
dev=dev/(sqrt(sum(dev^2)))
if(tukdepthsdemiineg(P,G,nombrepoints1,-normale,dev,
alpha/npoints)==1)
{
comptmauv=comptmauv+1
mauvais[comptmauv]=j
}
}
}
pointsintersection=pointsintersection[-mauvais,]
if(comptbonp>0 && comptbonp<=2)
{
pointsintersection=rbind(pointsintersection,P[bonpoints,])
pointsintersection=unique(round(pointsintersection,12))
}
W=schema.polyspher(pointsintersection)
faces=W[[1]]
norms=W[[2]]
return(list(pointsintersection,faces,1,norms))
}
pointshemi=function(P)
#La fonction pointshemi calcule le nombre de points de l'echantillon
#contenus dans chacun des hemispheres dont la frontiere est un grand
#cercle passant par deux points de l'echantillon.
#Entrees:
#-Une matrice dont chaque ligne contient un point de l'echantillon en
#coordonnees cartesiennes.
#Sorties:
#-Une liste a 4 composantes.
#1:Un vecteur dont la j ieme composante est le nombre de points contenus
#dans l'hemisphere dont la frontiere est le grand cercle j et se situant
#du cote du vecteur normal j.
#2:Un vecteur dont la j ieme composante est le nombre de points contenus
#dans l'hemisphere dont la frontiere est le grand cercle j et se situant
#du cote inverse du vecteur normal j.
#3:Une matrice dont la ligne j est un vecteur normal du plan contenant
#le grands cercles j.
#4:Une matrice (2 parmi n) x 2 contenant les 2 parmi n combinaisons
#formant les grands cercles delimitant les hemispheres.
{
npoints=nrow(P)
G=combn(1:npoints,2)
G=t(G)
ncomb=nrow(G)
nombrepoints1=rep(2,ncomb)
nombrepoints2=rep(2,ncomb)
normale=matrix(NA,ncol=3,nrow=ncomb)
for(i in 1:ncomb)
{
demi=P[G[i,],]
normale[i,]=prodvect2(rbind(demi,c(0,0,0)))
points=c(1:npoints)[-G[i,]]
for(k in points)
{
if(normale[i,]%*%P[k,]>0)
{
nombrepoints1[i]=nombrepoints1[i]+1
}
if(normale[i,]%*%P[k,]<0)
{
nombrepoints2[i]=nombrepoints2[i]+1
}
if(normale[i,]%*%P[k,]==0)
{
stop("There is a problem.")
}
}
}
return(list(nombrepoints1,nombrepoints2,normale,G))
}
verifextpoly=function(normales,point)
{
nfaces=nrow(normales)
for(i in 1:nfaces)
{
if((normales[i,]%*%point)<0)
{
return(1)
}
}
return(0)
}
tukdepthsdemi=function(P,demi,nombrepoints,normale,theta)
#Cette fonction calcule la profondeur de Tukey lorsque le nombre de
#points contenue dans chaque hemisphere (formee a partir de deux
#points de l'echantillon) a deja ete calcule.
#Entrees:
#1:P contient les points de l'echantillons
#2:demi contient les indices des points (en rapport avec P) a partir
#desquels sont formes les grands cercles formants le bord des
#hemispheres.
#3:nombrepoints est un vecteur dont l'entree j donne le nombre de
#points que contient l'hemisphere j.
#4:normale contient une matrice dont la ligne j contient un vecteur
#normal indiquant le sens de l'hemisphere j.
#5:theta est le point de la sphere auquel on veut calculer la
#profondeur.
#Sortie: La profondeur de Tukey.
{
npoints=nrow(P)
ndemi=nrow(demi)
nombrepointst=rep(NA,ndemi)
for(i in 1:ndemi)
{
if(round((theta[1]==P[demi[i,1],1]&&theta[2]==P[demi[i,1],2]&&
theta[3]==P[demi[i,1],3])||(theta[1]==P[demi[i,2],1]&&
theta[2]==P[demi[i,2],2]&&theta[3]==P[demi[i,2],3]),14))
{
nombrepointst[i]=min(nombrepoints[i]-2+1,npoints-nombrepoints[i]
+1)
}
else
{
if((normale[i,]%*%theta)>0)
{
nombrepointst[i]=nombrepoints[i]-2
}
else
{
nombrepointst[i]=npoints-nombrepoints[i]
}
}
}
tukdepthsdemi=min(nombrepointst)/npoints
}
tukdepthsdemiineg=function(P,demi,nombrepoints,normale,theta,alpha)
#Cette fonction verifie si la profondeur de Tukey est plus petite que
#alpha lorsque le nombre de points contenue dans chaque hemisphere a
#deja ete calcule.
#Retourne 1 si oui et 0 sinon.
#Entrees:
#1:P contient les points de l'echantillons
#2:demi contient les indices des points (en rapport avec P) a partir
#desquels sont formes les grands cercles formants le bord des
#hemispheres.
#3:nombrepoints est un vecteur dont l'entree j donne le nombre de
#points que contient l'hemisphere j.
#4:normale contient une matrice dont la ligne j contient un vecteur
#normal indiquant le sens de l'hemisphere j.
#5:theta est le point de la sphere auquel on veut tester la
#profondeur.
#6:critere de comparaison pour la profondeur de Tukey.
#Sorties: 1 ou 0.
{
npoints=nrow(P)
ndemi=nrow(demi)
nombrepointst=rep(NA,ndemi)
for(i in 1:ndemi)
{
if(round((theta[1]==P[demi[i,1],1]&&theta[2]==P[demi[i,1],2]&&
theta[3]==P[demi[i,1],3])||(theta[1]==P[demi[i,2],1]&&
theta[2]==P[demi[i,2],2]&&theta[3]==P[demi[i,2],3]),14))
{
nombrepointst[i]=min(nombrepoints[i]-2+1,npoints-nombrepoints[i]
+1)
}
else
{
if((normale[i,]%*%theta)>0)
{
nombrepointst[i]=nombrepoints[i]-2
}
else
{
nombrepointst[i]=npoints-nombrepoints[i]
}
}
if(round(nombrepointst[i]<(alpha*npoints),14))
{
return(1)
}
}
return(0)
}
intersect=function(N1,N2)
#Cette fonction calcule les deux points formes par l'intersection
#de la sphere unite et de deux plans secants passant par l'origine.
#Entrees:-N1 et N2 doivent etre les vecteurs normaux respectifs des
#plans.
#Sorties: Une liste a 3 composantes.
#1: Le premier point d'intersection.
#2: Le deuxieme point d'intersection.
#3: Le nombre 1.
{
vect=prodvect2(rbind(N1,N2,c(0,0,0)))
vect=vect/(sqrt(sum(vect^2)))
point1=vect
point2=(-1)*vect
return(list(point1,point2,1))
}
tracearc=function(P1,P2,colarc="red")
#La fonction trace l'arc de grand cercle reliant les points P1 et P2.
#Entrees:
#-Les deux points.
#-La couleur de l'arc.
{
vperp=P2-((t(P1)%*%P2)/(sqrt(t(P1)%*%P1)))*P1
vperp=vperp/(sqrt(vperp%*%vperp))
alpha=acos((t(P1)%*%P2)/(sqrt(t(P1)%*%P1)*sqrt(t(P2)%*%P2)))
theta=seq(from=0,to=alpha,by=0.001)
P=matrix(0,nrow=length(theta),ncol=3)
for(i in 1:length(theta))
{
P[i,]=P1*cos(theta[i])+vperp*sin(theta[i])
}
rgl.linestrips(P[,1],P[,2],P[,3],size=2,col=colarc)
}
schema.polyspher=function(P)
#Cette fonction calule les cotes d'un polygone spherique a partir
#des sommets.
#Entree:
#-Une matrice P contenant les sommets du polygone spherique.
#Sortie:
#-Une matrice ayant deux colonnes dont chaque ligne contient le numero
#des sommets formant une face du polygone spherique.
{
npoints=nrow(P)
G=combn(c(1:npoints),2)
G=t(G)
ncomb=nrow(G)
compteur=0
for(i in 1:ncomb)
{
arc=P[G[i,],]
normale=prodvect2(rbind(arc,c(0,0,0)))
points=c(1:npoints)[-G[i,]]
demi1=0
demi2=0
for(j in points)
{
if(0<(normale%*%P[j,]))
{
demi1=1
}
else
{
demi2=1
}
}
if(demi1*demi2==0)
{
compteur=compteur+1
if(compteur>1)
{
faces=rbind(faces,G[i,])
}
else
{
faces=G[i,]
}
ajust=1
if(demi2==1)
{
ajust=-1
}
if(compteur>1)
{
matnormale=rbind(matnormale,ajust*normale)
}
else
{
matnormale=ajust*normale
}
}
}
return(list(faces,matnormale))
}
tracepolyspher=function(P,faces,colarc="red")
#Cette fonction trace un polygone spherique de la couleur desiree
#a partir des sommets et des faces d'un polygone spherique.
#Une fenetre graphique du type "rgl" avec la sphere doit etre
#ouverte prealablement, a l'aide des lignes de codes >rgl.open(),
#>rgl.spheres(0,0,0,1)
#Entrees:
#-Une matrice contenant les sommets du polygones.
#-Une matrice contenant les faces.
#-La couleur desiree.
#Sortie:
#-Le graphique du polygone spherique.
{
nfaces=nrow(faces)
for(i in 1:nfaces)
{
tracearc(P[faces[i,1],],P[faces[i,2],],colarc=colarc)
}
}
traceregion=function(P,alpha,tracepoints=FALSE,rgl=FALSE,open=FALSE,
colpoints="black",colarc="red",sizepoints=3,demicalcul=FALSE,pointshemi=NA,
tukdepthscalcul=FALSE,tukdepths=NA,derniereregion=FALSE,normsdern=NA)
#Cette fonction trace le polygone spherique decoupant une region
#tronquee par la profondeur de Tukey sur la sphere.
#Entrees obligatoires:
#-Une matrice contenant les points de l'echantillon.
#-La profondeur minimale de l'espace tronquee.
#Entrees optionnelles:
#-tracepoints=F determine si l'on souhaite faire apparaitre les points
#de l'echantillon sur la sphere.
#-rgl=F determine si une fenetre rgl a deja ete preparee ou pas.
#-colpoints="black" determine la couleur des points de l'echantillon
#(si tracepoints=T)
#-colarc="red" determine la couleur des cotes des polygones spheriques
#apparaissant.
#-sizepoints=3 determine la grosseur des points.
#-demicalcul=T indique que pointshemi a deja ete calcule. La valeur F
#est donnee par defaut.
#-Dans le cas ou demicalcul=T, pointshemi doit contenir une liste a
#quatre composantes qui est calculee par pointshemi (voir la fonction
#pointshemi).
#-tukdepthscalcul=T indiquent si la profondeur en chaque points de
#l'echantillon a ete calculee.
#-Dans le cas ou tukdepthscalcul=T, tukdepths doit etre un vecteur ou
#l'entree j est la profondeur de Tukey au j ieme point de l'echantillon.
#-derniereregion=T indique que le calcul d'une region tronquee de
#profondeur inferieur a deja ete cacule.
#-Dans le cas ou derniereregion=T, normsdern doit contenir la matrice
#contenant les vecteurs normales associe aux cotes du polygone
#spherique correspondant a la region tronquee anterieur.
#Les 4 dernieres options sont utilises par defaut par la fonction STD()
#et il n'est pas essentiel pour l'utilisateur de bien comprendre leur
#fonctionnement.
#Sorties:
#-Le graphe du polygone spherique determinant la region tronquee.
#-Une liste a 4 composantes:
#1:Les points representant les sommets de la region tronquee (polygone
#spherique).
#2:Les faces du polygone spherique en question.
#3:Un scalaire valant 0,1 ou 2 selon le cas ou l'espace tronque est
#vide, un polygone spherique ou la sphere au complet respectivement.
#4:La matrice des vecteurs normales aux faces du polygone spherique.
#polygone spherique ou la sphere au complet respectivement.
{
A=espacetronque(P,alpha,demicalcul=demicalcul,pointshemi=pointshemi,
tukdepthscalcul=tukdepthscalcul,tukdepths=tukdepths,
derniereregion=derniereregion,normsdern=normsdern)
if(A[[3]]==1)
{
pointsregion=A[[1]]
if(open==FALSE)
{
rgl.open()
}
if(rgl==FALSE)
{
rgl.spheres(0,0,0,1)
}
if(is.matrix(pointsregion)==0)
{
rgl.points(pointsregion[1],pointsregion[2],pointsregion[3],size=3,
col="red")
}
else
{
faces=A[[2]]
if(tracepoints==TRUE)
{
rgl.points(P[,1],P[,2],P[,3],size=sizepoints,col=colpoints)
}
tracepolyspher(pointsregion,faces,colarc=colarc)
}
}
return(list(A[[1]],A[[2]],A[[3]],A[[4]]))
}
sphericaldevide=function(sommets,faces)
#Cette fonctions decoupe un polygone spherique en triangles
#spherique ayant tous des angles inferieur a pi.
#Entrees:
#1: sommets est une matrice n par 3 decrivant les n sommets du polygone.
#2: faces est une matrice n par 2 decrivant les n cotes du polygone
#spherique a l'aides d'indices referants aux entrees de la matrice de
#sommets.
#Sorties: Un tableau tridimentionnel de dimension m par 3 par 3
#decrivant m triangles spheriques subdivisant le polygone spherique.
{
npoints=nrow(sommets)
nfaces=nrow(faces)
pointsopp=c(1:npoints)[-1]
trianglespheric=array(NA,dim=c(npoints-2,3,3))
compt=0
for(i in 1:nfaces)
{
if(faces[i,1]!=1 && faces[i,2]!=1)
{
compt=compt+1
trianglespheric[compt,,]=rbind(sommets[1,],sommets[faces[i,],])
}
}
k=1
while(k <=compt)
{
tri=trianglespheric[k,,]
A=tri[1,]
B=tri[2,]
C=tri[3,]
alpha=acos(B%*%C)
beta=acos(A%*%C)
gamma=acos(A%*%B)
max=max(c(alpha,beta,gamma))
if(max>=pi/2)
{
if(alpha>=pi/2 && round(alpha==max,15))
{
milieu=((B+C)/2)
milieu=milieu/sqrt(milieu%*%milieu)
tri1=rbind(A,B,milieu)
tri2=rbind(A,C,milieu)
}
if(beta>=pi/2 && round(beta==max,15))
{
milieu=((A+C)/2)
milieu=milieu/sqrt(milieu%*%milieu)
tri1=rbind(B,C,milieu)
tri2=rbind(B,A,milieu)
}
if(gamma>=pi/2 && round(gamma==max,15))
{
milieu=((A+B)/2)
milieu=milieu/sqrt(milieu%*%milieu)
tri1=rbind(C,A,milieu)
tri2=rbind(C,B,milieu)
}
if(k==1)
{
if(compt!=1)
{
trianglespheric=abind(array(as.vector(tri1),dim=c(1,3,3)),
array(as.vector(tri2),dim=c(1,3,3)),
trianglespheric[2:compt,,],along=1)
}
else
{
trianglespheric=abind(array(as.vector(tri1),dim=c(1,3,3)),
array(as.vector(tri2),dim=c(1,3,3)),along=1)
}
}
if(k==compt && compt!=1)
{
trianglespheric=abind(trianglespheric[1:(compt-1),,],
array(as.vector(tri1),dim=c(1,3,3)),array(as.vector(tri2),
dim=c(1,3,3)),along=1)
}
if(k>1 && k<compt)
{
trianglespheric=abind(trianglespheric[1:(k-1),,],
array(as.vector(tri1),dim=c(1,3,3)),array(as.vector(tri2),
dim=c(1,3,3)),trianglespheric[(k+1):compt,,],along=1)
}
compt=compt+1
}
else
{
k=k+1
}
}
return(trianglespheric)
}
airtrispher=function(sommets)
#Cette fonction calcule l'aire d'un triangle spherique
#Entree: Une matrice 3 par 3 dont chaque ligne est un sommet du
#triangle spherique.
#Sortie: L'aire du triangle spherique.
{
angles=rep(NA,3)
arcs=rep(NA,3)
for(i in 1:3)
{
pointsopp=sommets[-i,]
v1=pointsopp[1,]
v2=pointsopp[2,]
theta=acos((v2%*%v1)/(sqrt(v1%*%v1)*sqrt(v2%*%v2)))
arcs[i]=theta
}
for(i in 1:3)
{
arcsadj=arcs[-i]
angles[i]=acos((cos(arcs[i])-cos(arcsadj[1])*
cos(arcsadj[2]))/(sin(arcsadj[1])*sin(arcsadj[2])))
}
airtrispher=sum(angles)-pi
return(airtrispher)
}
spherpolycentroide=function(sommets,faces)
#Cette fonction calcule le centre de gravite projete sur la sphere
#d'un polygone spherique
# Entrees:
#1:sommets est une matrice n par 3 decrivant les n sommets du polygone.
#2:faces est une matrice n par 2 decrivant les n cotes du polygone
#spherique a l'aides d'indices referants aux entrees de la matrice de
#sommets.
#Sortie: Un vecteur a trois entrees donnant le centre de gravite projete
#sur la sphere du polygone spherique.
{
sphertri=sphericaldevide(sommets,faces)
ntri=dim(sphertri)[1]
airtri=rep(NA,ntri)
cgtri=matrix(NA,ncol=3,nrow=ntri)
poids=rep(NA,ntri)
cgtripoids=matrix(NA,ncol=3,nrow=ntri)
for(i in 1:ntri)
{
airtri[i]=airtrispher(sphertri[i,,])
cgtri[i,]=cmsphertri(sphertri[i,,])
}
airpoly=sum(airtri)
for(i in 1:ntri)
{
poids[i]=airtri[i]/airpoly
cgtripoids[i,]=poids[i]*cgtri[i,]
}
cgpoly=rep(NA,3)
for(j in 1:3)
{
cgpoly[j]=sum(cgtripoids[,j])
}
cgpoly=cgpoly/(sqrt(cgpoly%*%cgpoly))
return(cgpoly)
}
cmsphertri=function(ABC)
#Cette fonction calcule le centre de gravite d'un triangle spherique.
#Les angles spherique du triangle doivent etre plus petit que pi/2 rad.
#Entrees: Une matrice 3 par 3 dont chaque ligne est un sommet du
#triangle spherique.
#Sortie: le centre de gravite du triangle spherique
{
A=ABC[1,]
B=ABC[2,]
C=ABC[3,]
AXB=prodvect2(rbind(A,B,c(0,0,0)))
AXB=AXB*(sign(AXB%*%C))
BXC=prodvect2(rbind(B,C,c(0,0,0)))
BXC=BXC*(sign(BXC%*%A))
CXA=prodvect2(rbind(C,A,c(0,0,0)))
CXA=CXA*(sign(CXA%*%B))
angles=rep(NA,3)
arcs=rep(NA,3)
for(i in 1:3)
{
pointsopp=ABC[-i,]
v1=pointsopp[1,]
v2=pointsopp[2,]
theta=acos((v2%*%v1)/(sqrt(v1%*%v1)*sqrt(v2%*%v2)))
arcs[i]=theta
}
for(i in 1:3)
{
arcsadj=arcs[-i]
angles[i]=acos((cos(arcs[i])-cos(arcsadj[1])*
cos(arcsadj[2]))/(sin(arcsadj[1])*sin(arcsadj[2])))
}
E=sum(angles)-pi
alpha=arcs[1]
beta=arcs[2]
gamma=arcs[3]
cg=(1/(2*E))*((AXB/(sqrt(AXB%*%AXB)))*gamma+
(BXC/(sqrt(BXC%*%BXC)))*alpha+(CXA/(sqrt(CXA%*%CXA)))*beta)
return(cg)
}
cmsphertri2=function(ABC)
#Cette fonction calcule le centre de gravite d'un triangle spherique.
#Les angles spherique du triangle doivent etre plus petit que pi/2 rad.
#Entrees: Une matrice 3 par 3 dont chaque ligne est un sommet du
#triangle spherique.
#Sortie: le centre de gravite du triangle spherique
{
A=ABC[1,]
B=ABC[2,]
C=ABC[3,]
angles=rep(NA,3)
arcs=rep(NA,3)
for(i in 1:3)
{
pointsopp=ABC[-i,]
v1=pointsopp[1,]
v2=pointsopp[2,]
theta=acos((v2%*%v1)/(sqrt(v1%*%v1)*sqrt(v2%*%v2)))
arcs[i]=theta
}
for(i in 1:3)
{
arcsadj=arcs[-i]
angles[i]=acos((cos(arcs[i])-cos(arcsadj[1])*
cos(arcsadj[2]))/(sin(arcsadj[1])*sin(arcsadj[2])))
}
excess=sum(angles)-pi
alpha=arcs[1]
beta=arcs[2]
gamma=arcs[3]
N1=prodvect2(rbind(C,A,c(0,0,0)))
N1=N1/sqrt(N1%*%N1)
N2=prodvect2(rbind(A,B,c(0,0,0)))
N2=N2/sqrt(N2%*%N2)
N3=prodvect2(rbind(C,B,c(0,0,0)))
N3=N3/sqrt(N3%*%N3)
E=C-((C%*%A)/(A%*%A))*A
E=E/sqrt(E%*%E)
if((E%*%C)<0)
{
E=-E
}
D=B-((B%*%A)/(A%*%A))*A
D=D/sqrt(D%*%D)
if((D%*%B)<0)
{
D=-D
}
F=C-((C%*%B)/(B%*%B))*B
F=F/sqrt(F%*%F)
if((F%*%C)<0)
{
F=-F
}
G=A-((A%*%B)/(B%*%B))*B
G=G/sqrt(G%*%G)
if((G%*%A)<0)
{
G=-G
}
H=B-((B%*%C)/(C%*%C))*C
H=H/sqrt(H%*%H)
if((H%*%B)<0)
{
H=-H
}
K=A-((A%*%C)/(C%*%C))*C
K=K/sqrt(K%*%K)
if((K%*%A)<0)
{
K=-K
}
K1=prodvect2(rbind(E,D,c(0,0,0)))
K1=K1/sqrt(K1%*%K1)
if((K1%*%A)<0)
{
K1=-K1
}
K2=prodvect2(rbind(F,G,c(0,0,0)))
K2=K2/sqrt(K2%*%K2)
if((K2%*%B)<0)
{
K2=-K2
}
K3=prodvect2(rbind(H,K,c(0,0,0)))
K3=K3/sqrt(K3%*%K3)
if((K3%*%C)<0)
{
K3=-K3
}
p1=pi/2-min(acos(K1%*%N3),pi-acos(K1%*%N3))
p2=pi/2-min(acos(K2%*%N1),pi-acos(K2%*%N1))
p3=pi/2-min(acos(K3%*%N2),pi-acos(K3%*%N2))
matrice=rbind(K1,K2,K3)
vecteur=c((1/2)*(alpha*sin(p1))/(excess),
(1/2)*(beta*sin(p2))/(excess),(1/2)*(gamma*sin(p3))/(excess))
cm=solve(matrice,vecteur)
return(cm)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.