rel.vert.tetraCM: The index of the CM-vertex region in a tetrahedron that...

View source: R/AuxDelaunay.R

rel.vert.tetraCMR Documentation

The index of the CM-vertex region in a tetrahedron that contains a point

Description

Returns the index of the vertex whose region contains point p in a tetrahedron th=T(A,B,C,D) and vertex regions are based on the center of mass CM=(A+B+C+D)/4 of th. (see the plots in the example for illustrations).

The vertices of the tetrahedron th are labeled as 1=A, 2=B, 3=C, and 4=C also according to the row number the vertex is recorded in th.

If the point, p, is not inside th, then the function yields NA as output. The corresponding vertex region is the simplex with the vertex, CM, and midpoints of the edges adjacent to the vertex.

See also (\insertCiteceyhan:Phd-thesis,ceyhan:comp-geo-2010;textualpcds).

Usage

rel.vert.tetraCM(p, th)

Arguments

p

A 3D point for which CM-vertex region it resides in is to be determined in the tetrahedron th.

th

A 4 \times 3 matrix with each row representing a vertex of the tetrahedron.

Value

A list with two elements

rv

Index of the CM-vertex region that contains point, p in the tetrahedron th

th

The vertices of the tetrahedron, where row number corresponds to the vertex index in rv.

Author(s)

Elvan Ceyhan

References

\insertAllCited

See Also

rel.vert.tetraCC and rel.vert.triCM

Examples

## Not run: 
A<-c(0,0,0); B<-c(1,0,0); C<-c(1/2,sqrt(3)/2,0);
D<-c(1/2,sqrt(3)/6,sqrt(6)/3)
tetra<-rbind(A,B,C,D)

n<-20  #try also n<-40

Xp<-runif.std.tetra(n)$g

rel.vert.tetraCM(Xp[1,],tetra)

Rv<-vector()
for (i in 1:n)
  Rv<-c(Rv, rel.vert.tetraCM(Xp[i,],tetra)$rv )
Rv

Xlim<-range(tetra[,1],Xp[,1])
Ylim<-range(tetra[,2],Xp[,2])
Zlim<-range(tetra[,3],Xp[,3])
xd<-Xlim[2]-Xlim[1]
yd<-Ylim[2]-Ylim[1]
zd<-Zlim[2]-Zlim[1]

CM<-apply(tetra,2,mean)

plot3D::scatter3D(tetra[,1],tetra[,2],tetra[,3], phi =0,theta=40, bty = "g",
xlim=Xlim+xd*c(-.05,.05),ylim=Ylim+yd*c(-.05,.05), zlim=Zlim+zd*c(-.05,.05),
          pch = 20, cex = 1, ticktype = "detailed")
L<-rbind(A,A,A,B,B,C); R<-rbind(B,C,D,C,D,D)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3], add=TRUE,lwd=2)
#add the data points
plot3D::points3D(Xp[,1],Xp[,2],Xp[,3],pch=".",cex=3, add=TRUE)

plot3D::text3D(tetra[,1],tetra[,2],tetra[,3],
labels=c("A","B","C","D"), add=TRUE)
plot3D::text3D(CM[1],CM[2],CM[3], labels=c("CM"), add=TRUE)

D1<-(A+B)/2; D2<-(A+C)/2; D3<-(A+D)/2; D4<-(B+C)/2; D5<-(B+D)/2; D6<-(C+D)/2;
L<-rbind(D1,D2,D3,D4,D5,D6); R<-matrix(rep(CM,6),ncol=3,byrow=TRUE)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3], add=TRUE,lty = 2)

F1<-intersect.line.plane(A,CM,B,C,D)
L<-matrix(rep(F1,4),ncol=3,byrow=TRUE); R<-rbind(D4,D5,D6,CM)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3],col=2,
add=TRUE,lty = 2)

F2<-intersect.line.plane(B,CM,A,C,D)
L<-matrix(rep(F2,4),ncol=3,byrow=TRUE); R<-rbind(D2,D3,D6,CM)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3],col=3,
add=TRUE,lty = 2)

F3<-intersect.line.plane(C,CM,A,B,D)
L<-matrix(rep(F3,4),ncol=3,byrow=TRUE); R<-rbind(D3,D5,D6,CM)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3],col=4,
add=TRUE,lty = 2)

F4<-intersect.line.plane(D,CM,A,B,C)
L<-matrix(rep(F4,4),ncol=3,byrow=TRUE); R<-rbind(D1,D2,D4,CM)
plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3],col=5,
add=TRUE,lty = 2)

plot3D::text3D(Xp[,1],Xp[,2],Xp[,3], labels=factor(Rv), add=TRUE)

## End(Not run)


elvanceyhan/pcds documentation built on June 29, 2023, 8:12 a.m.