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

Last change on this file since 4999 was 4999, checked in by ldelgass, 9 years ago

fix typo in comment
File size: 59.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    # Unirect3d: isosurface
1313    if {[info exists _comp2unirect3d($cname)]} {
1314        return [$_comp2unirect3d($cname) vtkdata]
1315    }
1316    # VTK file data:
1317    if { [info exists _comp2vtk($cname)] } {
1318        return $_comp2vtk($cname)
1319    }
1320    # Points on mesh:  Construct VTK file output.
1321    if { [info exists _comp2mesh($cname)] } {
1322        # Data is in the form mesh and vector
1323        foreach {mesh vector} $_comp2mesh($cname) break
1324        set label $cname
1325        regsub -all { } $label {_} label
1326        append out "# vtk DataFile Version 3.0\n"
1327        append out "[hints label]\n"
1328        append out "ASCII\n"
1329        append out [$mesh vtkdata]
1330
1331        if { $_comp2assoc($cname) == "pointdata" } {
1332            set vtkassoc "POINT_DATA"
1333        } elseif { $_comp2assoc($cname) == "celldata" } {
1334            set vtkassoc "CELL_DATA"
1335        } elseif { $_comp2assoc($cname) == "fielddata" } {
1336            set vtkassoc "FIELD"
1337        } else {
1338            error "unknown association \"$_comp2assoc($cname)\""
1339        }
1340        set elemSize [numComponents $cname]
1341        set numValues [expr [$vector length] / $elemSize]
1342        if { $_comp2assoc($cname) == "fielddata" } {
1343            append out "$vtkassoc FieldData 1\n"
1344            append out "$label $elemSize $numValues double\n"
1345        } else {
1346            append out "$vtkassoc $numValues\n"
1347            if { $_comp2type($cname) == "colorscalars" } {
1348                # Must be float for ASCII, unsigned char for BINARY
1349                append out "COLOR_SCALARS $label $elemSize\n"
1350            } elseif { $_comp2type($cname) == "normals" } {
1351                # elemSize must equal 3
1352                append out "NORMALS $label double\n"
1353            } elseif { $_comp2type($cname) == "scalars" } {
1354                # elemSize can be 1, 2, 3 or 4
1355                append out "SCALARS $label double $elemSize\n"
1356                append out "LOOKUP_TABLE default\n"
1357            } elseif { $_comp2type($cname) == "tcoords" } {
1358                # elemSize must be 1, 2, or 3
1359                append out "TEXTURE_COORDINATES $label $elemSize double\n"
1360            } elseif { $_comp2type($cname) == "tensors" } {
1361                # elemSize must equal 9
1362                append out "TENSORS $label double\n"
1363            } elseif { $_comp2type($cname) == "vectors" } {
1364                # elemSize must equal 3
1365                append out "VECTORS $label double\n"
1366            } else {
1367                error "unknown element type \"$_comp2type($cname)\""
1368            }
1369        }
1370        append out [$vector range 0 end]
1371        append out "\n"
1372        if 0 {
1373            VerifyVtkDataSet $out
1374        }
1375        return $out
1376    }
1377    error "can't find vtkdata for $cname"
1378}
1379
1380#
1381# BuildPointsOnMesh --
1382#
1383#       Parses the field XML description to build a mesh and values vector
1384#       representing the field.  Right now we handle the deprecated types
1385#       of "cloud", "unirect2d", and "unirect3d" (mostly for flows).
1386#
1387itcl::body Rappture::Field::BuildPointsOnMesh {cname} {
1388    #
1389    # More complex 2D/3D data is represented by a mesh
1390    # object and an associated vector for field values.
1391    #
1392    set path [$_field get $cname.mesh]
1393    if {[$_xmlobj element $path] == ""} {
1394        # Unknown mesh designated.
1395        return 0
1396    }
1397    set viewer [$_field get "about.view"]
1398    if { $viewer != "" } {
1399        set _viewer $viewer
1400    }
1401    set element [$_xmlobj element -as type $path]
1402    set name $cname
1403    regsub -all { } $name {_} name
1404    set _fld2Label($name) $name
1405    set label [hints zlabel]
1406    if { $label != "" } {
1407        set _fld2Label($name) $label
1408    }
1409    set _fld2Units($name) [hints zunits]
1410    set _fld2Components($name) $_comp2size($cname)
1411    lappend _comp2fldName($cname) $name
1412
1413    # Handle bizarre cases that hopefully will be deprecated.
1414    if { $element == "unirect3d" } {
1415        # Special case: unirect3d (should be deprecated) + flow.
1416        if { [$_field element $cname.extents] != "" } {
1417            set vectorsize [$_field get $cname.extents]
1418        } else {
1419            set vectorsize 1
1420        }
1421        set _type unirect3d
1422        set _dim 3
1423        if { $_viewer == "" } {
1424            set _viewer flowvis
1425        }
1426        set _comp2dims($cname) "3D"
1427        set _comp2unirect3d($cname) \
1428            [Rappture::Unirect3d \#auto $_xmlobj $_field $cname $vectorsize]
1429        set _comp2style($cname) [$_field get $cname.style]
1430        set limits {}
1431        foreach axis { x y z } {
1432            lappend limits $axis [$_comp2unirect3d($cname) limits $axis]
1433        }
1434        # Get the data limits
1435        set vector [$_comp2unirect3d($cname) values]
1436        set minmax [VectorLimits $vector $vectorsize]
1437        lappend limits $cname $minmax
1438        lappend limits v      $minmax
1439        set _comp2limits($cname) $limits
1440        if {[$_field element $cname.flow] != ""} {
1441            set _comp2flowhints($cname) \
1442                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1443        }
1444        incr _counter
1445        return 1
1446    }
1447    if { $element == "unirect2d" && [$_field element $cname.flow] != "" } {
1448        # Special case: unirect2d (normally deprecated) + flow.
1449        if { [$_field element $cname.extents] != "" } {
1450            set vectorsize [$_field get $cname.extents]
1451        } else {
1452            set vectorsize 1
1453        }
1454        set _type unirect2d
1455        set _dim 2
1456        if { $_viewer == "" } {
1457            set _viewer "flowvis"
1458        }
1459        set _comp2dims($cname) "2D"
1460        set _comp2unirect2d($cname) \
1461            [Rappture::Unirect2d \#auto $_xmlobj $path]
1462        set _comp2style($cname) [$_field get $cname.style]
1463        set _comp2flowhints($cname) \
1464            [Rappture::FlowHints ::\#auto $_field $cname $_units]
1465        set _values [$_field get $cname.values]
1466        set limits {}
1467        foreach axis { x y z } {
1468            lappend limits $axis [$_comp2unirect2d($cname) limits $axis]
1469        }
1470        set xv [blt::vector create \#auto]
1471        $xv set $_values
1472        set minmax [VectorLimits $xv $vectorsize]
1473        lappend limits $cname $minmax
1474        lappend limits v $minmax
1475        blt::vector destroy $xv
1476        set _comp2limits($cname) $limits
1477        incr _counter
1478        return 1
1479    }
1480    switch -- $element {
1481        "cloud" {
1482            set mesh [Rappture::Cloud::fetch $_xmlobj $path]
1483            set _type cloud
1484        }
1485        "mesh" {
1486            set mesh [Rappture::Mesh::fetch $_xmlobj $path]
1487            set _type mesh
1488        }           
1489        "unirect2d" {
1490            if { $_viewer == "" } {
1491                set _viewer "heightmap"
1492            }
1493            set mesh [Rappture::Unirect2d::fetch $_xmlobj $path]
1494            set _type unirect2d
1495        }
1496    }
1497    if { ![$mesh isvalid] } {
1498        puts stderr "Mesh is invalid"
1499        return 0
1500    }
1501    set _dim [$mesh dimensions]
1502    if { $_dim == 3 } {
1503        set dim 0
1504        foreach axis {x y z} {
1505            foreach {min max} [$mesh limits $axis] {
1506                if { $min < $max } {
1507                    incr dim
1508                }
1509            }
1510        }
1511        if { $dim != 3 } {
1512            set _dim $dim
1513        }
1514    }
1515
1516    if {$_dim < 2} {
1517        puts stderr "ERROR: Can't convert 1D cloud/mesh to curve.  Please use curve output for 1D meshes."
1518        return 0
1519
1520        # 1D data: Create vectors for graph widget.
1521        # The prophet tool currently outputs 1D clouds with fields
1522        # Band Structure Lab used to (see isosurface1 test in rappture-bat)
1523        #
1524        # Is there a natural growth path in generating output from 1D to
1525        # higher dimensions?  If there isn't, let's kill this in favor
1526        # or explicitly using a <curve> instead.  Otherwise, the features
1527        # (methods such as xmarkers) or the <curve> need to be added
1528        # to the <field>.
1529        #
1530        #set xv [blt::vector create x$_counter]
1531        #set yv [blt::vector create y$_counter]
1532
1533        # This only works with a Cloud mesh type, since the points method
1534        # is not implemented for the Mesh object
1535        #$xv set [$mesh points]
1536        # TODO: Put field values in yv
1537        #set _comp2dims($cname) "1D"
1538        #set _comp2xy($cname) [list $xv $yv]
1539        #incr _counter
1540        #return 1
1541    }
1542    if {$_dim == 2} {
1543        # 2D data: By default surface or contour plot using heightmap widget.
1544        set v [blt::vector create \#auto]
1545        $v set [$_field get $cname.values]
1546        if { [$v length] == 0 } {
1547            return 0
1548        }
1549        if { $_viewer == "" } {
1550            set _viewer "contour"
1551        }
1552        set numFieldValues [$v length]
1553        set numComponentsPerTuple [numComponents $cname]
1554        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1555            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1556            return 0
1557        }
1558        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1559        if { $_comp2assoc($cname) == "pointdata" } {
1560            set numPoints [$mesh numpoints]
1561            if { $numPoints != $numFieldTuples } {
1562                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1563                return 0
1564            }
1565        } elseif { $_comp2assoc($cname) == "celldata" } {
1566            set numCells [$mesh numcells]
1567            if { $numCells != $numFieldTuples } {
1568                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1569                return 0
1570            }
1571        }
1572        set _comp2dims($cname) "[$mesh dimensions]D"
1573        set _comp2mesh($cname) [list $mesh $v]
1574        set _comp2style($cname) [$_field get $cname.style]
1575        if {[$_field element $cname.flow] != ""} {
1576            set _comp2flowhints($cname) \
1577                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1578        }
1579        incr _counter
1580        array unset _comp2limits $cname
1581        foreach axis { x y z } {
1582            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1583        }
1584        set minmax [VectorLimits $v $_comp2size($cname)]
1585        lappend _comp2limits($cname) $cname $minmax
1586        lappend _comp2limits($cname) v $minmax
1587        return 1
1588    }
1589    if {$_dim == 3} {
1590        # 3D data: By default isosurfaces plot using isosurface widget.
1591        if { $_viewer == "" } {
1592            set _viewer "isosurface"
1593        }
1594        set v [blt::vector create \#auto]
1595        $v set [$_field get $cname.values]
1596        if { [$v length] == 0 } {
1597            return 0
1598        }
1599        set numFieldValues [$v length]
1600        set numComponentsPerTuple [numComponents $cname]
1601        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1602            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1603            return 0
1604        }
1605        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1606        if { $_comp2assoc($cname) == "pointdata" } {
1607            set numPoints [$mesh numpoints]
1608            if { $numPoints != $numFieldTuples } {
1609                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1610                return 0
1611            }
1612        } elseif { $_comp2assoc($cname) == "celldata" } {
1613            set numCells [$mesh numcells]
1614            if { $numCells != $numFieldTuples } {
1615                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1616                return 0
1617            }
1618        }
1619        set _comp2dims($cname) "[$mesh dimensions]D"
1620        set _comp2mesh($cname) [list $mesh $v]
1621        set _comp2style($cname) [$_field get $cname.style]
1622        if {[$_field element $cname.flow] != ""} {
1623            set _comp2flowhints($cname) \
1624                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1625        }
1626        incr _counter
1627        foreach axis { x y z } {
1628            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1629        }
1630        set minmax [VectorLimits $v $_comp2size($cname)]
1631        lappend _comp2limits($cname) $cname $minmax
1632        lappend _comp2limits($cname) v $minmax
1633        return 1
1634    }
1635    error "unhandled case in field dim=$_dim element=$element"
1636}
1637
1638itcl::body Rappture::Field::AvsToVtk { cname contents } {
1639    package require vtk
1640
1641    set reader $this-datasetreader
1642    vtkAVSucdReader $reader
1643
1644    # Write the contents to a file just in case it's binary.
1645    set tmpfile $cname[pid].ucd
1646    set f [open "$tmpfile" "w"]
1647    fconfigure $f -translation binary -encoding binary
1648    puts $f $contents
1649    close $f
1650    $reader SetFileName $tmpfile
1651    $reader Update
1652    file delete $tmpfile
1653
1654    set tmpfile $cname[pid].vtk
1655    set writer $this-datasetwriter
1656    vtkDataSetWriter $writer
1657    $writer SetInputConnection [$reader GetOutputPort]
1658    $writer SetFileName $tmpfile
1659    $writer Write
1660    rename $reader ""
1661    rename $writer ""
1662
1663    set f [open "$tmpfile" "r"]
1664    fconfigure $f -translation binary -encoding binary
1665    set vtkdata [read $f]
1666    close $f
1667    file delete $tmpfile
1668    return $vtkdata
1669}
1670
1671itcl::body Rappture::Field::DicomToVtk { cname path } {
1672    package require vtk
1673
1674    if { ![file exists $path] } {
1675        puts stderr "path \"$path\" doesn't exist."
1676        return 0
1677    }
1678
1679    if { [file isdir $path] } {
1680        set files [glob -nocomplain $path/*.dcm]
1681        if { [llength $files] == 0 } {
1682            set files [glob -nocomplain $path/*]
1683            if { [llength $files] == 0 } {
1684                puts stderr "no dicom files found in \"$path\""
1685                return 0
1686            }
1687        }
1688
1689        #array set data [Rappture::DicomToVtk files $files]
1690        array set data [Rappture::DicomToVtk dir $path]
1691    } else {
1692        array set data [Rappture::DicomToVtk files [list $path]]
1693    }
1694
1695    if 0 {
1696        foreach key [array names data] {
1697            if {$key == "vtkdata"} {
1698                if 0 {
1699                    set f [open /tmp/$cname.vtk "w"]
1700                    fconfigure $f -translation binary -encoding binary
1701                    puts -nonewline $f $data(vtkdata)
1702                    close $f
1703                }
1704            } else {
1705                puts stderr "$key = \"$data($key)\""
1706            }
1707        }
1708    }
1709
1710    # Save viewer choice
1711    set viewer $_viewer
1712    ReadVtkDataSet $cname $data(vtkdata)
1713    # Restore viewer choice (ReadVtkDataSet wants to set it to contour/isosurface)
1714    set _viewer $viewer
1715    return $data(vtkdata)
1716}
1717
1718itcl::body Rappture::Field::GetTypeAndSize { cname } {
1719    array set type2components {
1720        "colorscalars"         4
1721        "normals"              3
1722        "scalars"              1
1723        "tcoords"              2
1724        "tensors"              9
1725        "vectors"              3
1726    }
1727    set type [$_field get $cname.elemtype]
1728    if { $type == "" } {
1729        set type "scalars"
1730    }
1731    if { ![info exists type2components($type)] } {
1732        error "unknown <elemtype> \"$type\" in field"
1733    }
1734    set size [$_field get $cname.elemsize]
1735    if { $size == "" } {
1736        set size $type2components($type)
1737    }
1738    set _comp2type($cname) $type
1739    set _comp2size($cname) $size
1740}
1741
1742itcl::body Rappture::Field::GetAssociation { cname } {
1743    set assoc [$_field get $cname.association]
1744    if { $assoc == "" } {
1745        set _comp2assoc($cname) "pointdata"
1746        return
1747    }
1748    switch -- $assoc {
1749        "pointdata" - "celldata" - "fielddata" {
1750            set _comp2assoc($cname) $assoc
1751            return
1752        }
1753        default {
1754            error "unknown field association \"$assoc\""
1755        }
1756    }
1757}
1758
1759#
1760# Compute the per-component limits or limits of vector magnitudes
1761#
1762itcl::body Rappture::Field::VectorLimits {vector vectorsize {comp -1}} {
1763    if {$vectorsize == 1} {
1764        set minmax [$vector limits]
1765    } else {
1766        set len [$vector length]
1767        if {[expr $len % $vectorsize] != 0} {
1768            error "Invalid vectorsize: $vectorsize"
1769        }
1770        if {$comp > $vectorsize-1} {
1771            error "Invalid vector component: $comp"
1772        }
1773        set numTuples [expr ($len/$vectorsize)]
1774        for {set i 0} {$i < $numTuples} {incr i} {
1775            if {$comp >= 0} {
1776                set idx [expr ($i * $vectorsize + $comp)]
1777                set val [$vector index $idx]
1778            } else {
1779                set idx [expr ($i * $vectorsize)]
1780                set mag 0
1781                for {set j 0} {$j < $vectorsize} {incr j} {
1782                    set val [$vector index $idx]
1783                    set mag [expr ($mag + $val * $val)]
1784                    incr idx
1785                }
1786                set val [expr (sqrt($mag))]
1787            }
1788            if (![info exists minmax]) {
1789                set minmax [list $val $val]
1790            } else {
1791                if {$val < [lindex $minmax 0]} {
1792                    lset minmax 0 $val
1793                }
1794                if {$val > [lindex $minmax 1]} {
1795                    lset minmax 1 $val
1796                }
1797            }
1798        }
1799    }
1800    return $minmax
1801}
Note: See TracBrowser for help on using the repository browser.