---
title: "VS2.3 - Illustration of PCDs in One Tetrahedron"
author: "Elvan Ceyhan"
date: '`r Sys.Date()` '
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
bibliography: References.bib
vignette: >
  %\VignetteIndexEntry{VS2.3 - Illustration of PCDs in One Tetrahedron}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
.math {
  font-size: small;
}
</style>

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE,comment = "#>",fig.width=6, fig.height=4, fig.align = "center") 
```

\newcommand{\X}{\mathcal{X}}
\newcommand{\Y}{\mathcal{Y}}

First we load the `pcds` package:
```{r setup, message=FALSE, results='hide'}
library(pcds)
```

# Functions for PCDs with Vertices in One Tetrahedron {#sec:PCD-tetra}
For illustration of construction of PCDs and related quantities, 
we first focus on one Delaunay cell, which is a tetrahedron in the 3D setting.
Due to geometry invariance of PE- and CS-PCDs with vertices from uniform distribution in a tetrahedron in $\mathbb R^3$, 
most data generation and computations can be done in the *standard regular tetrahedron*, 
and for AS-PCD, one can restrict attention to the *basic tetrahedron*.
Here, 
we only consider the PE-PCD with vertices being 3D data points, 
extension of the other PCDs is left for future work.
The *standard regular tetrahedron* is $T_{reg}=T(A,B,C,D)$ 
with vertices $A=(0,0,0)$, $B=(1,0,0)$, and $C=\left(1/2,\sqrt{3}/2,0\right)$, and $D=\left(1/2,\sqrt{3}/6,\sqrt{6}/3\right)$.

For more detail on the extension of PCDs to higher dimensions,
see @ceyhan:Phd-thesis, @ceyhan:arc-density-PE, @ceyhan:arc-density-CS, and @ceyhan:masa-2007.

We first choose an arbitrary tetrahedron.
We choose the arbitrary tetrahedron $T=T(A,B,C,D)$ with vertices jittered around the vertices of the standard regular tetrahedron
(for better visualization).
```{r }
set.seed(1)
A<-c(0,0,0)+runif(3,-.25,.25);
B<-c(1,0,0)+runif(3,-.25,.25);
C<-c(1/2,sqrt(3)/2,0)+runif(3,-.25,.25);
D<-c(1/2,sqrt(3)/6,sqrt(6)/3)+runif(3,-.25,.25)
tetra<-rbind(A,B,C,D)
n<-5  #try also n<-10 or 20
```
And we generate $n=$ `r n` $\X$ points uniformly
in the tetrahedron $T$ using the function `runif.tetra` in `pcds`.
One can use the center of mass (default) or the circumcenter of the tetrahedron to construct the vertex regions.
$\X$ points are denoted as `Xp` and $\Y$ points correspond to the vertices of the tetrahedron $T$ (i.e. $\{A,B,C,D\}$).
```{r }
Xp<-runif.tetra(n,tetra)$g  #try also Xp[,1]<-Xp[,1]+1
```

Plot of the tetrahedron $T$ and the $\X$ points in it,
and we also add the vertex names using the `text3D` function from package `plot3D`.
```{r one-th, fig.cap="Scatterplot of the uniform $X$ points in the tetrahedron $T$."}
xlim<-range(tetra[,1],Xp[,1])
ylim<-range(tetra[,2],Xp[,2])
zlim<-range(tetra[,3],Xp[,3])

xr<-xlim[2]-xlim[1]
yr<-ylim[2]-ylim[1]
zr<-zlim[2]-zlim[1]

plot3D::scatter3D(Xp[,1],Xp[,2],Xp[,3], phi=0,theta=-60, bty = "g",main="Points in One Tetrahedron",
                  xlab="x", ylab="y", zlab="z", xlim=xlim+xr*c(-.05,.05), ylim=ylim+yr*c(-.05,.05),
                  zlim=zlim+zr*c(-.05,.05), pch = 20, cex = 1, ticktype = "detailed")
#add the vertices of the tetrahedron
plot3D::points3D(tetra[,1],tetra[,2],tetra[,3], add = TRUE)
A<-tetra[1,]; B<-tetra[2,]; C<-tetra[3,]; D<-tetra[4,]
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=1,lty=2)
#now we add the vertex names and annotation
txt<-tetra
xc<-txt[,1]+c(-.02,.02,-.15,.02)
yc<-txt[,2]+c(.02,.02,.02,.02)
zc<-txt[,3]+c(-.04,.02,.02,.02)
txt.str<-c("A","B","C","D")
plot3D::text3D(xc,yc,zc,txt.str, add = TRUE)
```

## Functions for Proportional Edge PCDs with Vertices in One Tetrahedron
We use the same tetrahedron $T$ and the generated data above, 
and choose the expansion parameter $r=1.5$
and center `M="CM"` to illustrate the PE proximity region construction.

The function `NPEtetra` is used for the construction of PE proximity regions taking the arguments
`p,th,r,M,rv` where

- `p`, a 3D point whose PE proximity region is to be computed,
- `r`, a positive real number which serves as the expansion parameter in PE proximity region; must be $\ge 1$.
- `th`, A $4 \times 3$ matrix with each row representing a vertex of the tetrahedron,
- `M`, the center to be used in the construction of the vertex regions in the tetrahedron, `th`,
Currently it only takes `"CC"` for circumcenter and `"CM"` for center of mass; default=`"CM"`,
- `rv`, the index of the vertex region containing the point, either `1,2,3,4` (default is `NULL`).

`NPEtetra` returns the proximity region as the the vertices of the tetrahedron proximity region.
```{r eval=F}
M<-"CM"  #try also M<-"CC"
r<-1.5
NPEtetra(Xp[1,],tetra,r)  #uses the default M="CM"
#>            [,1]      [,2]       [,3]
#> [1,] 0.72233763 0.9464243 0.06455702
#> [2,] 0.06169821 0.1514047 0.04242222
#> [3,] 1.10142305 0.0843472 0.17049892
#> [4,] 0.37498004 0.3131847 0.52897936
```

Indicator for the presence of an arc from a (data or $\X$) point to another for PE-PCDs is the function `IarcPEtetra`
which takes arguments `p1,p2,th,r,M,rv` where

- `p1`, a 3D point whose PE proximity region is constructed.
- `p2`, a 3D point. The function determines whether `p2` is inside the PE proximity region of `p1` or not.
- `th,r,M,rv` are as in `NPEtetra`.

This function returns $I(p_2 \in N_{PE}(p_1,r))$, 
that is, returns 1 if `p2` is in $N_{PE}(p_1,r)$, 0 otherwise.
One can use it for points in the data set or for arbitrary points (as if they were in the data set).
```{r eval=F}
IarcPEtetra(Xp[1,],Xp[2,],tetra,r)  #uses the default M="CM"
#> [1] 1
IarcPEtetra(Xp[2,],Xp[1,],tetra,r,M)
#> [1] 1
```

Number of arcs of the PE-PCD can be computed by the function `num.arcsPEtetra`,
which is an object of class "`NumArcs`" and takes arguments `Xp,th,r,M`
where `Xp` is the data set, and `th,r,M` are as in `NPEtetra`.
This function returns the list of

* `desc`: A description of the PCD and the output
* `num.arcs`: Number of arcs of the PE-PCD,
* `num.in.tetra`: Number of `Xp` points in the tetrahedron `tetra`,
* `ind.in.tetra`: The vector of indices of the `Xp` points that reside in the tetrahedron.
* `tess.points`: Vertices of the tetrahedron (i.e. `Yp` points)
* `vertices`: vertices of the PCD (i.e. `Xp` points)

The output is as in `num.arcsPEtri`.
```{r eval=F}
Narcs = num.arcsPEtetra(Xp,tetra,r,M)
summary(Narcs)
#> Call:
#> num.arcsPEtetra(Xp = Xp, th = tetra, r = r, M = M)
#>
#> Description of the output:
#> Number of Arcs of the PE-PCD with vertices Xp and Quantities Related to the Support Tetrahedron
#>
#> Number of data (Xp) points in the tetrahedron =  5
#> Number of arcs in the digraph =  10
#>
#> Indices of data points in the tetrahedron:
#> 1 2 3 4 5 
#> 
#plot(Narcs) #gives error
```

The arc density of the PE-PCD can be computed by the function `PEarc.dens.tetra`.
It takes the arguments `Xp,th,r,M,th.cor` where `Xp,th,r,M` are as in `num.arcsPEtetra`
and `th.cor` is a logical argument for computing the arc density for only the points inside the tetrahedron,
`th`. (default is `th.cor=FALSE`), i.e., if `th.cor=TRUE` only the induced digraph
with the vertices inside `th` are considered in the computation of arc density.
Arc density can be adjusted (or corrected) for the points outside the triangle by using the option `th.cor=TRUE`.
```{r eval=F }
PEarc.dens.tetra(Xp,tetra,r,M)
#> [1] 0.5
```

The incidence matrix of the PE-PCD for the one tetrahedron case can be found by `inci.matPEtetra`, 
using the `inci.matPEtetra(Xp,tetra,r,M)` command.
It has the same arguments as `num.arcsPEtetra` function.

<!-- Plots of the PE proximity regions, the first is for all points,  -->
<!-- the second is for one of the $\X$ points only (for better visualization). -->
<!-- ```{r 3DPEPR1, fig.cap="PE proximity regions for the 10 $X$ points in the tetrahedron $T$ used above."} -->
<!-- plotPEregs.tetra(Xp,tetra,r) #uses the default M="CM" -->
<!-- ``` -->

Plots of the PE proximity region for one of the $\X$ points only (for better visualization).
```{r 3DPEPR2, fig.cap="PE proximity regions for one of the $X$ points in the tetrahedron $T$ for better visualization."}
plotPEregs.tetra(Xp[1,],tetra,r=1.5)  #uses the default M="CM"
```

## Functions for Construction and Simulation of Proportional Edge PCDs in the Standard Regular Tetrahedron {#sec:PE-PCD-std-tetra}
We first define the standard regular tetrahedron $T_{reg}$ as in Section \@ref(sec:PCD-tetra).
```{r eval=F}
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<-10  #try also n<-20
```
And we generate $\X$ points of size $n=$ `r n` using the function `runif.std.tetra` in the `pcds` package, 
and use either center of mass (default) or the circumcenter of $T_{reg}$ to construct vertex regions.
```{r eval=F}
Xp<-runif.std.tetra(n)$g
```

Plot of the standard regular tetrahedron $T_{reg}$ and the $\X$ points in it are provided
using the below code.
We also add the vertex names using the `text3D` function from package `plot3D`.
```{r std-reg-th, eval=F, fig.cap="Scatterplot of 10 $X$ points in the standard regular tetrahedron $T_{reg}$."}
xlim<-range(tetra[,1],Xp[,1])
ylim<-range(tetra[,2],Xp[,2])
zlim<-range(tetra[,3],Xp[,3])

xr<-xlim[2]-xlim[1]
yr<-ylim[2]-ylim[1]
zr<-zlim[2]-zlim[1]

plot3D::scatter3D(Xp[,1],Xp[,2],Xp[,3], phi=0,theta=-60, bty = "g",main="Points in the Standard Regular Tetrahedron",
                  xlab="x", ylab="y", zlab="z", xlim=xlim+xr*c(-.05,.05), ylim=ylim+yr*c(-.05,.05),
                  zlim=zlim+zr*c(-.05,.05), pch = 20, cex = 1, ticktype = "detailed")
#add the vertices of the tetrahedron
plot3D::points3D(tetra[,1],tetra[,2],tetra[,3], add = TRUE)
A<-tetra[1,]; B<-tetra[2,]; C<-tetra[3,]; D<-tetra[4,]
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=1,lty=2)

#now we add the vertex names and annotation
txt<-tetra
xc<-txt[,1]+c(-.01,.02,-.12,.02)
yc<-txt[,2]+c(.02,.02,.02,-.01)
zc<-txt[,3]+c(-.04,.02,.02,.02)
txt.str<-c("A","B","C","D")
plot3D::text3D(xc,yc,zc,txt.str, add = TRUE)
```

The function `NPEstd.tetra` is used for the construction of PE proximity regions for points in $T_{reg}$
taking arguments `p,r,rv` which are same as in `NPEtetra`.
This function returns the proximity region as the the vertices of the tetrahedron proximity region.
Center of mass is equivalent to circumcenter for $T_{reg}$, so no argument is specified for the center.
```{r eval=F}
r<-1.5
NPEstd.tetra(Xp[1,],r)
#>           [,1]      [,2]      [,3]
#> [1,] 0.5000000 0.8660254 0.0000000
#> [2,] 0.1618881 0.2803983 0.0000000
#> [3,] 0.5000000 0.4756074 0.5521345
#> [4,] 0.8381119 0.2803983 0.0000000
NPEstd.tetra(Xp[5,],r)
#>            [,1]      [,2]      [,3]
#> [1,] 1.00000000 0.0000000 0.0000000
#> [2,] 0.52499658 0.8227301 0.0000000
#> [3,] 0.52499658 0.2742434 0.7756773
#> [4,] 0.04999316 0.0000000 0.0000000
```

Indicator for the presence of an arc from a (data or $\X$) point to another for PE-PCDs is the function `IarcPEstd.tetra`
which takes arguments `p1,p2,r,rv` and
returns $I(p_2 \in N_{PE}(p_1,r))$ which are same as in `IarcPEtetra`.
One can use it for points in the data set or for arbitrary points (as if they were in the data set).
```{r eval=F}
IarcPEstd.tetra(Xp[1,],Xp[3,],r)  #uses the default M="CM"
#> [1] 1
```

<!-- Plots of the PE proximity regions, the first is for all points, the second is for one of the $\X$ points only (for better visualization). -->
Plots of the PE proximity regions for all points can be obtained with the below code:
```{r 3DPEPRsth1, eval=F, fig.cap="PE proximity regions for the 10 uniform $X$ points in $T_{reg}$ used above."}
plotPEregs.std.tetra(Xp,r)
```

<!-- ```{r 3DPEPRsth2, fig.cap="PE proximity regions for one of the $X$ points in $T_{reg}$ for better visualization."} -->
<!-- plotPEregs.std.tetra(Xp[1,],r) -->
<!-- ``` -->

## Auxiliary Functions to Define Proximity Regions for Points in a Tetrahedron
We only use circumcenter (CC) or center of mass (CM) to construct vertex regions in a tetrahedron, 
and here we illustrate finding CC (as CM is just the componentwise arithmetic average of the vertices)
of a tetrahedron.
The function `circumcenter.tetra` takes the sole argument `th` for a $4 \times 3$ matrix with each row representing a vertex of the tetrahedron and returns the circumcenter of a general tetrahedron.
We use a tetrahedron $T$ whose vertices are jittered around the vertices of a standard regular tetrahedron for better visualization
in the plots.
```{r eval=F}
set.seed(123)
A<-c(0,0,0)+runif(3,-.2,.2);
B<-c(1,0,0)+runif(3,-.2,.2);
C<-c(1/2,sqrt(3)/2,0)+runif(3,-.2,.2);
D<-c(1/2,sqrt(3)/6,sqrt(6)/3)+runif(3,-.2,.2);
tetra<-rbind(A,B,C,D)

CC<-circumcenter.tetra(tetra)
CC
#> [1] 0.5516851 0.3386671 0.1212977
```

<!-- We illustrate the CC in a tetrahedron in the following plot.  -->
<!-- We annotate the vertices and the CC. -->
<!-- See `example(circumcenter.tetra)` for the code to generate this figure. -->
<!-- ```{r 3DCC, fig.cap="Circumcenter of the tetrahedron $T$."} -->
<!-- Xlim<-range(tetra[,1],CC[1]) -->
<!-- Ylim<-range(tetra[,2],CC[2]) -->
<!-- Zlim<-range(tetra[,3],CC[3]) -->
<!-- xd<-Xlim[2]-Xlim[1] -->
<!-- yd<-Ylim[2]-Ylim[1] -->
<!-- zd<-Zlim[2]-Zlim[1] -->

<!-- plot3D::scatter3D(tetra[,1],tetra[,2],tetra[,3], phi =0,theta=40, bty = "g", -->
<!--                   main="Illustration of the Circumcenter\n in a Tetrahedron", -->
<!--                   xlim=Xlim+xd*c(-.05,.05),ylim=Ylim+yd*c(-.05,.05), zlim=Zlim+zd*c(-.05,.05), -->
<!--                   pch = 20, cex = 1, ticktype = "detailed") -->
<!-- #add the vertices of the tetrahedron -->
<!-- plot3D::points3D(CC[1],CC[2],CC[3], add=TRUE) -->
<!-- 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) -->

<!-- plot3D::text3D(tetra[,1],tetra[,2],tetra[,3], labels=c("A","B","C","D"), 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(CC,6),byrow = TRUE,ncol=3) -->
<!--   plot3D::segments3D(L[,1], L[,2], L[,3], R[,1], R[,2],R[,3], add=TRUE,lty=2) -->

<!-- plot3D::text3D(CC[1],CC[2],CC[3], labels="CC", add=TRUE) -->
<!-- ``` -->

We can also illustrate the CC in a tetrahedron in a plot with vertex and CC annotation. 
Type `? circumcenter.tetra` for the code to generate the corresponding figure.

The function `rel.vert.tetraCC` takes arguments `p,th` and
returns the index of the $CC$-vertex region in a tetrahedron `th` where the point `p` resides.
```{r eval=F}
n<-10  #try also n<-20
Xp<-runif.tetra(n,tetra)$g
rel.vert.tetraCC(Xp[1,],tetra)
#> $rv
#> [1] 2
#> 
#> $tetra
#>                 [,1]      [,2]        [,3]
#> vertex 1 -0.08496899 0.1153221 -0.03640923
#> vertex 2  1.15320696 0.1761869 -0.18177740
#> vertex 3  0.51124220 1.0229930  0.02057401
#> vertex 4  0.48264589 0.4714085  0.79783024

Rv<-vector()
for (i in 1:n)
  Rv<-c(Rv,rel.vert.tetraCC(Xp[i,],tetra)$rv)
Rv
#>  [1] 2 2 1 3 2 1 2 3 2 1
```
We also generate $n=$ `r n` $\X$ points uniformly in the tetrahedron and 
find the indices of the vertex regions they reside in.

We illustrate the CC-vertex regions using the following code. 
We annotate the vertices and CC and provide the scatterplot of these points labeled according to the vertex region they reside in. 
Type also `? rel.vert.tetraCC` for the code to generate this figure.
```{r 3DCCVR, eval=F, fig.cap="CC-Vertex regions in the tetrahedron $T=(A,B,C,D)$."}
CC<-circumcenter.tetra(tetra)
CC

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

plot3D::scatter3D(tetra[,1],tetra[,2],tetra[,3], phi =0,theta=40, bty = "g",
                  main="Scatterplot of data points with CC-vertex regions",
                  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(CC[1],CC[2],CC[3], labels=c("CC"), 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(CC,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,CC,B,C,D)
L<-matrix(rep(F1,4),ncol=3,byrow=TRUE); R<-rbind(D4,D5,D6,CC)
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,CC,A,C,D)
L<-matrix(rep(F2,4),ncol=3,byrow=TRUE); R<-rbind(D2,D3,D6,CC)
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,CC,A,B,D)
L<-matrix(rep(F3,4),ncol=3,byrow=TRUE); R<-rbind(D3,D5,D6,CC)
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,CC,A,B,C)
L<-matrix(rep(F4,4),ncol=3,byrow=TRUE); R<-rbind(D1,D2,D4,CC)
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)
```

The function `rel.vert.tetraCM` takes the same arguments as `rel.vert.tetraCC` and 
returns the index of the $CM$-vertex region in a tetrahedron `th` where the point `p` resides.
We use standard regular tetrahedron in this example.
```{r eval=F}
#The index of the $CM$-vertex region in a tetrahedron that contains a point
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<-10  #try also n<-20
Xp<-runif.std.tetra(n)$g
rel.vert.tetraCM(Xp[1,],tetra)
#> $rv
#> [1] 4
#> 
#> $tetra
#>          [,1]      [,2]      [,3]
#> vertex 1  0.0 0.0000000 0.0000000
#> vertex 2  1.0 0.0000000 0.0000000
#> vertex 3  0.5 0.8660254 0.0000000
#> vertex 4  0.5 0.2886751 0.8164966

Rv<-vector()
for (i in 1:n)
  Rv<-c(Rv, rel.vert.tetraCM(Xp[i,],tetra)$rv )
Rv
#>  [1] 4 4 3 4 4 1 1 3 3 2
```
We also generate $n=$ `r n` $\X$ points in the tetrahedron and find the indices of the vertex regions they reside in.

We can also illustrate the CC in a tetrahedron in a plot with vertex and CC annotation. 
Type `? rel.vert.tetraCM` for the code to generate the corresponding figure.

<!-- We illustrate the CM-vertex regions in the following plot. -->
<!-- We annotate the vertices and CM and provide the scatterplot of these points labeled according to the vertex region they reside in.  -->
<!-- See `example(rel.vert.tetraCM)` for the code to generate this figure. -->
<!-- ```{r 3DCMVR, fig.cap="CM-Vertex regions in the tetrahedron $T=(A,B,C,D)$."} -->
<!-- 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", -->
<!--                   main="Scatterplot of data points \n and CM-vertex regions", -->
<!--                   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) -->
<!-- ``` -->

# References