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