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

Last change on this file since 6651 was 6537, checked in by ldelgass, 8 years ago

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