四元数とベクトルとメッシュと面積と

3次元物体をあつかっている。
3次元座標を四元数で扱うと回転や行列計算などが簡単らしい。
いま、頂点が三角メッシュ化されて表面が形成されている.obj のうさぎがいるとする。

いま、表面の各三角形の面積の合計が、四元数と疎行列を用いてこんな感じで計算される。

# 特に疎行列において、面単位で不定回数の値加算をするためのユーティリティ関数
my.vector.access <- function(v,a,func=sum,zero=0){
  if(is.vector(v)){
        v <- matrix(v,ncol=1)
    }
    ord <- order(a)
    rle.out <- rle(a[ord])
    num.row <- length(rle.out[[1]])
    num.col <- max(rle.out[[1]])
    tmp1 <- rep(1:num.row,rle.out[[1]])
    tmp2 <- c()
    for(i in 1:num.row){
        tmp2 <- c(tmp2,1:rle.out[[1]][i])
    }
    addr <- tmp1 + num.row*(tmp2-1)
    ret.v <- matrix(0,num.row,ncol(v))
    for(i in 1:ncol(v)){
        if(zero==0){
            tmp.V <- sparseVector(v[ord,i],i=addr,length=num.row*num.col)
            M <- Matrix(tmp.V,num.row,num.col)
        }else{
            M <- matrix(zero,num.row,num.col)
            M[addr] <- v[ord,i]

        }
        ret.v[,i] <- apply(M,1,func)

    }
    return(list(A = rle.out[[2]],V = ret.v))
}

また、三角形の面積は、四元数の積の虚部について絶対値をとって2で割ればよい。

Mod(Im(quanterion1- quanterion0))

　
xyz 座標での三角形が空間上になす平面の面積は、

で与えられるから、8つの三角形を構成している左足の裏を見てみると

i <- 1
bun.v.new <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
plot3d(bun.v.new)
mesh.tri <- tmesh3d(t(bun.v.new),t(bun.f),homogeneous=FALSE)
# 三角形を灰色で塗る
shade3d(mesh.tri,col="gray")
wire3d(mesh.tri)
j <- 337
text3d(bun.v.new[j, ], texts=j, col="blue")

# j 番目の点を構成成分とする三角形を抽出
k <- setdiff(unique(unlist(bun.f[apply(bun.f==j, 1, any),])), j)
# 点に目印をつける
for(l in k) text3d(bun.v.new[l, ], texts=l, col="red")
# 三角形を抽出する
plane <- as.matrix(bun.f[apply(bun.f==j, 1, any),])
# 三角形の面積の計算
tri_area <- rep(0, nrow(plane))
for(p in seq(nrow(plane))){
  hogemat <- as.matrix(bun.v[plane[p,],])
  tri_area[p] <- sqrt(sum(mapply(function(r) det(cbind(hogemat[, c(r, r%%3+1)], 1))^2 ,seq(3))))/2
}

M[j,]

[1] 0.0002846982 0.0002755177 0.0003673656 0.0003241680 0.0004353401
[6] 0.0004162023 0.0003144500 0.0004536900

tri_area

[1] 0.0002846982 0.0002755177 0.0003673656 0.0003241680 0.0004353401
[6] 0.0004162023 0.0003144500 0.0004536900

　
合ってる。
だがしかし、いま行っているC++ によるものとは計算時間ががが