VTKで球を描く (手作り)

vtkSphereSourceは便利だが,勉強にならないので,球をセルから構成してみる.

定数

package require vtk

set nlon 100
set nlat [expr {$nlon/2 + 1}]
set pi [expr {acos(-1.)}]
set pi2 [expr {$pi * 2.}]
set npts [expr {$nlon * $nlat}]
#puts stdout "pi=$pi pi2=$pi2 npts=$npts"

経度φと余緯度θからデカルト座標に変換するコマンドを作る.
VTKでは,経度がθ余緯度がφのようだが気にしない.

proc phitheta2xyz args {
set phi [lindex $args 0]
set theta [lindex $args 1]
set xyz {}
set sin_theta [expr {sin($theta)}]
lappend xyz [expr {cos($phi) * $sin_theta}]
lappend xyz [expr {sin($phi) * $sin_theta}]
lappend xyz [expr {cos($theta)}]
}

経度,余緯度を生成.

set lon {}
set dlon [expr {$pi/$nlon}]
for {set i 0} {$i<$nlon} {incr i} {
lappend lon [expr {$pi2 * $i / $nlon}]
}
#puts stdout $lon
set lat {}
set dlat [expr {$pi/($nlat-1)}]
for {set j 0} {$j<$nlat} {incr j} {
lappend lat [expr {$dlat * $j}]
}
#puts stdout $lat

座標点のデカルト座標を生成してvtkPointsに登録.

vtkPoints pts
pts SetNumberOfPoints $npts
set i 0
foreach theta $lat {
foreach phi $lon {
# puts stdout "$phi $theta"
# puts [phitheta2xyz $phi $theta]
eval pts SetPoint [concat $i [phitheta2xyz $phi $theta]]
incr i
}
}

4点からセルを作る.経度方向に周期境界.$nlon-1番目の隣は,0番目に戻る.

vtkCellArray clls
set nlonm1 [expr {$nlon - 1}]
for {set j 0} {$j<$nlat-1} {incr j} {
for {set i 0} {$i<$nlon} {incr i} {
set k0 [expr {$nlon*$j + $i}]
if {$i!=$nlonm1} {
set k1 [expr {$k0 + 1}]
} else {
set k1 [expr {$k0 - $nlonm1}]
}
set k2 [expr {$k1 + $nlon}]
set k3 [expr {$k0 + $nlon}]
clls InsertNextCell 4
clls InsertCellPoint $k0
clls InsertCellPoint $k1
clls InsertCellPoint $k2
clls InsertCellPoint $k3
}
}

点とセルをvtkPolyDataに登録.

vtkPolyData sfc
sfc SetPoints pts
sfc SetPolys clls
pts Delete
clls Delete
ファイルに保存して完了.

vtkPolyDataWriter writer
writer SetInput sfc
writer SetFileName sphere2.vtk
writer Write
writer Delete
sfc Delete

コメント

このブログの人気の投稿

静的サイト生成ツールの導入

uptex

TeXのソースを分割してTeXShopでタイプセット