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

Last change on this file since 4186 was 4177, checked in by ldelgass, 11 years ago

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