離散曲面上のガウス曲率 - 驚異のアニヲタ社会復帰の予備

3次元メッシュ上の点のガウス曲率を求める。
これとかこれ (PDF)とかに載っている。
　
あるvertex $v$ の周辺に $v_i$ の点がある。これはよく1-ring と書かれている。
結論から書くと、青で示された面積 $A_{mixed}$ , $\{v, v_i, v_{i+1}\}$ で囲まれる三角形の $v$ の角度 $\theta_{i}$ を用いて
$G=\frac{2\pi-\sum_{i}\theta_{i}}{A_{mixed}}$
で定められる。曲面についてオイラー数とかガウス・ボネの定理とか使うと求まる。
$G=0$ のとき、 $\sum_{i}\theta_{i}=2\pi$ である。これは、もともと基本曲率 $k_1,k_2$ によってガウス曲率が $G=\frac{k_1+k_2}{2}$ と書けることにより、 $k_1$ , $k_2$ のどちらかが0である。このとき、1-ringの $v$ は平面となるのに一致してそう。
　
具体的な計算は、 $\{v, v_i, v_{i+1}\}$ で囲まれる三角形について、 $\theta_i$ が鋭角(non-obtuse)ならば外心点 $c_i$ を、 $\theta_i$ が鈍角(obtuse)ならばedge $v_i, v_{i+1}$ の中点 $c_i$ とする。 $c_i$ がすべて外心ならばボロノイ分割となるが、中点が混じればmixed cell という。
このとき、青色の領域の面積は
$A_{mixed}=\sum_{v_i\in N(v)}\{A(c_i,v,\frac{v+v_i}{2})+A(c_{i+1},v,\frac{v+v_i}{2})\}$
図はPDFより借用。

mesh オブジェクトの四元数座標、三角形index の入った三角面情報を入力として、全頂点についてガウス曲率を計算する。

ガウス曲率
口の部分や足先など出っ張っている部分は大きい(赤)ガウス曲率、背中や足の裏は0に近い(黒)、ももの付け根や鼻のシワなど凹んでいるところは小さい(緑)ガウス曲率になった。

　
平均曲率も計算しておいた。ガウス曲率と似たような、凸は正、凹は負、な感じだが、足の裏は負である。

　
平均曲率は三角形平面についても計算してあるので、面についてプロットした。ももの付け根や鼻のシワなど、境界がダイナミックに変化しそうな部分の三角形ははっきり色が分かれる。が、ももの付け根(凹)と足先(凸)が同じ色。
たぶん、ガウス曲率の三角形面ならば、凹凸は一致しそう。三角形面のガウス曲率は平均曲率と同様に、頂点の曲率から求まりそうなのだが、定義式がいずこ?

　
plot3dの技術的な話として、三角形面をshade3d でプロットするときは、面の情報は3つの頂点で決まるので、色情報が3倍必要である。
そして遅い。R でやると遅い。8スレッド並列化しているけど遅い。
C++ で書き換え中。
また、 $c_i$ の選び方は円順列にしないと $A_{mixed}$ の領域がかけないが、 $v_i$ を順番にたどるのがハミルトン回路になるはずなのにigrpah にないのこれ?

# cumpute Gaussian curvature of j_th vertex
# input: index of vertices of quanterion vertices object
gaussian_curvature <- function(j, vertices, faces.v){
  #for(j in seq(n_vertices)){
    plane <- t(faces.v[, apply(faces.v==j, 2, any)])
    # shuffle the order of summation of triangle areas
    p <- p1 <- 1 # start
    s <- c(j, setdiff(plane[p,], j)[1]) # 2nd vertex
    p1 <- c(p1, setdiff(which(apply(plane == tail(s, 1), 1, any)), p))
    while(p < nrow(plane) - 1){
      s <- c(s, setdiff(plane[tail(p1, 1),], s))
      p1 <- c(p1, which(sapply(apply(plane, 1, setdiff, s), length) == 1)[2])
      p <- p + 1
    }
    plane <- plane[p1,] # ordered sequentially

    circum_mat <- matrix(1, 3, 3)
    diag(circum_mat) <- -1
    theta_j <- rep(0, nrow(plane))
    circumcenter <- matrix(0, nrow(plane), 3)
    for(p in seq(nrow(plane))){
      hogemat <- t(as.matrix(vertices[plane[p,]])[-1,])
      rname <- plane[p,]
      e12 <- sweep(hogemat[!(rname %in% j),], 2, hogemat[ (rname %in% j),], "-")
      theta_j[p] <- acos( e12[1,] %*% e12[2,]/prod(sqrt(rowSums(e12^2))) )
      if(theta_j[p] < pi/2){ # non-obtuse angle compute circumcenter point
        abc <- mapply(function(k) sum(apply(hogemat[c(k%%3+1, (k+1)%%3+1),], 2, diff)^2), seq(3))
        # abc <- rowSums(rbind(apply(e12, 2, diff), e12)^2)
        abc2 <- circum_mat %*% abc * abc
        circumcenter[p,] <- colSums(diag(c(abc2)) %*% hogemat)/sum(abc2)
        # circumcenter <- (abc2[1]*hogemat[1,] + abc2[2]*hogemat[2,] + abc2[3]*hogemat[3,])/sum(abc2)
      } else { # obtuse angle
        circumcenter[p,] <- colSums(hogemat[!(rname %in% j),])/2 # center point
      }
    }
    mat <- sweep(circumcenter, 2, as.numeric(vertices[j])[-1], "-") # vector from vertex j_th
    n <- nrow(circumcenter)
    Amixed <- mapply(function(k) sqrt(prod(rowSums(mat[c(k, k%%n+1),]^2)) - (mat[k,]%*%mat[k%%n+1,])^2) / 2, seq(n))
    Gcurvature <- (2*pi - sum(theta_j))/sum(Amixed)
  #}
  return(Gcurvature)
}