Curvature decomposition on surface

Curvature decomposition https://rpubs.com/lizhi1800/Surface_Normal_Tangent_Plane

There is a surface \(S\) and curve \(C\) the on surface. At point \(P\), the surface has a unit normal vector \(N\) and tangent plane \(\bar r_u,\bar r_v\). At point \(P\), the curve has a unit tangent vector \(t\) and unit normal vector \(n\). The curve’s curvature vector \(\bar k\) is in direction of \(n\), \(1/\bar k\) is a radius of curvature circle for curve.

Normal curvature \(\bar k_n\) is projection of \(\bar k\) onto plane containing vector \(N\). Geodesic curvature \(\bar k_g\) is projection vector \(\bar k\) onto \(\bar r_u, \bar r_v\) tangent plane.

The 3 vectors \(\bar k, \bar k_n, \bar k_g\) has relation \[\bar k=\frac{dt}{ds}= \bar k_n+ \bar k_g\] The 3 magnitudes \(k, k_n, k_g\) has relation \[k=||\bar k||,k_n=||\bar k_n||, k_g=||\bar k_g||\] \[k^2= k_n^2+ k_g^2\] So we say the curvature \(\bar k\) has geodesic and normal components.

We can easily verify \[\bar k=kn,\bar k_n=k_nN,\bar k_g=?\]

Normal Curvature

https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node29.html#eqn:surf_kappa_n_lambda

The fundamental coefficients are 6 quantities, note I use \(N\) for both unit normal and a fundamental coefficient. \[E=\bar r_u \cdot \bar r_u,F=\bar r_u \cdot \bar r_v,G=\bar r_v \cdot \bar r_v\] \[L=\bar r_{uu}\cdot N,M=\bar r_{uv}\cdot N,N=\bar r_{vv}\cdot N\] https://tutorial.math.lamar.edu/Classes/CalcII/DotProduct.aspx

So the normal curvature \(k_n\) can be computed as a projection onto vector \(N\) by dot product, where \(\theta\) is direction of tangent line \(t\) to \(C\) at \(P\). \[k_n(\theta)=\frac{dt}{ds}\cdot N=\frac{L\cos^2(\theta)+2M\cos(\theta)\sin(\theta)+N\sin^2(\theta)}{E\cos^2(\theta)+2F\cos(\theta)\sin(\theta)+G\sin^2(\theta)}\]

For example, cylinder \(\bar r=\{x=\cos(u),y=\sin(u),z=v\}\)

library(Deriv)
## Warning: package 'Deriv' was built under R version 4.0.5
#rbar
x="cos(u)"
y="sin(u)"
z="v"

r=c(x,y,z)

First order partial derivatives, \(\bar r_u\), \(\bar r_v\)

# 1st order partial
drdu=sapply(r,function(m) Deriv(m,"u"))
drdv=sapply(r,function(m) Deriv(m,"v"))

Second order partial derivatives

# 2nd order partial
drdudu=sapply(drdu,function(m) Deriv(m,"u"))
drdudv=sapply(drdu,function(m) Deriv(m,"v"))
drdvdv=sapply(drdv,function(m) Deriv(m,"u"))
drdvdu=sapply(drdv,function(m) Deriv(m,"v"))

In this case, tangent vectors \(\bar r_u\) \(\bar r_v\) are unit vectors, dot product we have \[E=1,F=0,G=1\] and \[E\cos^2(\theta)+2F\cos(\theta)\sin(\theta)+G\sin^2(\theta)=1\]

dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}
dot_prod(drdu,drdu)
## [1] "(cos(u)^2 + sin(u)^2)"
dot_prod(drdu,drdv)
## [1] "(0)"
dot_prod(drdv,drdv)
## [1] "(1)"
E=1
F=0
G=1

From our previous functions, we compute the principle normal vector as \(N=\bar r_u \times \bar r_v=\{\cos(\theta),\sin(\theta),0\}\).

N=c("cos(u)","sin(u)","0")

And the Fundamental Coefficients \(L ,M,N\)

l = dot_prod(drdudu,N)
m = dot_prod(drdudv,N)
n = dot_prod(drdvdv,N)

Quadratic form product

quadratic_prod <-function(v,m){
 x=expand.grid(v,v)
 res=matrix(paste("(",as.vector(m),")*(",x[,1],")*(",x[,2],")"),nrow=length(v))
 return(res)
}

We use \(\cos(u)^2 + \sin(u)^2=1\), therefore \(k_n(\theta)= -\cos^2(\theta)\).

numerator=paste0(l,"*cos(theta)^2+2*",m,"*cos(theta)*sin(theta)+",n,"*sin(theta)^2")
kn=Simplify(numerator)
kn
## [1] "-(cos(theta)^2 * (cos(u)^2 + sin(u)^2))"

Euler’s theorem

#https://weiqinggu.github.io/Math142/schedule.html #https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node34.html

Euler’s theorem of curvature provides more properties for normal curvature at point \(P\).

It says, at point p, normal curvature in all directions can be positive, negative or zero. The surface normal vector \(N\)’s direction is the positive direction.

At point p, normal curvature in all directions have minimum \(k_1\) and maximum \(k_2\). These 2 normal curvatures are perpendicular to each other.

For a \(k_n\), let \(\phi\) be angle between \(k_n\) and \(k_1\), we can compute \[k_n=k_1\cos^2(\phi)+k_2\sin^2(\phi)\] \(\theta\) and \(\phi\) are angles to rotate normal plane which containing vector \(N\).

Dupin’s indicatrix

https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node34.html

\[x=\cos(t)/\sqrt{|k_n|}=\cos(t)/\sqrt{k_1\cos^2(\phi)+k_2\sin^2(\phi)}\] \[y=\sin(t)/\sqrt{|k_n|}=\sin(t)/\sqrt{k_1\cos^2(\phi)+k_2\sin^2(\phi)}\] In case 2, there needs to be a sign change. There is an angle \(q\), when \(t\) reaches \(q\) signs of \(k_1,k_2\) need to be changed. \[q=\arctan(\sqrt{|k_1/k_2|})\] Sign=1, if \(t \in (q , \pi-q)\) or \(( \pi+q , 2\pi-q)\), otherwise, sign=-1.

Case 1. \(k_1\) and \(k_2\) have same sign.

t <- seq(0, 2*pi, 0.1) 

k1= 1/15
k2=1/5
x=cos(t)/sqrt(k1*cos(t)^2+k2*sin(t)^2)
y=sin(t)/sqrt(k1*cos(t)^2+k2*sin(t)^2)

plot(x,y,pch=19, col='blue')

Case 2. \(k_1\) and \(k_2\) have opposite signs.

t <- seq(0, 2*pi, 0.1) 

k1=-1/15
k2= 1/5

ang=atan(sqrt(abs(k1/k2)))
sign=ifelse((t>ang & t<(pi-ang))|(t>pi+ang & t<(2*pi-ang)),1,-1)

x=cos(t)/sqrt(sign*k1*cos(t)^2+sign*k2*sin(t)^2)
y=sin(t)/sqrt(sign*k1*cos(t)^2+sign*k2*sin(t)^2)

plot(x,y,pch=19, col='blue')

Case 3. One of the numbers k1 and k2 is equal to zero.

t <- seq(0, 2*pi, 0.1) 

k1= 0
k2=1/5
x=cos(t)/sqrt(k1*cos(t)^2+k2*sin(t)^2)
y=sin(t)/sqrt(k1*cos(t)^2+k2*sin(t)^2)

plot(x,y,pch=19, col='blue')

Meusnier’s theorem

https://www.math.brown.edu/tbanchof/balt/ma106/lab8.html?dtext86.html

A sphere coincides surface normal vector \(N\) at point \(P\) on surface \(S\) with radius \(r_n=1/k_n\), which has a curve tangent vector \(t\) at \(P\), is called a Meusnier’s sphere.

We need to know the curve tangent \(t\), because as rotating \(t\), the \(k_n\) would vary between \(k_1\), \(k_2\).

So we cut the surface \(S\) perpendicularly to vector \(t\), the sphere becomes circle, two osculating planes along the vector \(t\) become lines \(N\) and \(n\).

Meusnier’s theorem

Vector \(n\) is the curve’s normal at \(P\), it has a plane with vector \(t\), slices the sphere and creates a circle \(PQ\), but from our way of looking \(PQ\) is a line in our picture.

Circle \(PQ\) has radius \(r_L\), we can show \(r_L=1/k\). \[r_L=\cos(\theta)r_n=\cos(\theta)\frac{1}{k_n}=\frac{1}{k}\] So \(PQ\) is the osculating circle at \(P\) for the curve.

Principal directions

https://math.stackexchange.com/questions/2571902/find-principal-curvature-and-gaussian-curvature

Consider cylinder, we do 1st and 2nd degree partial derivative on surface point \[\bar r=\{\cos(v), \sin(v), u\}\]

library(Deriv)
rbar=c("cos(v)", "sin(v)", "u")

# 1st order partial
drdv=sapply(rbar,function(m) Deriv(m,"v"))
drdu=sapply(rbar,function(m) Deriv(m,"u"))
 
# 2nd order partial
drdvdv=sapply(drdv,function(m) Deriv(m,"v"))
drdvdu=sapply(drdv,function(m) Deriv(m,"u"))
drdudu=sapply(drdu,function(m) Deriv(m,"u"))
drdudv=sapply(drdu,function(m) Deriv(m,"v"))

Form 1 \[Form_1=\left| {\begin{array} *{E}&{F}\\ {F}&{G} \end{array}} \right|=\left| {\begin{array} *{1}&{0}\\ {0}&{1} \end{array}} \right|\] Form 2 \[Form_2=\left| {\begin{array} *{L}&{M}\\ {M}&{N} \end{array}} \right|=\left| {\begin{array} *{1}&{0}\\ {0}&{0} \end{array}} \right|\]

dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}

#form1
dot_prod(drdu,drdu)
## [1] "(1)"
dot_prod(drdu,drdv)
## [1] "(0)"
dot_prod(drdv,drdv)
## [1] "(cos(v)^2 + sin(v)^2)"
E=1
F=0
G=1

cross_prod<-function(x,y){
a1=paste0("(",x[2],")*(",y[3],")-(",y[2],")*(",x[3],")")
a2=paste0("(",x[3],")*(",y[1],")-(",y[3],")*(",x[1],")")
a3=paste0("(",x[1],")*(",y[2],")-(",y[1],")*(",x[2],")")

return(c(Simplify(a1),Simplify(a2),Simplify(a3)))
}

N=cross_prod(drdu,drdv)

#form2
dot_prod(drdvdv,N)
## [1] "(cos(v)^2 + sin(v)^2)"
dot_prod(drdvdu,N)
## [1] "(0)"
dot_prod(drdudu,N)
## [1] "(0)"
l = 1
m = 0
n = 0

The principal curvatures are the roots of \[det(Form_2-kForm_1)=0\] Hence \[det(\left| {\begin{array} *{1}&{0}\\ {0}&{0} \end{array}} \right|-k\left| {\begin{array} *{1}&{0}\\ {0}&{1} \end{array}} \right|)=0\] \[det(\left| {\begin{array} *{1-k}&{0}\\ {0}&{0-k} \end{array}} \right|)=0\] \[(1-k)(0-k)=0\] Hence, \(k_1=0\) and \(k_2=1\)

https://math.hmc.edu/calculus/hmc-mathematics-calculus-online-tutorials/linear-algebra/eigenvalues-and-eigenvectors/

For the eigenvectors \(T_i\) we solve \[(Form_2-k_iForm_1)T_i=0\] When \(k_1=0\) \[(\left| {\begin{array} *{1}&{0}\\ {0}&{0} \end{array}} \right|-k_1\left| {\begin{array} *{1}&{0}\\ {0}&{1} \end{array}} \right|)T_1=\left| {\begin{array} *{0}\\ {0} \end{array}} \right|\] Let \(T_1=\{e,n\}\), from equation in this case we see \(e=0\) and \(n\) can be any real number that vector direction would not change. So we choose \(n=1\) making \(T_1=\{0,1\}\) an unit vector.

When \(k_2=1\) \[(\left| {\begin{array} *{1}&{0}\\ {0}&{0} \end{array}} \right|-k_2\left| {\begin{array} *{1}&{0}\\ {0}&{1} \end{array}} \right|)T_2=\left| {\begin{array} *{0}\\ {0} \end{array}} \right|\] Let \(T_2=\{e,n\}\), from equation in this case we see \(n=0\) and \(e\) can be any real number that vector direction would not change. So we choose \(e=1\) making \(T_2=\{1,0\}\) an unit vector.

k1=0
k2=1

#unit eigenvector
T1=c(0,1)
T2=c(1,0)

The principal directions are \[t_i=\{\bar r_u,\bar r_v\}T_i\] I use \(\{,,...\}\) to denote c() or cbind() in R.

Principal vector \(t_1\) \[t_1=\{0,0,1\}\] Principal vector \(t_2\) \[t_2=\{−\sin (v), \cos (v), 0\}\]

tangent_plane=cbind(drdv,drdu)

print("t1")
## [1] "t1"
dot_prod(tangent_plane[1,],T1)
## [1] "(0)"
dot_prod(tangent_plane[2,],T1)
## [1] "(0)"
dot_prod(tangent_plane[3,],T1)
## [1] "(1)"
print("t2")
## [1] "t2"
dot_prod(tangent_plane[1,],T2)
## [1] "(-sin(v))"
dot_prod(tangent_plane[2,],T2)
## [1] "(cos(v))"
dot_prod(tangent_plane[3,],T2)
## [1] "(0)"

Mean curvature

\[K_{average}=\frac{k_1+k_2}{2}\]

Total curvature / Gaussian curvature

\[K=k_1k_2\]