vtk를 이용한 contour plot 구현

set R 0
set Z 0

set in [open vtka.in r]

gets $in line
scan $line "%d %d" R Z
close $in

vtkFloatArray rcoords
vtkFloatArray zcoords
vtkFloatArray scalars

vtkRectilinearGrid rgrid
        rgrid SetDimensions $R $Z 1
vtkScalarBarActor scalarBar
vtkCubeAxesActor2D axes

set f [open ne.dat r]

for {set i 1} {$i <= [expr ($R+$Z)*($Z-1)]} {incr i 1} {
                gets $f line
                foreach j {1 2 3} {
                        set arr($i,$j) [lindex $line [expr $j - 1]]
                }      
        }
       
set n 1
set m 1
for {set a 1} {$a <= $Z} {incr a 1} {
        for {set b 1} {$b <= $R} {incr b 1} {          
               
                if {$b == 1} {
                        foreach y $arr($m,1) {
                                puts $y
                        }
                        zcoords InsertNextValue $y
                }
                if {$a == 1} {
                        foreach x $arr($b,2) {
                                puts $x
                        }
                        rcoords InsertNextValue $x
                }
                foreach s $arr($n,3) {
                        puts $s
                }
                scalars InsertNextValue $s
                incr n 1
               
        }
        incr m [expr $R+1]
        incr n 1
       
}                              
                               
close $f

set min [scalars GetValue 0]
        set max [scalars GetValue 0]

        for {set i 0} {$i < [expr $R*$Z]} { incr i 1 } {
                if { [scalars GetValue $i] < $min } {          
                        set min [scalars GetValue $i]
                }
        }
        for {set i 0} {$i < [expr $R*$Z]} { incr i 1 } {
                if { [scalars GetValue $i] > $max } {          
                        set max [scalars GetValue $i]
                }
        }

        set rmin [rcoords GetValue 0]
        set rmax [rcoords GetMaxNorm]

        for {set i 0} {$i < $R} {incr i 1} {
                if { [rcoords GetValue $i] < $rmin } {
                        set rmin [rcoords GetValue $i]
                }
        }

        set zmin [zcoords GetValue 0]
        for {set i 0} {$i < $Z} {incr i 1} {
                if { [zcoords GetValue $i] < $zmin } {
                        set zmin [zcoords GetValue $i]
                }
        }
        set zmax [zcoords GetValue 0]
        for {set i 0} {$i < $Z} {incr i 1} {
                if { [zcoords GetValue $i] > $zmax } {
                        set zmax [zcoords GetValue $i]
                }
        }

rgrid SetXCoordinates rcoords
rgrid SetYCoordinates zcoords
[rgrid GetPointData] SetScalars scalars

vtkLookupTable lut
        lut SetNumberOfColors 256
        lut SetHueRange 0.667 0.0
        lut Build

vtkDataSetMapper mapper
        mapper SetInput rgrid
        mapper SetLookupTable lut
        mapper SetScalarRange $min $max
vtkActor actor
        actor SetMapper mapper

#

vtkRenderer renderer
        renderer SetBackground 0 0 0
        renderer AddActor actor
        renderer AddActor2D scalarBar
        renderer AddActor2D axes
        [renderer GetActiveCamera] SetPosition 13.0 5.0 50.0
        [renderer GetActiveCamera] SetFocalPoint 13.0 5.0 0.0

#create ScalarBar
scalarBar SetLookupTable lut
scalarBar SetTitle "Density"
[scalarBar GetPositionCoordinate] SetCoordinateSystemToNormalizedViewport
[scalarBar GetPositionCoordinate] SetValue 0.1 0.01
scalarBar SetOrientationToHorizontal
scalarBar SetWidth 0.8
scalarBar SetHeight 0.1

#create CubeAxes
axes SetInput rgrid
axes SetFlyModeToOuterEdges
axes SetFontFactor 0.7
axes ScalingOff
axes SetXLabel "R cm"
axes SetYLabel "Z cm"
axes SetCornerOffset 0.0
axes SetNumberOfLabels 5
axes SetCamera [renderer GetActiveCamera]
axes SetLabelFormat "%6.4g"
axes ZAxisVisibilityOff

vtkRenderWindow renwin
        renwin AddRenderer renderer
#       renwin SetSize 500 500

vtkRectilinearGridWriter writer
        writer SetInput rgrid
        writer SetFileName "rectilinear.vtk"
        writer Write

set vtkw [vtkTkRenderWidget .ren -width 500 -height 400 -rw renwin]

#Setup Tk bindings and vtk obserbers for that widget.
::vtk::bind_tk_render_widget $vtkw

frame .params

ttk::button .params.quit -text "Quit" -command ::vtk::cb_exit
pack .params.quit -side top -anchor nw -fill both
pack $vtkw .params -side top -fill both -expand yes

wm protocol . WM_DELETE_WINDOW ::vtk::cb_exit

tkwait window .