source: trunk/gui/scripts/field.tcl @ 3960

Last change on this file since 3960 was 3960, checked in by gah, 11 years ago

initial take on dicom reader in field.tcl
File size: 58.4 KB
Line 
1# -*- mode: tcl; indent-tabs-mode: nil -*-
2
3# ----------------------------------------------------------------------
4#  COMPONENT: field - extracts data from an XML description of a field
5#
6#  This object represents one field in an XML description of a device.
7#  It simplifies the process of extracting data vectors that represent
8#  the field.
9# ======================================================================
10#  AUTHOR:  Michael McLennan, Purdue University
11#  Copyright (c) 2004-2012  HUBzero Foundation, LLC
12#
13#  See the file "license.terms" for information on usage and
14#  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
15# ======================================================================
16
17# TODO:
18#    Vector field limits are wrong: need to compute magnitude limits and
19#    component-wise limits.
20
21#
22# Possible field dataset types:
23#
24# 2D Datasets
25#       vtk             (range of z-axis is zero).
26#       unirect2d       (deprecated)
27#       cloud           (x,y point coordinates) (deprecated)
28#       mesh
29# 3D Datasets
30#       vtk
31#       unirect3d       (deprecated)
32#       cloud           (x,y,z coordinates) (deprecated)
33#       mesh
34#       dx              (FIXME: make dx-to-vtk converter work)
35#       ucd avs
36#
37# Viewers:
38#       Format     Dim  Description                     Viewer          Server
39#       vtk         2   vtk file data.                  contour         vtkvis
40#       vtk         3   vtk file data.                  isosurface      vtkvis
41#       mesh        2   points-on-mesh                  heightmap       vtkvis
42#       mesh        3   points-on-mesh                  isosurface      vtkvis
43#       dx          3   DX                              volume          nanovis
44#       unirect2d   2   unirect3d + extents > 1 flow    flow            nanovis
45#       unirect3d   3   unirect2d + extents > 1 flow    flow            nanovis
46#       
47# With <views>, can specify which viewer for specific datasets.  So it's OK
48# for the same dataset to be viewed in more than one way.
49#  o Any 2D dataset can be viewed as a contour/heightmap.
50#  o Any 3D dataset can be viewed as a isosurface.
51#  o Any 2D dataset with vector data can be streamlines or flow. 
52#  o Any 3D uniform rectilinear dataset can be viewed as a volume.
53#  o Any 3D dataset with vector data can be streamlines or flow.
54#
55# Need <views> to properly do things like qdot: volume with a polydata
56# transparent shell.  The view will combine the two objects <field> and
57# <drawing> (??) into a single viewer.
58#
59package require Itcl
60package require BLT
61
62namespace eval Rappture {
63    # forward declaration
64}
65
66itcl::class Rappture::Field {
67    protected variable _dim 0;          # Dimension of the mesh
68    private variable _xmlobj "";        # ref to XML obj with field data
69    private variable _limits;           # maps axis name => {z0 z1} limits
70    private variable _field ""
71    private variable _comp2fldName ;    # cname => field names.
72    private variable _comp2type ;       # cname => type (e.g. "vectors")
73    private variable _comp2size ;       # cname => # of components in element
74    private variable _comp2assoc;       # cname => association (e.g. pointdata)
75    private variable _fld2Components;   # field name => number of components
76    private variable _fld2Label;        # field name => label
77    private variable _fld2Units;        # field name => units
78    private variable _hints
79    private variable _viewer "";        # Hints which viewer to use
80    private variable _xv "";            # For 1D meshes only.  Holds the points
81    private variable _isValid 0;        # Indicates if the field contains
82                                        # valid data.
83    private variable _isValidComponent; #  Array of valid components found
84    constructor {xmlobj path} {
85        # defined below
86    }
87    destructor {
88        # defined below
89    }
90    public method blob { cname }
91    public method components {args}
92    public method controls {option args}
93    public method extents {{cname -overall}}
94    public method numComponents {cname}
95    public method fieldlimits {}
96    public method valueLimits { cname }
97    public method flowhints { cname }
98    public method hints {{key ""}}
99    public method isunirect2d {}
100    public method isunirect3d {}
101    public method limits {axis}
102    public method mesh {{cname -overall}}
103    public method style { cname }
104    public method type {}
105    public method values { cname }
106    public method vtkdata {cname}
107
108    public method fieldnames { cname } {
109        if { ![info exists _comp2fldName($cname)] } {
110            return ""
111        }
112        return $_comp2fldName($cname)
113    }
114    public method fieldinfo { fname } {
115        lappend out $_fld2Label($fname)
116        lappend out $_fld2Units($fname)
117        lappend out $_fld2Components($fname)
118        return $out
119    }
120    public method isvalid {} {
121        return $_isValid
122    }
123    public method viewer {} {
124        return $_viewer
125    }
126    protected method Build {}
127    protected method _getValue {expr}
128
129    private variable _path "";          # Path of this object in the XML
130    private variable _units ""   ;      # system of units for this field
131    private variable _zmax 0     ;# length of the device
132
133    private variable _comp2dims  ;# maps component name => dimensionality
134    private variable _comp2xy    ;# maps component name => x,y vectors
135    private variable _comp2vtk   ;# maps component name => vtk file data
136    private variable _comp2dx    ;# maps component name => OpenDX data
137    private variable _comp2unirect2d ;# maps component name => unirect2d obj
138    private variable _comp2unirect3d ;# maps component name => unirect3d obj
139    private variable _comp2style ;# maps component name => style settings
140    private variable _comp2cntls ;# maps component name => x,y control points
141    private variable _comp2extents
142    private variable _comp2limits;      #  Array of limits per component
143    private variable _type ""
144    private variable _comp2flowhints
145    private variable _comp2mesh
146    private common _counter 0    ;# counter for unique vector names
147
148    private method AvsToVtk { cname contents }
149    private method DicomSeriesToVtk { cname contents }
150    private method BuildPointsOnMesh { cname }
151    private method ConvertToVtkData { cname }
152    protected method GetAssociation { cname }
153    protected method GetTypeAndSize { cname }
154    protected method ReadVtkDataSet { cname contents }
155    private method InitHints {}
156
157    private method VerifyVtkDataSet { contents }
158    private variable _values ""
159}
160
161# ----------------------------------------------------------------------
162# CONSTRUCTOR
163# ----------------------------------------------------------------------
164itcl::body Rappture::Field::constructor {xmlobj path} {
165    package require vtk
166    if {![Rappture::library isvalid $xmlobj]} {
167        error "bad value \"$xmlobj\": should be Rappture::library"
168    }
169    set _xmlobj $xmlobj
170    set _path $path
171    set _field [$xmlobj element -as object $path]
172    set _units [$_field get units]
173
174    set xunits [$xmlobj get units]
175    if {"" == $xunits || "arbitrary" == $xunits} {
176        set xunits "um"
177    }
178
179    # determine the overall size of the device
180    set z0 [set z1 0]
181    foreach elem [$_xmlobj children components] {
182        switch -glob -- $elem {
183            box* {
184                if {![regexp {[0-9]$} $elem]} {
185                    set elem "${elem}0"
186                }
187                set z0 [$_xmlobj get components.$elem.corner0]
188                set z0 [Rappture::Units::convert $z0 \
189                    -context $xunits -to $xunits -units off]
190
191                set z1 [$_xmlobj get components.$elem.corner1]
192                set z1 [Rappture::Units::convert $z1 \
193                    -context $xunits -to $xunits -units off]
194
195                set _limits($elem) [list $z0 $z1]
196            }
197        }
198    }
199    set _zmax $z1
200
201    # build up vectors for various components of the field
202    Build
203    InitHints
204}
205
206# ----------------------------------------------------------------------
207# DESTRUCTOR
208# ----------------------------------------------------------------------
209itcl::body Rappture::Field::destructor {} {
210    itcl::delete object $_field
211    # don't destroy the _xmlobj! we don't own it!
212
213    foreach name [array names _comp2xy] {
214        # Destroys both x and y vectors.
215        eval blt::vector destroy $_comp2xy($name)
216    }
217    foreach name [array names _comp2dx] {
218        set vector [lindex $_comp2dx($name) 1]
219        eval blt::vector destroy $vector
220    }
221    foreach name [array names _comp2unirect2d] {
222        itcl::delete object $_comp2unirect2d($name)
223    }
224    foreach name [array names _comp2unirect3d] {
225        itcl::delete object $_comp2unirect3d($name)
226    }
227    foreach name [array names _comp2flowhints] {
228        itcl::delete object $_comp2flowhints($name)
229    }
230    foreach name [array names _comp2mesh] {
231        # Data is in the form of a mesh and a vector.
232        foreach { mesh vector } $_comp2mesh($name) break
233        # Release the mesh (may be shared)
234        set class [$mesh info class]
235        ${class}::release $mesh
236        # Destroy the vector
237        blt::vector destroy $vector
238    }
239}
240
241# ----------------------------------------------------------------------
242# USAGE: components ?-name|-dimensions|-style? ?<pattern>?
243#
244# Returns a list of names or types for the various components of
245# this field.  If the optional glob-style <pattern> is specified,
246# then it returns only the components with names matching the pattern.
247# ----------------------------------------------------------------------
248itcl::body Rappture::Field::components {args} {
249    Rappture::getopts args params {
250        flag what -name default
251        flag what -dimensions
252        flag what -style
253        flag what -particles
254        flag what -flow
255        flag what -box
256    }
257
258    set pattern *
259    if {[llength $args] > 0} {
260        set pattern [lindex $args 0]
261        set args [lrange $args 1 end]
262    }
263    if {[llength $args] > 0} {
264        error "wrong # args: should be \"components ?switches? ?pattern?\""
265    }
266
267    # There's only one dimension of the field.  Components can't have
268    # different dimensions in the same field.  They would by definition be
269    # using different meshes and viewers.
270    if { $params(what) == "-dimensions" } {
271        return "${_dim}D"
272    }
273    # BE CAREFUL: return component names in proper order
274    set rlist ""
275    set components {}
276    # First compile a list of valid components that match the pattern
277    foreach cname [$_field children -type component] {
278        if { ![info exists _isValidComponent($cname)] } {
279            continue
280        }
281        if { [string match $pattern $cname] } {
282            lappend components $cname
283        }
284    }
285    # Now handle the tests.
286    switch -- $params(what) {
287        -name {
288            set rlist $components
289        }
290        -style {
291            foreach cname $components {
292                if { [info exists _comp2style($cname)] } {
293                    lappend rlist $_comp2style($cname)
294                }
295            }
296        }
297    }
298    return $rlist
299}
300
301# ----------------------------------------------------------------------
302# USAGE: mesh ?<name>?
303#
304# Returns a list {xvec yvec} for the specified field component <name>.
305# If the name is not specified, then it returns the vectors for the
306# overall field (sum of all components).
307# ----------------------------------------------------------------------
308itcl::body Rappture::Field::mesh {{cname -overall}} {
309    if {$cname == "-overall" || $cname == "component0"} {
310        set cname [lindex [components -name] 0]
311    }
312    if {[info exists _comp2xy($cname)]} {
313        return [lindex $_comp2xy($cname) 0]  ;# return xv
314    }
315    if { [info exists _comp2vtk($cname)] } {
316        # FIXME: extract mesh from VTK file data.
317        if { $_comp2dims($cname) == "1D" } {
318            return $_xv
319        }
320        error "method \"mesh\" is not implemented for VTK file data"
321    }
322    if {[info exists _comp2dx($cname)]} {
323        return ""  ;# no mesh -- it's embedded in the value data
324    }
325    if {[info exists _comp2mesh($cname)]} {
326        return ""  ;# no mesh -- it's embedded in the value data
327    }
328    if {[info exists _comp2unirect2d($cname)]} {
329        set mobj [lindex $_comp2unirect2d($cname) 0]
330        return [$mobj mesh]
331    }
332    if {[info exists _comp2unirect3d($cname)]} {
333        set mobj [lindex $_comp2unirect3d($cname) 0]
334        return [$mobj mesh]
335    }
336    error "can't get field mesh: Unknown component \"$cname\": should be one of [join [lsort [array names _comp2dims]] {, }]"
337}
338
339# ----------------------------------------------------------------------
340# USAGE: values ?<name>?
341#
342# Returns a list {xvec yvec} for the specified field component <name>.
343# If the name is not specified, then it returns the vectors for the
344# overall field (sum of all components).
345# ----------------------------------------------------------------------
346itcl::body Rappture::Field::values {cname} {
347    if {$cname == "component0"} {
348        set cname "component"
349    }
350    if {[info exists _comp2xy($cname)]} {
351        return [lindex $_comp2xy($cname) 1]  ;# return yv
352    }
353    # VTK file data
354    if { [info exists _comp2vtk($cname)] } {
355        # FIXME: extract the values from the VTK file data
356        if { $_comp2dims($cname) == "1D" } {
357            return $_values
358        }
359        error "method \"values\" is not implemented for vtk file data"
360    }
361    # Points-on-mesh
362    if { [info exists _comp2mesh($cname)] } {
363        set vector [lindex $_comp2mesh($cname) 1]
364        return [$vector range 0 end]
365    }
366    if {[info exists _comp2dx($cname)]} {
367        set vector [lindex $_comp2dx($cname) 1]
368        return [$vector range 0 end]
369    }
370    if {[info exists _comp2unirect2d($cname)]} {
371        return [$_comp2unirect2d($cname) values]
372    }
373    if {[info exists _comp2unirect3d($cname)]} {
374        return [$_comp2unirect3d($cname) blob]
375    }
376    error "can't get field values. Unknown component \"$cname\": should be [join [lsort [array names _comp2dims]] {, }]"
377}
378
379# ----------------------------------------------------------------------
380# USAGE: blob ?<name>?
381#
382# Returns a string representing the blob of data for the mesh and values.
383# ----------------------------------------------------------------------
384itcl::body Rappture::Field::blob {cname} {
385    if {$cname == "component0"} {
386        set cname "component"
387    }
388    if {[info exists _comp2xy($cname)]} {
389        return ""
390    }
391    if { [info exists _comp2vtk($cname)] } {
392        error "blob not implemented for VTK file data"
393    }
394    if {[info exists _comp2dx($cname)]} {
395        return [lindex $_comp2dx($cname) 0]
396    }
397    if {[info exists _comp2unirect2d($cname)]} {
398        set blob [$_comp2unirect2d($cname) blob]
399        lappend blob "values" $_values
400        return $blob
401    }
402    if {[info exists _comp2unirect3d($cname)]} {
403        return [$_comp2unirect3d($cname) blob]
404    }
405    error "can't get field blob: Unknown component \"$cname\": should be one of [join [lsort [array names _comp2dims]] {, }]"
406}
407
408itcl::body Rappture::Field::valueLimits { cname } {
409    if { [info exists _comp2limits($cname)] } {
410        return $_comp2limits($cname)
411    }
412    return ""
413}
414
415# ----------------------------------------------------------------------
416# USAGE: limits <axis>
417#
418# Returns a list {min max} representing the limits for the specified
419# axis.
420# ----------------------------------------------------------------------
421itcl::body Rappture::Field::limits {which} {
422    set min ""
423    set max ""
424    blt::vector tmp zero
425
426    foreach cname [array names _comp2dims] {
427        switch -- $_comp2dims($cname) {
428            1D {
429                switch -- $which {
430                    x - xlin {
431                        set pos 0; set log 0; set axis x
432                    }
433                    xlog {
434                        set pos 0; set log 1; set axis x
435                    }
436                    y - ylin - v - vlin {
437                        set pos 1; set log 0; set axis y
438                    }
439                    ylog - vlog {
440                        set pos 1; set log 1; set axis y
441                    }
442                    default {
443                        error "bad axis \"$which\": should be x, xlin, xlog, y, ylin, ylog, v, vlin, vlog"
444                    }
445                }
446
447                set vname [lindex $_comp2xy($cname) $pos]
448                $vname variable vec
449
450                if {$log} {
451                    # on a log scale, use abs value and ignore 0's
452                    $vname dup tmp
453                    $vname dup zero
454                    zero expr {tmp == 0}            ;# find the 0's
455                    tmp expr {abs(tmp)}             ;# get the abs value
456                    tmp expr {tmp + zero*max(tmp)}  ;# replace 0's with abs max
457                    set axisMin [blt::vector expr min(tmp)]
458                    set axisMax [blt::vector expr max(tmp)]
459                } else {
460                    set axisMin $vec(min)
461                    set axisMax $vec(max)
462                }
463
464                if {"" == $min} {
465                    set min $axisMin
466                } elseif {$axisMin < $min} {
467                    set min $axisMin
468                }
469                if {"" == $max} {
470                    set max $axisMax
471                } elseif {$axisMax > $max} {
472                    set max $axisMax
473                }
474            }
475            2D - 3D {
476                if {[info exists _comp2limits($cname)]} {
477                    array set limits $_comp2limits($cname)
478                    switch -- $which {
479                        x - xlin - xlog {
480                            set axis x
481                            foreach {axisMin axisMax} $limits(x) break
482                        }
483                        y - ylin - ylog {
484                            set axis y
485                            foreach {axisMin axisMax} $limits(y) break
486                        }
487                        z - zlin - zlog {
488                            set axis z
489                            foreach {axisMin axisMax} $limits(z) break
490                        }
491                        v - vlin - vlog {
492                            set axis v
493                            foreach {axisMin axisMax} $limits(v) break
494                        }
495                        default {
496                            if { ![info exists limits($which)] } {
497                                error "limits: unknown axis \"$which\""
498                            }
499                            set axis v
500                            foreach {axisMin axisMax} $limits($which) break
501                        }
502                    }
503                } else {
504                    set axisMin 0  ;# HACK ALERT! must be OpenDX data
505                    set axisMax 1
506                    set axis v
507                }
508            }
509        }
510        if { "" == $min || $axisMin < $min } {
511            set min $axisMin
512        }
513        if { "" == $max || $axisMax > $max } {
514            set max $axisMax
515        }
516    }
517    blt::vector destroy tmp zero
518    set val [$_field get "${axis}axis.min"]
519    if {"" != $val && "" != $min} {
520        if {$val > $min} {
521            # tool specified this min -- don't go any lower
522            set min $val
523        }
524    }
525    set val [$_field get "${axis}axis.max"]
526    if {"" != $val && "" != $max} {
527        if {$val < $max} {
528            # tool specified this max -- don't go any higher
529            set max $val
530        }
531    }
532    return [list $min $max]
533}
534
535
536# ----------------------------------------------------------------------
537# USAGE: fieldlimits
538#
539# Returns a list {min max} representing the limits for the specified
540# axis.
541# ----------------------------------------------------------------------
542itcl::body Rappture::Field::fieldlimits {} {
543    foreach cname [array names _comp2limits] {
544        array set limits $_comp2limits($cname)
545        foreach fname $_comp2fldName($cname) {
546            if { ![info exists limits($fname)] } {
547                puts stderr "ERROR: field \"$fname\" unknown in \"$cname\""
548                continue
549            }
550            foreach {min max} $limits($fname) break
551            if { ![info exists overall($fname)] } {
552                set overall($fname) $limits($fname)
553                continue
554            }
555            foreach {omin omax} $overall($fname) break
556            if { $min < $omin } {
557                set omin $min
558            }
559            if { $max > $omax } {
560                set omax $max
561            }
562            set overall($fname) [list $min $max]
563        }
564    }
565    if { [info exists overall] } {
566        return [array get overall]
567    }
568    return ""
569}
570 
571# ----------------------------------------------------------------------
572# USAGE: controls get ?<name>?
573# USAGE: controls validate <path> <value>
574# USAGE: controls put <path> <value>
575#
576# Returns a list {path1 x1 y1 val1  path2 x2 y2 val2 ...} representing
577# control points for the specified field component <name>.
578# ----------------------------------------------------------------------
579itcl::body Rappture::Field::controls {option args} {
580    switch -- $option {
581        get {
582            set cname [lindex $args 0]
583            if {[info exists _comp2cntls($cname)]} {
584                return $_comp2cntls($cname)
585            }
586            return ""
587        }
588        validate {
589            set path [lindex $args 0]
590            set value [lindex $args 1]
591            set units [$_xmlobj get $path.units]
592
593            if {"" != $units} {
594                set nv [Rappture::Units::convert \
595                    $value -context $units -to $units -units off]
596            } else {
597                set nv $value
598            }
599            if {![string is double $nv]
600                  || [regexp -nocase {^(inf|nan)$} $nv]} {
601                error "Value out of range"
602            }
603
604            set rawmin [$_xmlobj get $path.min]
605            if {"" != $rawmin} {
606                set minv $rawmin
607                if {"" != $units} {
608                    set minv [Rappture::Units::convert \
609                        $minv -context $units -to $units -units off]
610                    set nv [Rappture::Units::convert \
611                        $value -context $units -to $units -units off]
612                }
613                # fix for the case when the user tries to
614                # compare values like minv=-500 nv=-0600
615                set nv [format "%g" $nv]
616                set minv [format "%g" $minv]
617
618                if {$nv < $minv} {
619                    error "Minimum value allowed here is $rawmin"
620                }
621            }
622
623            set rawmax [$_xmlobj get $path.max]
624            if {"" != $rawmax} {
625                set maxv $rawmax
626                if {"" != $units} {
627                    set maxv [Rappture::Units::convert \
628                        $maxv -context $units -to $units -units off]
629                    set nv [Rappture::Units::convert \
630                        $value -context $units -to $units -units off]
631                }
632                # fix for the case when the user tries to
633                # compare values like maxv=-500 nv=-0600
634                set nv [format "%g" $nv]
635                set maxv [format "%g" $maxv]
636
637                if {$nv > $maxv} {
638                    error "Maximum value allowed here is $rawmax"
639                }
640            }
641
642            return "ok"
643        }
644        put {
645            set path [lindex $args 0]
646            set value [lindex $args 1]
647            $_xmlobj put $path.current $value
648            Build
649        }
650        default {
651            error "bad field controls option \"$option\": should be get or put"
652        }
653    }
654}
655
656# ----------------------------------------------------------------------
657# USAGE: hints ?<keyword>?
658#
659# Returns a list of key/value pairs for various hints about plotting
660# this field.  If a particular <keyword> is specified, then it returns
661# the hint for that <keyword>, if it exists.
662# ----------------------------------------------------------------------
663itcl::body Rappture::Field::hints {{keyword ""}} {
664    if { $keyword != "" } {
665        if {[info exists _hints($keyword)]} {
666            return $_hints($keyword)
667        }
668        return ""
669    }
670    return [array get _hints]
671}
672
673
674# ----------------------------------------------------------------------
675# USAGE: InitHints
676#
677# Returns a list of key/value pairs for various hints about plotting
678# this field.  If a particular <keyword> is specified, then it returns
679# the hint for that <keyword>, if it exists.
680# ----------------------------------------------------------------------
681itcl::body Rappture::Field::InitHints {} {
682    foreach {key path} {
683        camera          camera.position
684        color           about.color
685        default         about.default
686        group           about.group
687        label           about.label
688        scale           about.scale
689        seeds           about.seeds
690        style           about.style
691        type            about.type
692        xlabel          about.xaxis.label
693        ylabel          about.yaxis.label
694        zlabel          about.zaxis.label
695        xunits          about.xaxis.units
696        yunits          about.yaxis.units
697        zunits          about.zaxis.units
698        units           units
699        updir           updir
700        vectors         about.vectors
701    } {
702        set str [$_field get $path]
703        if { "" != $str } {
704            set _hints($key) $str
705        }
706    }
707    foreach cname [components] {
708        if { ![info exists _comp2mesh($cname)] } {
709            continue
710        }
711        set mesh [lindex $_comp2mesh($cname) 0]
712        foreach axis {x y z} {
713            if { ![info exists _hints(${axis}units)] } {
714                set _hints(${axis}units) [$mesh units $axis]
715            }
716            if { ![info exists _hints(${axis}label)] } {
717                set _hints(${axis}label) [$mesh label $axis]
718            }
719        }
720    }
721    foreach {key path} {
722        toolid          tool.id
723        toolname        tool.name
724        toolcommand     tool.execute
725        tooltitle       tool.title
726        toolrevision    tool.version.application.revision
727    } {
728        set str [$_xmlobj get $path]
729        if { "" != $str } {
730            set _hints($key) $str
731        }
732    }
733    # Set toolip and path hints
734    set _hints(path) $_path
735    if { [info exists _hints(group)] && [info exists _hints(label)] } {
736        # pop-up help for each curve
737        set _hints(tooltip) $_hints(label)
738    }
739}
740
741# ----------------------------------------------------------------------
742# USAGE: Build
743#
744# Used internally to build up the vector representation for the
745# field when the object is first constructed, or whenever the field
746# data changes.  Discards any existing vectors and builds everything
747# from scratch.
748# ----------------------------------------------------------------------
749itcl::body Rappture::Field::Build {} {
750
751    # Discard any existing data
752    foreach name [array names _comp2xy] {
753        eval blt::vector destroy $_comp2xy($name)
754    }
755    foreach name [array names _comp2dx] {
756        set vector [lindex $_comp2dx($name) 1]
757        eval blt::vector destroy $vector
758    }
759    array unset _comp2vtk
760    foreach name [array names _comp2unirect2d] {
761        eval itcl::delete object $_comp2unirect2d($name)
762    }
763    foreach name [array names _comp2unirect3d] {
764        eval itcl::delete object $_comp2unirect3d($name)
765    }
766    catch {unset _comp2xy}
767    catch {unset _comp2dx}
768    catch {unset _comp2dims}
769    catch {unset _comp2style}
770    array unset _comp2unirect2d
771    array unset _comp2unirect3d
772    array unset _comp2extents
773    array unset _dataobj2type
774    #
775    # Scan through the components of the field and create
776    # vectors for each part.
777    #
778    array unset _isValidComponent
779    foreach cname [$_field children -type component] {
780        set type ""
781        if { ([$_field element $cname.constant] != "" &&
782              [$_field element $cname.domain] != "") ||
783              [$_field element $cname.xy] != "" } {
784            set type "1D"
785        } elseif { [$_field element $cname.mesh] != "" &&
786                   [$_field element $cname.values] != ""} {
787            set type "points-on-mesh"
788        } elseif { [$_field element $cname.vtk] != ""} {
789            set type "vtk"
790            set viewer [$_field get "about.view"]
791            if { $viewer != "" } {
792                set _viewer $viewer
793            }
794        } elseif {[$_field element $cname.opendx] != ""} {
795            global env
796            if { [info exists env(VTKVOLUME)] } {
797                set _viewer "vtkvolume"
798            }
799            set type "opendx"
800        } elseif {[$_field element $cname.dx] != ""} {
801            global env
802            if { [info exists env(VTKVOLUME)] } {
803                set _viewer "vtkvolume"
804            }
805            set type "dx"
806        } elseif {[$_field element $cname.ucd] != ""} {
807            set type "ucd"
808        } elseif {[$_field element $cname.dicom] != ""} {
809            set type "dicom"
810        }
811        set _comp2style($cname) ""
812        if { $type == "" } {
813            puts stderr "WARNING: ignoring field component \"$_path.$cname\": no data found."
814            continue
815        }
816        # Save the extents of the component
817        if { [$_field element $cname.extents] != "" } {
818            set extents [$_field get $cname.extents]
819        } else {
820            set extents 1
821        }
822        set _comp2extents($cname) $extents
823        set _type $type
824
825        GetTypeAndSize $cname
826        GetAssociation $cname
827        if { $_comp2size($cname) > 1 } {
828            set viewer [$_field get "about.view"]
829            if { $viewer == "" } {
830                set _viewer "glyphs"
831            }
832        }
833        if {$type == "1D"} {
834            #
835            # 1D data can be represented as 2 BLT vectors,
836            # one for x and the other for y.
837            #
838            set xv ""
839            set yv ""
840            set _dim 1
841            set val [$_field get $cname.constant]
842            if {$val != ""} {
843                set domain [$_field get $cname.domain]
844                if {$domain == "" || ![info exists _limits($domain)]} {
845                    set z0 0
846                    set z1 $_zmax
847                } else {
848                    foreach {z0 z1} $_limits($domain) { break }
849                }
850                set xv [blt::vector create x$_counter]
851                $xv append $z0 $z1
852
853                foreach {val pcomp} [_getValue $val] break
854                set yv [blt::vector create y$_counter]
855                $yv append $val $val
856
857                if {$pcomp != ""} {
858                    set zm [expr {0.5*($z0+$z1)}]
859                    set _comp2cntls($cname) \
860                        [list $pcomp $zm $val "$val$_units"]
861                }
862            } else {
863                set xydata [$_field get $cname.xy]
864                if {"" != $xydata} {
865                    set xv [blt::vector create x$_counter]
866                    set yv [blt::vector create y$_counter]
867                    set tmp [blt::vector create \#auto]
868                    $tmp set $xydata
869                    $tmp split $xv $yv
870                    blt::vector destroy $tmp
871                }
872            }
873
874            if {$xv != "" && $yv != ""} {
875                # sort x-coords in increasing order
876                $xv sort $yv
877                set _dim 1
878                set _comp2dims($cname) "1D"
879                set _comp2xy($cname) [list $xv $yv]
880                incr _counter
881            }
882        } elseif {$type == "points-on-mesh"} {
883            if { ![BuildPointsOnMesh $cname] } {
884                continue;               # Ignore this component
885            }
886        } elseif {$type == "vtk"} {
887            set contents [$_field get $cname.vtk]
888            if { $contents == "" } {
889                puts stderr "WARNING: no data fo \"$_path.$cname.vtk\""
890                continue;               # Ignore this component
891            }
892            ReadVtkDataSet $cname $contents
893            set _comp2vtk($cname) $contents
894            set _comp2style($cname) [$_field get $cname.style]
895            incr _counter
896        } elseif {$type == "dx" || $type == "opendx" } {
897            #
898            # HACK ALERT!  Extract gzipped, base64-encoded OpenDX
899            # data.  Assume that it's 3D.  Pass it straight
900            # off to the NanoVis visualizer.
901            #
902            set viewer [$_field get "about.view"]
903            if { $viewer != "" } {
904                set _viewer $viewer
905            }
906            if { $_viewer == "" } {
907                set _viewer "nanovis"
908            }
909            set _dim 3
910            set _comp2dims($cname) "3D"
911            set data [$_field get -decode no $cname.$type]
912            set contents [Rappture::encoding::decode -as zb64 $data]
913            if { $contents == "" } {
914                puts stderr "WARNING: no data for \"$_path.$cname.$type\""
915                continue;               # Ignore this component
916            }
917            set vector ""
918            if 0 {
919                set f [open /tmp/$_path.$cname.dx "w"]
920                puts -nonewline $f $contents
921                close $f
922            }
923            # This is temporary.  I put a check for this in the DxToVtk
924            # parser. 
925            if { [string range $contents  0 3] == "<DX>" } {
926                set contents [string range $contents 4 end]
927            }
928            if { [catch { Rappture::DxToVtk $contents } vtkdata] == 0 } {
929                set vector [ReadVtkDataSet $cname $vtkdata]
930            } else {
931                puts stderr "Can't parse dx data: $vtkdata"
932            }
933            if 0 {
934                set f [open /tmp/$_path.$cname.vtk "w"]
935                puts -nonewline $f $vtkdata
936                close $f
937            }
938            if { $_viewer != "nanovis" && $_viewer != "flowvis" } {
939                set _type "vtk"
940                set _comp2vtk($cname) $vtkdata
941            } else {
942                set _type "dx"
943                set _comp2dx($cname) [list $data $vector]
944            }
945            set _comp2style($cname) [$_field get $cname.style]
946            if {[$_field element $cname.flow] != ""} {
947                set _comp2flowhints($cname) \
948                    [Rappture::FlowHints ::\#auto $_field $cname $_units]
949            }
950            set _dim 3
951            incr _counter
952        } elseif { $type == "dicom"} {
953            set contents [$_field get $cname.dicom]
954            if { $contents == "" } {
955                continue;               # Ignore this compoennt
956            }
957            if { ![file isdir $contents] } {
958                puts stderr "dicom path \"$contents\" is not a directory."
959                return 0
960            }
961            set vtkdata [DicomSeriesToVtk $cname $contents]
962            ReadVtkDataSet $cname $vtkdata
963            set _comp2vtk($cname) $vtkdata
964            set _comp2style($cname) [$_field get $cname.style]
965            incr _counter
966            if { $_viewer == "" } {
967                set _viewer "nanovis"
968            }
969        } elseif { $type == "ucd"} {
970            set contents [$_field get $cname.ucd]
971            if { $contents == "" } {
972                continue;               # Ignore this compoennt
973            }
974            set vtkdata [AvsToVtk $cname $contents]
975            ReadVtkDataSet $cname $vtkdata
976            set _comp2vtk($cname) $vtkdata
977            set _comp2style($cname) [$_field get $cname.style]
978            incr _counter
979        }
980        set _isValidComponent($cname) 1
981    }
982    if { [array size _isValidComponent] == 0 } {
983        puts stderr "WARNING: no valid components for field \"$_path\""
984        return 0
985    }
986    # Sanity check.  Verify that all components of the field have the same
987    # dimension.
988    set dim ""
989    foreach cname [array names _comp2dims] {
990        if { $dim == "" } {
991            set dim $_comp2dims($cname)
992            continue
993        }
994        if { $dim != $_comp2dims($cname) } {
995            puts stderr "WARNING: field can't have components of different dimensions: [join [array get _comp2dims] ,]"
996            return 0
997        }
998    }
999    # FIXME: about.scalars and about.vectors are temporary.  With views
1000    #        the label and units for each field will be specified there.
1001    #
1002    # FIXME: Test that every <field><component> has the same field names,
1003    #        units, components.
1004    #
1005    # Override what we found in the VTK file with names that the user
1006    # selected.  We override the field label and units.
1007    foreach { fname label units } [$_field get about.scalars] {
1008        if { ![info exists _fld2Name($fname)] } {
1009            set _fld2Name($fname) $fname
1010            set _fld2Components($fname) 1
1011        }
1012        set _fld2Label($fname) $label
1013        set _fld2Units($fname) $units
1014    }
1015    foreach { fname label units } [$_field get about.vectors] {
1016        if { ![info exists _fld2Name($fname)] } {
1017            set _fld2Name($fname) $fname
1018            # We're just marking the field as vector (> 1) for now.
1019            set _fld2Components($fname) 3
1020        }
1021        set _fld2Label($fname) $label
1022        set _fld2Units($fname) $units
1023    }
1024    set _isValid 1
1025    return 1
1026}
1027
1028# ----------------------------------------------------------------------
1029# USAGE: _getValue <expr>
1030#
1031# Used internally to get the value for an expression <expr>.  Returns
1032# a list of the form {val parameterPath}, where val is the numeric
1033# value of the expression, and parameterPath is the XML path to the
1034# parameter representing the value, or "" if the <expr> does not
1035# depend on any parameters.
1036# ----------------------------------------------------------------------
1037itcl::body Rappture::Field::_getValue {expr} {
1038    #
1039    # First, look for the expression among the <parameter>'s
1040    # associated with the device.
1041    #
1042    set found 0
1043    foreach pcomp [$_xmlobj children parameters] {
1044        set id [$_xmlobj element -as id parameters.$pcomp]
1045        if {[string equal $id $expr]} {
1046            set val [$_xmlobj get parameters.$pcomp.current]
1047            if {"" == $val} {
1048                set val [$_xmlobj get parameters.$pcomp.default]
1049            }
1050            if {"" != $val} {
1051                set expr $val
1052                set found 1
1053                break
1054            }
1055        }
1056    }
1057    if {$found} {
1058        set pcomp "parameters.$pcomp"
1059    } else {
1060        set pcomp ""
1061    }
1062
1063    if {$_units != ""} {
1064        set expr [Rappture::Units::convert $expr \
1065            -context $_units -to $_units -units off]
1066    }
1067
1068    return [list $expr $pcomp]
1069}
1070
1071#
1072# isunirect2d  --
1073#
1074# Returns if the field is a unirect2d object. 
1075#
1076itcl::body Rappture::Field::isunirect2d { } {
1077    return [expr [array size _comp2unirect2d] > 0]
1078}
1079
1080#
1081# isunirect3d  --
1082#
1083# Returns if the field is a unirect3d object. 
1084#
1085itcl::body Rappture::Field::isunirect3d { } {
1086    return [expr [array size _comp2unirect3d] > 0]
1087}
1088
1089#
1090# flowhints  --
1091#
1092# Returns the hints associated with a flow vector field. 
1093#
1094itcl::body Rappture::Field::flowhints { cname } {
1095    if { [info exists _comp2flowhints($cname)] } {
1096        return $_comp2flowhints($cname)
1097    }
1098    return ""
1099}
1100
1101#
1102# style  --
1103#
1104# Returns the style associated with a component of the field. 
1105#
1106itcl::body Rappture::Field::style { cname } {
1107    if { [info exists _comp2style($cname)] } {
1108        return $_comp2style($cname)
1109    }
1110    return ""
1111}
1112
1113#
1114# type  --
1115#
1116# Returns the data storage type of the field.
1117#
1118# FIXME: What are the valid types?
1119#
1120itcl::body Rappture::Field::type {} {
1121    return $_type
1122}
1123
1124#
1125# numComponents --
1126#
1127# Returns if the number of components in the field component.
1128#
1129itcl::body Rappture::Field::numComponents {cname} {
1130    set name $cname
1131    regsub -all { } $name {_} name
1132    if {[info exists _fld2Components($name)] } {
1133        return $_fld2Components($name)
1134    }
1135    return 1;                           # Default to scalar.
1136}
1137
1138#
1139# extents --
1140#
1141# Returns if the field is a unirect2d object. 
1142#
1143itcl::body Rappture::Field::extents {{cname -overall}} {
1144    if {$cname == "-overall" } {
1145        set max 0
1146        foreach cname [$_field children -type component] {
1147            if { ![info exists _comp2unirect3d($cname)] &&
1148                 ![info exists _comp2extents($cname)] } {
1149                continue
1150            }
1151            set value $_comp2extents($cname)
1152            if { $max < $value } {
1153                set max $value
1154            }
1155        }
1156        return $max
1157    }
1158    if { $cname == "component0"} {
1159        set cname [lindex [components -name] 0]
1160    }
1161    return $_comp2extents($cname)
1162}
1163
1164itcl::body Rappture::Field::ConvertToVtkData { cname } {
1165    set ds ""
1166    switch -- [typeof $cname] {
1167        "unirect2d" {
1168            foreach { x1 x2 xN y1 y2 yN } [$dataobj mesh $cname] break
1169            set spacingX [expr {double($x2 - $x1)/double($xN - 1)}]
1170            set spacingY [expr {double($y2 - $y1)/double($yN - 1)}]
1171           
1172            set ds [vtkImageData $this-grdataTemp]
1173            $ds SetDimensions $xN $yN 1
1174            $ds SetOrigin $x1 $y1 0
1175            $ds SetSpacing $spacingX $spacingY 0
1176            set arr [vtkDoubleArray $this-arrTemp]
1177            foreach {val} [$dataobj values $cname] {
1178                $arr InsertNextValue $val
1179            }
1180            [$ds GetPointData] SetScalars $arr
1181        }
1182        "unirect3d" {
1183            foreach { x1 x2 xN y1 y2 yN z1 z2 zN } [$dataobj mesh $cname] break
1184            set spacingX [expr {double($x2 - $x1)/double($xN - 1)}]
1185            set spacingY [expr {double($y2 - $y1)/double($yN - 1)}]
1186            set spacingZ [expr {double($z2 - $z1)/double($zN - 1)}]
1187           
1188            set ds [vtkImageData $this-grdataTemp]
1189            $ds SetDimensions $xN $yN $zN
1190            $ds SetOrigin $x1 $y1 $z1
1191            $ds SetSpacing $spacingX $spacingY $spacingZ
1192            set arr [vtkDoubleArray $this-arrTemp]
1193            foreach {val} [$dataobj values $cname] {
1194                $arr InsertNextValue $val
1195            }
1196            [$ds GetPointData] SetScalars $val
1197        }
1198        "contour" {
1199            return [$dataobj blob $cname]
1200        }
1201        "dx" {
1202            return [Rappture::DxToVtk $_comp2dx($cname)]
1203        }
1204        default {
1205            set mesh [$dataobj mesh $cname]
1206            switch -- [$mesh GetClassName] {
1207                vtkPoints {
1208                    # handle cloud of points
1209                    set ds [vtkPolyData $this-polydataTemp]
1210                    $ds SetPoints $mesh
1211                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1212                }
1213                vtkPolyData {
1214                    set ds [vtkPolyData $this-polydataTemp]
1215                    $ds ShallowCopy $mesh
1216                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1217                }
1218                vtkUnstructuredGrid {
1219                    # handle 3D grid with connectivity
1220                    set ds [vtkUnstructuredGrid $this-grdataTemp]
1221                    $ds ShallowCopy $mesh
1222                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1223                }
1224                vtkRectilinearGrid {
1225                    # handle 3D grid with connectivity
1226                    set ds [vtkRectilinearGrid $this-grdataTemp]
1227                    $ds ShallowCopy $mesh
1228                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1229                }
1230                default {
1231                    error "don't know how to handle [$mesh GetClassName] data"
1232                }
1233            }
1234        }
1235    }
1236
1237    if {"" != $ds} {
1238        set writer [vtkDataSetWriter $this-dsWriterTmp]
1239        $writer SetInput $ds
1240        $writer SetFileTypeToASCII
1241        $writer WriteToOutputStringOn
1242        $writer Write
1243        set out [$writer GetOutputString]
1244        $ds Delete
1245        $writer Delete
1246    } else {
1247        set out ""
1248        error "No DataSet to write"
1249    }
1250
1251    append out "\n"
1252    return $out
1253}
1254
1255itcl::body Rappture::Field::VerifyVtkDataSet { contents } {
1256    package require vtk
1257
1258    set reader $this-datasetreader
1259    vtkDataSetReader $reader
1260
1261    # Write the contents to a file just in case it's binary.
1262    set tmpfile file[pid].vtk
1263    set f [open "$tmpfile" "w"]
1264    fconfigure $f -translation binary -encoding binary
1265    puts $f $contents
1266    close $f
1267
1268    $reader SetFileName $tmpfile
1269    $reader ReadAllNormalsOn
1270    $reader ReadAllTCoordsOn
1271    $reader ReadAllScalarsOn
1272    $reader ReadAllColorScalarsOn
1273    $reader ReadAllVectorsOn
1274    $reader ReadAllTensorsOn
1275    $reader ReadAllFieldsOn
1276    $reader Update
1277    file delete $tmpfile
1278
1279    set dataset [$reader GetOutput]
1280    set dataAttrs [$dataset GetPointData]
1281    if { $dataAttrs == ""} {
1282        puts stderr "WARNING: no point data found in \"$_path\""
1283        rename $reader ""
1284        return 0
1285    }
1286    rename $reader ""
1287}
1288
1289itcl::body Rappture::Field::ReadVtkDataSet { cname contents } {
1290    package require vtk
1291
1292    set reader $this-datasetreader
1293    vtkDataSetReader $reader
1294
1295    # Write the contents to a file just in case it's binary.
1296    set tmpfile file[pid].vtk
1297    set f [open "$tmpfile" "w"]
1298    fconfigure $f -translation binary -encoding binary
1299    puts $f $contents
1300    close $f
1301
1302    $reader SetFileName $tmpfile
1303    $reader ReadAllNormalsOn
1304    $reader ReadAllTCoordsOn
1305    $reader ReadAllScalarsOn
1306    $reader ReadAllColorScalarsOn
1307    $reader ReadAllVectorsOn
1308    $reader ReadAllTensorsOn
1309    $reader ReadAllFieldsOn
1310    $reader Update
1311    file delete $tmpfile
1312
1313    set dataset [$reader GetOutput]
1314    set limits {}
1315    foreach {xmin xmax ymin ymax zmin zmax} [$dataset GetBounds] break
1316    # Figure out the dimension of the mesh from the bounds.
1317    set _dim 0
1318    if { $xmax > $xmin } {
1319        incr _dim
1320    }
1321    if { $ymax > $ymin } {
1322        incr _dim
1323    }
1324    if { $zmax > $zmin } {
1325        incr _dim
1326    }
1327    if { $_viewer == "" } {
1328        if { $_dim == 2 } {
1329            set _viewer contour
1330        } else {
1331            set _viewer isosurface
1332        }
1333    }
1334    set _comp2dims($cname) ${_dim}D
1335    if { $_dim < 2 } {
1336        set numPoints [$dataset GetNumberOfPoints]
1337        set xv [blt::vector create \#auto]
1338        for { set i 0 } { $i < $numPoints } { incr i } {
1339            set point [$dataset GetPoint $i]
1340            $xv append [lindex $point 0]
1341        }
1342        set yv [blt::vector create \#auto]
1343        set dataAttrs [$dataset GetPointData]
1344        if { $dataAttrs == ""} {
1345            puts stderr "WARNING: no point data found in \"$_path\""
1346            rename $reader ""
1347            return 0
1348        }
1349        set array [$dataAttrs GetScalars]
1350        if { $array == ""} {
1351            puts stderr "WARNING: no scalar point data found in \"$_path\""
1352            rename $reader ""
1353            return 0
1354        }
1355        set numTuples [$array GetNumberOfTuples]
1356        for { set i 0 } { $i < $numTuples } { incr i } {
1357            $yv append [$array GetComponent $i 0]
1358        }
1359        $xv sort $yv
1360        set _comp2xy($cname) [list $xv $yv]
1361    }
1362    lappend limits x [list $xmin $xmax]
1363    lappend limits y [list $ymin $ymax]
1364    lappend limits z [list $zmin $zmax]
1365    set dataAttrs [$dataset GetPointData]
1366    if { $dataAttrs == ""} {
1367        puts stderr "WARNING: no point data found in \"$_path\""
1368        rename $reader ""
1369        return 0
1370    }
1371    set vmin 0
1372    set vmax 1
1373    set numArrays [$dataAttrs GetNumberOfArrays]
1374    set vector [blt::vector create \#auto]
1375    if { $numArrays > 0 } {
1376        set array [$dataAttrs GetArray 0]
1377        # Calling GetRange with component set to -1 will return
1378        # either the scalar range or vector magnitude range
1379        foreach {vmin vmax} [$array GetRange -1] break
1380
1381        set numTuples [$array GetNumberOfTuples]
1382        for { set i 0 } { $i < $numTuples } { incr i } {
1383            $vector append [$array GetComponent $i 0]
1384        }
1385        for {set i 0} {$i < [$dataAttrs GetNumberOfArrays] } {incr i} {
1386            set array [$dataAttrs GetArray $i]
1387            set fname  [$dataAttrs GetArrayName $i]
1388            foreach {min max} [$array GetRange -1] break
1389            lappend limits $fname [list $min $max]
1390            set _fld2Units($fname) ""
1391            set _fld2Label($fname) $fname
1392            # Let the VTK file override the <type> designated.
1393            set _fld2Components($fname) [$array GetNumberOfComponents]
1394            lappend _comp2fldName($cname) $fname
1395        }
1396    }
1397   
1398    lappend limits v [list $vmin $vmax]
1399    set _comp2limits($cname) $limits
1400
1401    rename $reader ""
1402    return $vector
1403}
1404
1405#
1406# vtkdata --
1407#
1408#       Returns a string representing the mesh and field data for a specific
1409#       component in the legacy VTK file format.
1410#
1411itcl::body Rappture::Field::vtkdata {cname} {
1412    if {$cname == "component0"} {
1413        set cname "component"
1414    }
1415    # DX: Convert DX to VTK
1416    if {[info exists _comp2dx($cname)]} {
1417        set data $_comp2dx($cname)
1418        set data [Rappture::encoding::decode $data]
1419        return [Rappture::DxToVtk $data]
1420    }
1421    # Unirect3d: isosurface
1422    if {[info exists _comp2unirect3d($cname)]} {
1423        return [$_comp2unirect3d($cname) vtkdata]
1424    }
1425    # VTK file data:
1426    if { [info exists _comp2vtk($cname)] } {
1427        return $_comp2vtk($cname)
1428    }
1429    # Points on mesh:  Construct VTK file output.
1430    if { [info exists _comp2mesh($cname)] } {
1431        # Data is in the form mesh and vector
1432        foreach {mesh vector} $_comp2mesh($cname) break
1433        set label $cname
1434        regsub -all { } $label {_} label
1435        append out "# vtk DataFile Version 3.0\n"
1436        append out "[hints label]\n"
1437        append out "ASCII\n"
1438        append out [$mesh vtkdata]
1439
1440        if { $_comp2assoc($cname) == "pointdata" } {
1441            set vtkassoc "POINT_DATA"
1442        } elseif { $_comp2assoc($cname) == "celldata" } {
1443            set vtkassoc "CELL_DATA"
1444        } elseif { $_comp2assoc($cname) == "fielddata" } {
1445            set vtkassoc "FIELD"
1446        } else {
1447            error "unknown association \"$_comp2assoc($cname)\""
1448        }
1449        set elemSize [numComponents $cname]
1450        set numValues [expr [$vector length] / $elemSize]
1451        if { $_comp2assoc($cname) == "fielddata" } {
1452            append out "$vtkassoc FieldData 1\n"
1453            append out "$label $elemSize $numValues double\n"
1454        } else {
1455            append out "$vtkassoc $numValues\n"
1456            if { $_comp2type($cname) == "colorscalars" } {
1457                # Must be float for ASCII, unsigned char for BINARY
1458                append out "COLOR_SCALARS $label $elemSize\n"
1459            } elseif { $_comp2type($cname) == "normals" } {
1460                # elemSize must equal 3
1461                append out "NORMALS $label double\n"
1462            } elseif { $_comp2type($cname) == "scalars" } {
1463                # elemSize can be 1, 2, 3 or 4
1464                append out "SCALARS $label double $elemSize\n"
1465                append out "LOOKUP_TABLE default\n"
1466            } elseif { $_comp2type($cname) == "tcoords" } {
1467                # elemSize must be 1, 2, or 3
1468                append out "TEXTURE_COORDINATES $label $elemSize double\n"
1469            } elseif { $_comp2type($cname) == "tensors" } {
1470                # elemSize must equal 9
1471                append out "TENSORS $label double\n"
1472            } elseif { $_comp2type($cname) == "vectors" } {
1473                # elemSize must equal 3
1474                append out "VECTORS $label double\n"
1475            } else {
1476                error "unknown element type \"$_comp2type($cname)\""
1477            }
1478        }
1479        append out [$vector range 0 end]
1480        append out "\n"
1481        if 0 {
1482            VerifyVtkDataSet $out
1483        }
1484        return $out
1485    }
1486    error "can't find vtkdata for $cname. This method should only be called by the vtkheightmap widget"
1487}
1488
1489#
1490# BuildPointsOnMesh --
1491#
1492#       Parses the field XML description to build a mesh and values vector
1493#       representing the field.  Right now we handle the deprecated types
1494#       of "cloud", "unirect2d", and "unirect3d" (mostly for flows).
1495#
1496itcl::body Rappture::Field::BuildPointsOnMesh {cname} {
1497    #
1498    # More complex 2D/3D data is represented by a mesh
1499    # object and an associated vector for field values.
1500    #
1501    set path [$_field get $cname.mesh]
1502    if {[$_xmlobj element $path] == ""} {
1503        # Unknown mesh designated.
1504        return 0
1505    }
1506    set viewer [$_field get "about.view"]
1507    if { $viewer != "" } {
1508        set _viewer $viewer
1509    }
1510    set element [$_xmlobj element -as type $path]
1511    set name $cname
1512    regsub -all { } $name {_} name
1513    set _fld2Label($name) $name
1514    set label [hints zlabel]
1515    if { $label != "" } {
1516        set _fld2Label($name) $label
1517    }
1518    set _fld2Units($name) [hints zunits]
1519    set _fld2Components($name) $_comp2size($cname)
1520    lappend _comp2fldName($cname) $name
1521
1522    # Handle bizarre cases that hopefully will be deprecated.
1523    if { $element == "unirect3d" } {
1524        # Special case: unirect3d (should be deprecated) + flow.
1525        if { [$_field element $cname.extents] != "" } {
1526            set vectorsize [$_field get $cname.extents]
1527        } else {
1528            set vectorsize 1
1529        }
1530        set _type unirect3d
1531        set _dim 3
1532        if { $_viewer == "" } {
1533            set _viewer flowvis
1534        }
1535        set _comp2dims($cname) "3D"
1536        set _comp2unirect3d($cname) \
1537            [Rappture::Unirect3d \#auto $_xmlobj $_field $cname $vectorsize]
1538        set _comp2style($cname) [$_field get $cname.style]
1539        set limits {}
1540        foreach axis { x y z } {
1541            lappend limits $axis [$_comp2unirect3d($cname) limits $axis]
1542        }
1543        # Get the data limits
1544        set vector [$_comp2unirect3d($cname) valuesObj]
1545        lappend limits $cname [$vector limits]
1546        lappend limits v      [$vector limits]
1547        set _comp2limits($cname) $limits
1548        if {[$_field element $cname.flow] != ""} {
1549            set _comp2flowhints($cname) \
1550                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1551        }
1552        incr _counter
1553        return 1
1554    }
1555    if { $element == "unirect2d" && [$_field element $cname.flow] != "" } {
1556        # Special case: unirect2d (normally deprecated) + flow.
1557        if { [$_field element $cname.extents] != "" } {
1558            set vectorsize [$_field get $cname.extents]
1559        } else {
1560            set vectorsize 1
1561        }
1562        set _type unirect2d
1563        set _dim 2
1564        if { $_viewer == "" } {
1565            set _viewer "flowvis"
1566        }
1567        set _comp2dims($cname) "2D"
1568        set _comp2unirect2d($cname) \
1569            [Rappture::Unirect2d \#auto $_xmlobj $path]
1570        set _comp2style($cname) [$_field get $cname.style]
1571        set _comp2flowhints($cname) \
1572            [Rappture::FlowHints ::\#auto $_field $cname $_units]
1573        set _values [$_field get $cname.values]
1574        set limits {}
1575        foreach axis { x y } {
1576            lappend limits $axis [$_comp2unirect2d($cname) limits $axis]
1577        }
1578        set xv [blt::vector create \#auto]
1579        $xv set $_values
1580        lappend limits $cname [$xv limits]
1581        lappend limits v [$xv limits]
1582        lappend limits z [$xv limits]
1583        blt::vector destroy $xv
1584        set _comp2limits($cname) $limits
1585        incr _counter
1586        return 1
1587    }
1588    switch -- $element {
1589        "cloud" {
1590            set mesh [Rappture::Cloud::fetch $_xmlobj $path]
1591            set _type cloud
1592        }
1593        "mesh" {
1594            set mesh [Rappture::Mesh::fetch $_xmlobj $path]
1595            set _type mesh
1596        }           
1597        "unirect2d" {
1598            if { $_viewer == "" } {
1599                set _viewer "heightmap"
1600            }
1601            set mesh [Rappture::Unirect2d::fetch $_xmlobj $path]
1602            set _type unirect2d
1603        }
1604    }
1605    if { ![$mesh isvalid] } {
1606        puts stderr "Mesh is invalid"
1607        return 0
1608    }
1609    set _dim [$mesh dimensions]
1610    if { $_dim == 3 } {
1611        set dim 0
1612        foreach axis {x y z} {
1613            foreach {min max} [$mesh limits $axis] {
1614                if { $min < $max } {
1615                    incr dim
1616                }
1617            }
1618        }
1619        if { $dim  != 3 } {
1620            set _dim $dim
1621        }
1622    }
1623       
1624    if {$_dim == 1} {
1625        # 1D data: Create vectors for graph widget.
1626        # Is this used anywhere?
1627        #
1628        # OOPS!  This is 1D data
1629        # Forget the cloud/field -- store BLT vectors
1630        #
1631        # Is there a natural growth path in generating output from 1D to
1632        # higher dimensions?  If there isn't, let's kill this in favor
1633        # or explicitly using a <curve> instead.  Otherwise, the features
1634        # (methods such as xmarkers) or the <curve> need to be added
1635        # to the <field>.
1636        #
1637        set xv [blt::vector create x$_counter]
1638        set yv [blt::vector create y$_counter]
1639       
1640        $yv set [$mesh points]
1641        $xv seq 0 1 [$yv length]
1642        # sort x-coords in increasing order
1643        $xv sort $yv
1644       
1645        set _comp2dims($cname) "1D"
1646        set _comp2xy($cname) [list $xv $yv]
1647        incr _counter
1648        return 1
1649    }
1650    if {$_dim == 2} {
1651        # 2D data: By default surface or contour plot using heightmap widget.
1652        set v [blt::vector create \#auto]
1653        $v set [$_field get $cname.values]
1654        if { [$v length] == 0 } {
1655            return 0
1656        }
1657        if { $_viewer == "" } {
1658            set _viewer "contour"
1659        }
1660        set _comp2dims($cname) "[$mesh dimensions]D"
1661        set _comp2mesh($cname) [list $mesh $v]
1662        set _comp2style($cname) [$_field get $cname.style]
1663        if {[$_field element $cname.flow] != ""} {
1664            set _comp2flowhints($cname) \
1665                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1666        }
1667        incr _counter
1668        array unset _comp2limits $cname
1669        lappend _comp2limits($cname) x [$mesh limits x]
1670        lappend _comp2limits($cname) y [$mesh limits y]
1671        lappend _comp2limits($cname) $cname [$v limits]
1672        lappend _comp2limits($cname) v [$v limits]
1673        return 1
1674    }
1675    if {$_dim == 3} {
1676        # 3D data: By default isosurfaces plot using isosurface widget.
1677        if { $_viewer == "" } {
1678            set _viewer "isosurface"
1679        }
1680        set v [blt::vector create \#auto]
1681        $v set [$_field get $cname.values]
1682        if { [$v length] == 0 } {
1683            return 0
1684        }
1685        set _comp2dims($cname) "[$mesh dimensions]D"
1686        set _comp2mesh($cname) [list $mesh $v]
1687        set _comp2style($cname) [$_field get $cname.style]
1688        if {[$_field element $cname.flow] != ""} {
1689            set _comp2flowhints($cname) \
1690                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1691        }
1692        incr _counter
1693        lappend _comp2limits($cname) x [$mesh limits x]
1694        lappend _comp2limits($cname) y [$mesh limits y]
1695        lappend _comp2limits($cname) z [$mesh limits z]
1696        lappend _comp2limits($cname) $cname [$v limits]
1697        lappend _comp2limits($cname) v [$v limits]
1698        return 1
1699    }
1700    error "unhandled case in field dim=$_dim element=$element"
1701}
1702
1703itcl::body Rappture::Field::AvsToVtk { cname contents } {
1704    package require vtk
1705
1706    set reader $this-datasetreader
1707    vtkAVSucdReader $reader
1708
1709    # Write the contents to a file just in case it's binary.
1710    set tmpfile $cname[pid].ucd
1711    set f [open "$tmpfile" "w"]
1712    fconfigure $f -translation binary -encoding binary
1713    puts $f $contents
1714    close $f
1715    $reader SetFileName $tmpfile
1716    $reader Update
1717    file delete $tmpfile
1718
1719    set tmpfile $cname[pid].vtk
1720    set writer $this-datasetwriter
1721    vtkDataSetWriter $writer
1722    $writer SetInputConnection [$reader GetOutputPort]
1723    $writer SetFileName $tmpfile
1724    $writer Write
1725    rename $reader ""
1726    rename $writer ""
1727
1728    set f [open "$tmpfile" "r"]
1729    fconfigure $f -translation binary -encoding binary
1730    set vtkdata [read $f]
1731    close $f
1732    file delete $tmpfile
1733    return $vtkdata
1734}
1735
1736itcl::body Rappture::Field::DicomSeriesToVtk { cname path } {
1737    package require vtk
1738
1739    set reader $this-datasetreader
1740    vtkDICOMImageReader $reader
1741
1742    set files [glob -nocomplain $path/*.dcm]
1743    if { [llength $files] == 0 } {
1744        puts stderr "no dicom files found in \"$path\""
1745        return 0
1746    }       
1747    $reader SetDirectoryName $path
1748    $reader Update
1749
1750    set tmpfile $cname[pid].vtk
1751    set writer $this-datasetwriter
1752    vtkDataSetWriter $writer
1753    $writer SetInputConnection [$reader GetOutputPort]
1754    $writer SetFileName $tmpfile
1755    $writer Write
1756    rename $reader ""
1757    rename $writer ""
1758
1759    set f [open "$tmpfile" "r"]
1760    fconfigure $f -translation binary -encoding binary
1761    set vtkdata [read $f]
1762    close $f
1763    #file delete $tmpfile
1764    return $vtkdata
1765}
1766
1767itcl::body Rappture::Field::GetTypeAndSize { cname } {
1768    array set type2components {
1769        "colorscalars"         4
1770        "normals"              3
1771        "scalars"              1
1772        "tcoords"              2
1773        "tensors"              9
1774        "vectors"              3
1775    }
1776    set type [$_field get $cname.elemtype]
1777    if { $type == "" } {
1778        set type "scalars"
1779    }
1780    if { ![info exists type2components($type)] } {
1781        error "unknown <elemtype> \"$type\" in field"
1782    }
1783    set size [$_field get $cname.elemsize]
1784    if { $size == "" } {
1785        set size $type2components($type)
1786    }
1787    set _comp2type($cname) $type
1788    set _comp2size($cname) $size
1789}
1790
1791itcl::body Rappture::Field::GetAssociation { cname } {
1792    set assoc [$_field get $cname.association]
1793    if { $assoc == "" } {
1794        set _comp2assoc($cname) "pointdata"
1795        return
1796    }
1797    switch -- $assoc {
1798        "pointdata" - "celldata" - "fielddata" {
1799            set _comp2assoc($cname) $assoc
1800            return
1801        }
1802        default {
1803            error "unknown field association \"$assoc\""
1804        }
1805    }
1806}
Note: See TracBrowser for help on using the repository browser.