その際に、MapWindows GIS というオープンソースGISのAPIを拡張し、
VB.NET で、自分のための「なんちゃって自作GIS」を作ってみたわけです。
しかし…今でも後悔している。Active Xとか、Map Windows GIS のバージョンとか…。
あの時には、快適に動いていたのだが今は動かない。今は、「ごみ」コード。
もちろん、MapWindow GISは悪く無い。確かに、使いやすかった。
敗因は、見た目とユーザビリティに拘り過ぎたこと。商売する訳でもないのに…。
そして、私は結論に達した「いっそのこと、GUIやめちゃえば？」と。
最近は、さらに、考えが進んで「いっそ、可視化も捨てちゃえば？」になりつつある。
「可視化するから惑わされるのでは？」。そもそも、コンピュータは、目が見えないし。
場合によっては可視化を捨てるという選択はある。いや、独り言だから気にしない。
まぁ、良いや。形状分析の話だった。とにかく、私は頑張った。
VBを駆使して、楕円フーリエ記述子を用いて二次元形状を分析し、
それぞれの形状の類似度を比較する方法を実装した。
この方法、元々は「Morpho Metrics」という分野で発達した方法。
この方法、それほど特殊なものでは無いのだけれど、
地理空間データ（Shapefile）にも使えるようにしてみたのである。
っで、最初に書いたように、VBで書いたものは今はクズコード同然。
ちょっとした野暮用で、Rに再コーディングしたので、ここに掲載。
本当は、どこかで発表するつもりだったのだけれど、その機（気）を逃しので…。
詳しい（？）仕組みに関しては、以下のスライドにある通り。
このスライド、自分用のノートでなぐり書き状態。英語が可笑しいのは気にしない。
まぁ、何も無いよりかはマシという程度で…。
そうそう、ライセンスを気にする人もいるのだった。独占防止の意味で[CC(0)]。
私は、細かいことが嫌いだ。許可とかいらないので適当に。
後のソースコードに関しても。引用元を表示するかは「良心協定」ということで。
Simply shape from Yu Fujimoto
さて、前置きが長くなったが、以下にコードをペタペタする。
最初のコードは、楕円フーリエ変換を行い、楕円フーリエ記述子を取得する関数。
この結果を主成分分析やクラスタ分析にかけると形状の類似度が分析できる。
引数は、フーリエ級数から取り出すべき係数の数。既定は「20」。
増やした方が、形状の再現性は高くなるが、後ろの方は高周波なのでノイズかも。
無駄に多くしても計算量が増えるだけ。状況によるが、類似度の分析なら20〜50程度で十分。
以下が実際のコード。ここで全部書くと大変だからまずは、メインの部分だけ。
実際にフーリエ記述子を取り出す関数は、すぐ後の「baika.efd」の部分。
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ # Name: # baikaOne # Description: # Elliptic Fourier Descripter for shapefile. # Parameters: # harmo: Number of harmonix for elliptic fourier descriptors. Default is 20. # poly: Input polygon dataframes. # Values: # Returns elliptic fouier descriptors matrix. # Example: # shp <- readShapePoly("path/to/the/shapefilename.shp") # baika(shp,harmo=20) # _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ baikaOne <- function(polys, harmo=20, normalized=TRUE){ require(sp) require(maptools) # Counts number of polygons n <- length(polys) # Create Prograss Bar. pb <- txtProgressBar(min=0,max=n) # Row labels are acquired from shape IDs. nam <- character(n) # Column names. clb <- character(harmo*4) res <- list() res$efds <- list() res$matrix <- matrix(numeric(0),n,(harmo*4)) for (i in seq(1:n)){ # Update progress bar. setTxtProgressBar(pb,i) # Extracts control points from the polygon object. p <- length(polys@polygons[[i]]@Polygons) for (j in seq(1:p)){ if(polys@polygons[[i]]@Polygons[[j]]@hole==FALSE){ ctrlPts <- polys@polygons[[i]]@Polygons[[j]]@coords break() } } # Gets Identifier of the polygon IDvar <- polys@polygons[[i]]@ID nam[i] <- IDvar # Calculates EFD. if (normalized==TRUE){ res$efds[[i]] <- baika.efd(IDvar, ctrlPts,harmo=harmo,normalized=TRUE) } else { res$efds[[i]] <- baika.efd(IDvar, ctrlPts,harmo=harmo,normalized=FALSE) } # Take a name from the shape ID. for (j in seq(1:harmo)){ # Input each EFD to the matrix. k <- j+(j-1)*3 res$matrix[i,k] <- res$efds[[i]]$descriptors$A[j] res$matrix[i,(k+1)] <- res$efds[[i]]$descriptors$B[j] res$matrix[i,(k+2)] <- res$efds[[i]]$descriptors$C[j] res$matrix[i,(k+3)] <- res$efds[[i]]$descriptors$D[j] # Make a list of column labels if i is 1. if (i==1){ clb[k] <- paste("A",j,sep="") clb[k+1] <- paste("B",j,sep="") clb[k+2] <- paste("C",j,sep="") clb[k+3] <- paste("D",j,sep="") } } } # Asigns columun and row names to the matrix. rownames(res$matrix) <-nam colnames(res$matrix) <- clb # Returns matrix. return(invisible(res)) }
baika.efd <- function(IDvar, ctrlPts, harmo=20, normalized=TRUE){ require(sp) require(maptools) res <- list() if (harmo<1) { # If harmonix is less than 1 then exit. exit() } n <- nrow(ctrlPts) ti <- baika.lenFromStart(ctrlPts) l <- ti[n] res$id <- IDvar res$NumsOfPoints <- n res$totalLength <- l res$EdgeLengthFromStart <- ti res$centroid <- c(0,0) res$harmonix <- harmo res$parameters <- list() res$parameters$theta <- numeric(length=harmo) res$parameters$alpha <- numeric(length=harmo) res$parameters$beta <- numeric(length=harmo) res$parameters$phi <- numeric(length=harmo) res$parameters$ee <- numeric(length=harmo) res$descriptors <- list() res$descriptors$IsNormalized <- normalized res$descriptors$N <- harmo res$descriptors$A <- numeric(length=harmo) res$descriptors$B <- numeric(length=harmo) res$descriptors$C <- numeric(length=harmo) res$descriptors$D <- numeric(length=harmo) for (i in seq(1:harmo)){ Ai <- 0 Bi <- 0 Ci <- 0 Di <- 0 for (j in seq(1:(n-1))){ t1 <- ti[j] t2 <- ti[j+1] CosT1 <- cos((2*pi*i*t1)/l) CosT2 <- cos((2*pi*i*t2)/l) SinT1 <- sin((2*pi*i*t1)/l) SinT2 <- sin((2*pi*i*t2)/l) # Calculate Elliptic Fourier Descriptor for F(x) thisXi = ctrlPts[j,1] NextXi <- ctrlPts[j+1,1] Ai <- Ai+((CosT2-CosT1)*(((NextXi-thisXi)*l)/(t2-t1)))/(4*(pi^2)*(i^2)) Bi <- Bi+((SinT2-SinT1)*(((NextXi-thisXi)*l)/(t2-t1)))/(4*(pi^2)*(i^2)) # Calculate Elliptic Fourier Descriptor for G(x) thisYi <- ctrlPts[j,2] nextYi <- ctrlPts[j+1,2] Ci <- Ci+((CosT2-CosT1) * (((nextYi-thisYi)*l)/(t2-t1)))/(4*(pi^2)*(i^2)) Di <- Di+((SinT2-SinT1) * (((nextYi-thisYi)*l)/(t2-t1)))/(4*(pi^2)*(i^2)) } if (normalized==TRUE){ X <- Ai^2+Ci^2-Bi^2-Di^2 Y <- 2*(Ai*Bi+Ci*Di) pram_theta <- 0.5*(atan2(Y, X)) pram_alpha <- Ai*cos(pram_theta)+Bi*sin(pram_theta) pram_beta <- Ci*cos(pram_theta)+Di*sin(pram_theta) param_phi <- atan2(pram_beta,pram_alpha) param_ee <- sqrt(pram_alpha^2+pram_beta^2) res$parameters$theta[i] <- pram_theta res$parameters$alpha[i] <- pram_alpha res$parameters$beta[i] <- pram_beta res$parameters$phi[i] <- param_phi res$parameters$ee[i] <- param_ee SinTheta <- sin(pram_theta*i) CosTheta <- cos(pram_theta*i) SinPhi <- sin(param_phi) CosPhi <- cos(param_phi) res$descriptors$A[i]<-((CosPhi*Ai+SinPhi*Ci)*(CosTheta)+(CosPhi*Bi+SinPhi*Di)*SinTheta)/param_ee res$descriptors$B[i]<-((CosPhi*Ai+SinPhi*Ci)*(-1*SinTheta)+(CosPhi*Bi+SinPhi*Di)*CosTheta)/param_ee res$descriptors$C[i]<-((-1*SinPhi*Ai+CosPhi*Ci)*(CosTheta)+(-1*SinPhi*Bi+CosPhi*Di)*SinTheta)/param_ee res$descriptors$D[i]<-((-1*SinPhi*Ai+CosPhi*Ci)*(-1*SinTheta)+(-1*SinPhi*Bi+CosPhi*Di)*CosTheta)/param_ee } else { res$centroid[1] <- baika.centroid(ctrlPts, IDvar=IDvar)$coord[1] res$centroid[2] <- baika.centroid(ctrlPts, IDvar=IDvar)$coord[2] res$descriptors$A[i]<-Ai res$descriptors$B[i]<-Bi res$descriptors$C[i]<-Ci res$descriptors$D[i]<-Di } } return(invisible(res)) }
っで、次は、楕円フーリエ記述子から元の形状を復元する関数。
関数から復元するから、形状の構成点の配置は等間隔。引数には、構成点の数を指定。
実際に形状を復元しているのは、すぐ後ろの「baika.reconstruct」の部分。
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ # Name: # baikaOneinv # Description: # Build approximated polygon from elliptic fourier descriptors. # Parameters: # efds: efds are acquired by baika. # approx: Number of points for approximates. # Values: # Returnes SpatialPolygonDataFrame class # _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ baikaOneinv <- function(efds, approx){ require(sp) require(maptools) # Gets number of polygons n <- length(efds$efds) polys <- vector("list", n) # Creates the Prograss Bar. pb <- txtProgressBar(min=0,max=n) for (i in seq(1:n)){ # Update progress bar. setTxtProgressBar(pb,i) # Creates Polygon instance form coordinates of the approximated shape. poly <- Polygon(baika.reconstruct(efds$efds[i],approx),hole=FALSE) # Gets ID label from original IDs. fid <- as.character(efds$efds[i][[1]]$id) # Build SpatialPolygon instance from the shape. polys[[i]] <- Polygons(list(Polygon(poly)), ID=fid) } # Convert from SpatialPolygon dataframe to SpatialPolygonDataFrame. spPoly <- SpatialPolygons(polys) ids <- sapply(slot(spPoly, "polygons"), function(x) slot(x, "ID")) length(ids) df <- data.frame(rep(0,length(ids)), row.names=ids) res <- SpatialPolygonsDataFrame(spPoly,df) # Return SpatialPolygonDataFrame. return(invisible(res)) }
baika.reconstruct <- function(efd, approx){ h <- efd[[1]]$harmonix res <- matrix(nrow=(approx+1),ncol=2) A1 <- efd[[1]]$descriptors$A[1] B1 <- efd[[1]]$descriptors$B[1] C1 <- efd[[1]]$descriptors$C[1] D1 <- efd[[1]]$descriptors$D[1] for (i in seq(1:approx)){ theta <- (i*2*pi)/approx sX1 <- A1*cos(theta)+B1*sin(theta) sY1 <- C1*cos(theta)+D1*sin(theta) sXi <- 0 sYi <- 0 for (j in seq(2:h)){ Aj <- efd[[1]]$descriptors$A[j] Bj <- efd[[1]]$descriptors$B[j] Cj <- efd[[1]]$descriptors$C[j] Dj <- efd[[1]]$descriptors$D[j] sXi = sXi+Aj*cos(j*theta)+Bj*sin(j*theta) sYi = sYi+Cj*cos(j*theta)+Dj*sin(j*theta) } res[i,1] <- sX1+sXi+efd[[1]]$centroid[1] res[i,2] <- sY1+sYi+efd[[1]]$centroid[2] } res[(approx+1),1] <- res[1,1] res[(approx+1),2] <- res[1,2] return(invisible(res)) }
上のコード、何に使うの？って言う人もいるかも。実は、楕円フーリエによる方法は、
アメーバみたいな不整形な形にはあまり向かない。第一近似楕円の軸回転の問題も。
っで、そういったケースで上の方法を応用するとプロクラテス変換が使える。
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ # Name: # baikaTwo # Description: # Procrustes Analysis for spatial polygon dataframe created by baikaOneInv. # Parameters: # poly: Input polygon dataframes. # Values: # Returns elliptic fouier descriptors matrix. # Example: # shp <- readShapePoly("path/to/the/shapefilename.shp") # baika(shp,harmo=20) # _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ baikaTwo <- function(polys){ require(vegan) require(sp) require(maptools) # Counts number of polygons n <- length(polys) # Create Prograss Bar. pb <- txtProgressBar(min=0,max=n) # Row labels are acquired from shape IDs. rNams <- character(n) cNams <- character(n) res <- matrix(nrow=n,ncol=n) for (i in seq(1:n)){ # Update progress bar. setTxtProgressBar(pb,i) # Extracts control points from the polygon object. p <- length(polys@polygons[[i]]@Polygons) for (k in seq(1:p)){ if(polys@polygons[[i]]@Polygons[[k]]@hole==FALSE){ A <- polys@polygons[[i]]@Polygons[[k]]@coords break() } } # Gets Identifier of the polygon rNams[i] <- polys@polygons[[i]]@ID for (j in seq(1:n)){ # Extracts control points from the polygon object. p <- length(polys@polygons[[j]]@Polygons) for (k in seq(1:p)){ if(polys@polygons[[j]]@Polygons[[k]]@hole==FALSE){ B <- polys@polygons[[j]]@Polygons[[k]]@coords break() } } # Gets Identifier of the polygon if (i==1){ cNams[j] <- polys@polygons[[j]]@ID } res[i,j] <- procrustes(A,B,symmetric=TRUE)$ss } } # Asigns columun and row names to the matrix. rownames(res) <- rNams colnames(res) <- cNams # Returns matrix. return(invisible(res)) }
っで、以下には、中途に登場した自作関数を順番に。
実は、細々とした所で使っているはず。多分。
一年以上も前の話なので、少々、忘れていることが多い。
まずは、フーリエ記述子をつかって面積を求める部分。
baika.area <- function(efds){ res <- 0 n <- length(efds$A[i]) for (i in seq(1:n)){ Ai <- efds$A[i] Bi <- efds$B[i] Ci <- efds$C[i] Di <- efds$D[i] res <- res+((Ai*Di)-(Ci*Bi))*i*pi } return(invisible(res)) }
っで、形状の中心を求める関数。
baika.centroid <- function(ctrlPts, IDvar=1){ n <- nrow(ctrlPts) ti <- baika.lenFromStart(ctrlPts) l <- ti[n] res <- list() res$id <- IDvar res$coord <- c(0,0) for (i in seq(1:(n-1))){ thisXi <- ctrlPts[i,1] thisYi <- ctrlPts[i,2] nextXi <- ctrlPts[i+1,1] nextYi <- ctrlPts[i+1,2] dist <- ti[i+1]-ti[i] res$coord[1] <- res$coord[1]+((nextXi+thisXi)/2*(dist)) res$coord[2] <- res$coord[2]+((nextYi+thisYi)/2*(dist)) } res$coord[1] <- res$coord[1]/l res$coord[2] <- res$coord[2]/l return(invisible(res)) }
これは、何かの中間処理で使っている関数。起点となる点を取得する。
baika.startPoint <- function(ctrlPts, IDvar=1){ res <- list() res$id <- IDvar res$coord <- c(0,0) res$coord[1] <- ctrlPts[1,1] res$coord[2] <- ctrlPts[1,2] return(invisible(res)) }
次の関数は、上記の起点からの累積距離を得るための関数。
baika.lenFromStart <- function(ctrlPts){ # Count the number of control points. n <- nrow(ctrlPts) # Get the first set of coordinates of the first point. sx <- ctrlPts[1,1] sy <- ctrlPts[1,2] res <- {} res[1] <- 0 # Get length from start point at each control point. for (i in seq(2:n)){ thisXi <- ctrlPts[i,1] thisYi <- ctrlPts[i,2] nextXi <- ctrlPts[i+1,1] nextYi <- ctrlPts[i+1,2] distX <- thisXi-nextXi distY <- thisYi-nextYi res[i+1] <- res[i]+(sqrt(distX^2+distY^2)) } return(invisible(res)) }
まぁ、コードの説明に関しては非常に中途半端ではあるけれど、
解らなければ個別に聞いてください。その都度、説明を更新する予定。
ところで、ひとつずつコピペするのは大変だろうし、ソースコードを共有。
https://docs.google.com/open?id=0B9OtIs5BjRVhMFIxR2VmXzRPUDQ
0 件のコメント:
コメントを投稿