source: branches/1.4/gui/scripts/field.tcl @ 5010

Last change on this file since 5010 was 5010, checked in by ldelgass, 6 years ago

merge r4997 from trunk

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