読んでる。Comput Vis Image Underst. 2013 Sep 1;117(9):1107-1118.
Ricci flow という、Willmore flow とはまた違う条件で微分幾何学的物体変換を行う。
その前に下準備。
以前やった(離散)ガウス曲率は、1-ring下の三角形の角度と、と呼ばれる傘の面積で計算していたが、これはMeyer method というらしい。
本来は、ガウス曲率は角度だけで計算できて
となる。しかし、この場合は、論文の図にもあるけれども、理想的な球体のうえにランダムに点を発生させて三角形を作った時に、ガウス曲率の濃淡ができてしまう。というわけで、疎密に合わせてガウス曲率を補正したい。
というわけでで割る、らしい。
半径 の理想的な球を作っているので、ガウス曲率は
である。
ほとんどが1の周辺に分布している。
1からかけ離れたガウス曲率をもつ頂点をプロットした。赤が高いガウス曲率、青が低い。
ズームしてみると、すごい鋭角な三角形をもつ頂点が、低いガウス曲率になっている様子。
計算はあっているようである。
理想的な球を作り、点を配置し、三角形を作る。
library(rgl) library(onion) my.sphere.tri.mesh <- function(n.psi=30){ thetas <- list() psis <- seq(from=-pi/2,to=pi/2,length=n.psi) #psis <- sort(runif(n.psi, min=-pi/2, max=pi/2)) d.psis <- psis[2]-psis[1] hs <- sin(psis) rs <- sqrt(1-hs^2) ls <- 2*pi*rs n.thetas <- floor(ls/d.psis) thetas[[1]] <- c(2*pi) for(i in 2:(n.psi-1)){ thetas[[i]] <- seq(from=0,to=2*pi,length=n.thetas[i]+1) thetas[[i]] <- thetas[[i]][-(n.thetas[i]+1)] } thetas[[n.psi]] <- c(2*pi) sapply(thetas,length) bridge <- list() for(i in 1:(n.psi-1)){ a <- c(thetas[[i]],2*pi) b <- c(thetas[[i+1]],2*pi) bridge[[i]] <- matrix(c(1,1),1,2) loop <- TRUE while(loop){ n.r <- nrow(bridge[[i]]) id.a <- bridge[[i]][n.r,1] + 1 id.b <- bridge[[i]][n.r,2] + 1 if(id.a > length(thetas[[i]]) & id.b > length(thetas[[i+1]])){ if(id.a-1!=1 & id.b-1!=1){ bridge[[i]] <- rbind(bridge[[i]],c(1,id.b-1)) } loop <- FALSE }else{ if(id.a > length(thetas[[i]])){ tmp <- c(id.a-1,id.b) }else if(id.b > length(thetas[[i+1]])){ tmp <- c(id.a,id.b-1) }else{ if(a[id.a] < b[id.b]){ tmp <- c(id.a,id.b-1) }else{ tmp <- c(id.a-1,id.b) } } bridge[[i]] <- rbind(bridge[[i]],tmp) } } } xyz <- matrix(0,0,3) edge <- matrix(0,0,2) triangles <- matrix(0,0,3) n.triangles <- rep(0,n.psi-1) for(i in 1:n.psi){ n.r <- nrow(xyz) if(i > 1){ pre <- (n.r-length(thetas[[i-1]])+1):n.r post <- (n.r+1):(n.r+length(thetas[[i]])) edge <- rbind(edge,cbind(post,c(post[-1],post[1]))) br <- bridge[[i-1]] new.edge <- cbind(pre[br[,1]],post[br[,2]]) edge <- rbind(edge,new.edge) tmp.tri <- cbind(new.edge,rbind(new.edge[-1,],new.edge[1,])) tmp <- apply(tmp.tri,1,unique) triangles <- rbind(triangles,t(tmp)) n.triangles[i-1] <- length(tmp[1,]) } psi <- psis[i] theta <- thetas[[i]] xyz <- rbind(xyz,cbind(cos(psi) * cos(theta),cos(psi)*sin(theta),sin(psi))) } return(list(xyz=xyz,edge=edge,triangles=triangles,n.psi=n.psi,thetas=thetas,n.triangles=n.triangles)) } knit_hooks$set(rgl = hook_rgl) n.psi <- 25 sp.mesh <- my.sphere.tri.mesh(n.psi) spxyz <- sp.mesh$xyz apply(spxyz^2,1,sum) xx <- 0.8 yy <- 0.1 large.xy1 <- which(spxyz[,1]>xx & spxyz[,2]<=yy) large.xy2 <- which(spxyz[,1]>xx & spxyz[,2]>yy) spxyz[large.xy1,1] <- xx + ((spxyz[large.xy1,1]-xx)*10)^2 spxyz[large.xy2,1] <- xx + ((spxyz[large.xy2,1]-xx)*10)^2 spxyz[large.xy2,2] <- yy + ((spxyz[large.xy2,1]-yy)*10)^5 plot3d(spxyz) segments3d(spxyz[c(t(sp.mesh$edge)),]) mesh.tri <- tmesh3d(t(spxyz),t(sp.mesh$triangles),homogeneous=FALSE) # 打点して plot3d(sp.mesh$xyz) # 三角形メッシュオブジェクトを作り mesh.tri <- tmesh3d(t(sp.mesh$xyz),t(sp.mesh$triangles),homogeneous=FALSE) # 三角形を灰色で塗る shade3d(mesh.tri,col="gray")
四元数でごそごそしていたので、ガウス曲率の計算に四元数を入力する仕様になっている。
# 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, na.rm=TRUE) #} return(Gcurvature) } # すべての頂点について計算 b <- mapply(function(j) gaussian_curvature(j, vertices, faces.v), seq(length(vertices))) cut0 <- cut(b, quantile(b, c(0, 0.02, 0.98, 1))) cols <- bluered(length(levels(cut0)))[cut0] plot3d(mesh.tri, type="dits", axes=FALSE, xlab="", ylab="", zlab="") shade3d(mesh.tri, col="grey") wire3d(mesh.tri) spheres3d(t(mesh.tri$vb), col=cols, radius=0.02)
