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

Last change on this file since 5312 was 5174, checked in by ldelgass, 10 years ago

merge r5158:r5160 from trunk
File size: 60.8 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                unset contents
918                ReadVtkDataSet $cname $vtkdata
919                if 0 {
920                    set f [open /tmp/$_path.$cname.vtk "w"]
921                    puts -nonewline $f $vtkdata
922                    close $f
923                }
924            } else {
925                unset contents
926                puts stderr "Can't parse dx data"
927            }
928            if { $_alwaysConvertDX ||
929                 ($_viewer != "nanovis" && $_viewer != "flowvis") } {
930                set _type "vtk"
931                set _comp2vtk($cname) $vtkdata
932            } else {
933                set _type "dx"
934                set _comp2dx($cname) $data
935            }
936            unset data
937            unset vtkdata
938            set _comp2style($cname) [$_field get $cname.style]
939            if {[$_field element $cname.flow] != ""} {
940                set _comp2flowhints($cname) \
941                    [Rappture::FlowHints ::\#auto $_field $cname $_units]
942            }
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 "ERROR: All components of field \"$_path\" are invalid."
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: A 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    # 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        puts stderr "ERROR: Unknown mesh \"$path\""
1396        return 0
1397    }
1398    set viewer [$_field get "about.view"]
1399    if { $viewer != "" } {
1400        set _viewer $viewer
1401    }
1402    set element [$_xmlobj element -as type $path]
1403    set name $cname
1404    regsub -all { } $name {_} name
1405    set _fld2Label($name) $name
1406    set label [hints zlabel]
1407    if { $label != "" } {
1408        set _fld2Label($name) $label
1409    }
1410    set _fld2Units($name) [hints zunits]
1411    set _fld2Components($name) $_comp2size($cname)
1412    lappend _comp2fldName($cname) $name
1413
1414    # Handle bizarre cases that hopefully will be deprecated.
1415    if { $element == "unirect3d" } {
1416        # Special case: unirect3d (should be deprecated) + flow.
1417        if { [$_field element $cname.extents] != "" } {
1418            set vectorsize [$_field get $cname.extents]
1419        } else {
1420            set vectorsize 1
1421        }
1422        set _type unirect3d
1423        set _dim 3
1424        if { $_viewer == "" } {
1425            set _viewer flowvis
1426        }
1427        set _comp2dims($cname) "3D"
1428        set _comp2unirect3d($cname) \
1429            [Rappture::Unirect3d \#auto $_xmlobj $_field $cname $vectorsize]
1430        set _comp2style($cname) [$_field get $cname.style]
1431        set limits {}
1432        foreach axis { x y z } {
1433            lappend limits $axis [$_comp2unirect3d($cname) limits $axis]
1434        }
1435        # Get the data limits
1436        set vector [$_comp2unirect3d($cname) values]
1437        set minmax [VectorLimits $vector $vectorsize]
1438        lappend limits $cname $minmax
1439        lappend limits v      $minmax
1440        set _comp2limits($cname) $limits
1441        if {[$_field element $cname.flow] != ""} {
1442            set _comp2flowhints($cname) \
1443                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1444        }
1445        incr _counter
1446        return 1
1447    }
1448    if { $element == "unirect2d" && [$_field element $cname.flow] != "" } {
1449        # Special case: unirect2d (normally deprecated) + flow.
1450        if { [$_field element $cname.extents] != "" } {
1451            set vectorsize [$_field get $cname.extents]
1452        } else {
1453            set vectorsize 1
1454        }
1455        set _type unirect2d
1456        set _dim 2
1457        if { $_viewer == "" } {
1458            set _viewer "flowvis"
1459        }
1460        set _comp2dims($cname) "2D"
1461        set _comp2unirect2d($cname) \
1462            [Rappture::Unirect2d \#auto $_xmlobj $path]
1463        set _comp2style($cname) [$_field get $cname.style]
1464        set _comp2flowhints($cname) \
1465            [Rappture::FlowHints ::\#auto $_field $cname $_units]
1466        set _values [$_field get $cname.values]
1467        set limits {}
1468        foreach axis { x y z } {
1469            lappend limits $axis [$_comp2unirect2d($cname) limits $axis]
1470        }
1471        set xv [blt::vector create \#auto]
1472        $xv set $_values
1473        set minmax [VectorLimits $xv $vectorsize]
1474        lappend limits $cname $minmax
1475        lappend limits v $minmax
1476        blt::vector destroy $xv
1477        set _comp2limits($cname) $limits
1478        incr _counter
1479        return 1
1480    }
1481    switch -- $element {
1482        "cloud" {
1483            set mesh [Rappture::Cloud::fetch $_xmlobj $path]
1484            set _type cloud
1485        }
1486        "mesh" {
1487            set mesh [Rappture::Mesh::fetch $_xmlobj $path]
1488            set _type mesh
1489        }
1490        "unirect2d" {
1491            if { $_viewer == "" } {
1492                set _viewer "heightmap"
1493            }
1494            set mesh [Rappture::Unirect2d::fetch $_xmlobj $path]
1495            set _type unirect2d
1496        }
1497    }
1498    if { ![$mesh isvalid] } {
1499        puts stderr "Mesh is invalid"
1500        return 0
1501    }
1502    set _dim [$mesh dimensions]
1503    if { $_dim == 3 } {
1504        set dim 0
1505        foreach axis {x y z} {
1506            foreach {min max} [$mesh limits $axis] {
1507                if { $min < $max } {
1508                    incr dim
1509                }
1510            }
1511        }
1512        if { $dim != 3 } {
1513            set _dim $dim
1514        }
1515    }
1516
1517    if {$_dim < 2} {
1518        puts stderr "ERROR: Can't convert 1D cloud/mesh to curve.  Please use curve output for 1D meshes."
1519        return 0
1520
1521        # 1D data: Create vectors for graph widget.
1522        # The prophet tool currently outputs 1D clouds with fields
1523        # Band Structure Lab used to (see isosurface1 test in rappture-bat)
1524        #
1525        # Is there a natural growth path in generating output from 1D to
1526        # higher dimensions?  If there isn't, let's kill this in favor
1527        # or explicitly using a <curve> instead.  Otherwise, the features
1528        # (methods such as xmarkers) of the <curve> need to be added
1529        # to the <field>.
1530        #
1531        #set xv [blt::vector create x$_counter]
1532        #set yv [blt::vector create y$_counter]
1533
1534        # This only works with a Cloud mesh type, since the points method
1535        # is not implemented for the Mesh object
1536        #$xv set [$mesh points]
1537        # TODO: Put field values in yv
1538        #set _comp2dims($cname) "1D"
1539        #set _comp2xy($cname) [list $xv $yv]
1540        #incr _counter
1541        #return 1
1542    }
1543    if {$_dim == 2} {
1544        # 2D data: By default surface or contour plot using heightmap widget.
1545        set v [blt::vector create \#auto]
1546        $v set [$_field get $cname.values]
1547        if { [$v length] == 0 } {
1548            return 0
1549        }
1550        if { $_viewer == "" } {
1551            set _viewer "contour"
1552        }
1553        set numFieldValues [$v length]
1554        set numComponentsPerTuple [numComponents $cname]
1555        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1556            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1557            return 0
1558        }
1559        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1560        if { $_comp2assoc($cname) == "pointdata" } {
1561            set numPoints [$mesh numpoints]
1562            if { $numPoints != $numFieldTuples } {
1563                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1564                return 0
1565            }
1566        } elseif { $_comp2assoc($cname) == "celldata" } {
1567            set numCells [$mesh numcells]
1568            if { $numCells != $numFieldTuples } {
1569                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1570                return 0
1571            }
1572        }
1573        set _comp2dims($cname) "[$mesh dimensions]D"
1574        set _comp2mesh($cname) [list $mesh $v]
1575        set _comp2style($cname) [$_field get $cname.style]
1576        if {[$_field element $cname.flow] != ""} {
1577            set _comp2flowhints($cname) \
1578                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1579        }
1580        incr _counter
1581        array unset _comp2limits $cname
1582        foreach axis { x y z } {
1583            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1584        }
1585        set minmax [VectorLimits $v $_comp2size($cname)]
1586        lappend _comp2limits($cname) $cname $minmax
1587        lappend _comp2limits($cname) v $minmax
1588        return 1
1589    }
1590    if {$_dim == 3} {
1591        # 3D data: By default isosurfaces plot using isosurface widget.
1592        if { $_viewer == "" } {
1593            set _viewer "isosurface"
1594        }
1595        set v [blt::vector create \#auto]
1596        $v set [$_field get $cname.values]
1597        if { [$v length] == 0 } {
1598            puts stderr "ERROR: No field values"
1599            return 0
1600        }
1601        set numFieldValues [$v length]
1602        set numComponentsPerTuple [numComponents $cname]
1603        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1604            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1605            return 0
1606        }
1607        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1608        if { $_comp2assoc($cname) == "pointdata" } {
1609            set numPoints [$mesh numpoints]
1610            if { $numPoints != $numFieldTuples } {
1611                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1612                return 0
1613            }
1614        } elseif { $_comp2assoc($cname) == "celldata" } {
1615            set numCells [$mesh numcells]
1616            if { $numCells != $numFieldTuples } {
1617                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1618                return 0
1619            }
1620        }
1621        set _comp2dims($cname) "[$mesh dimensions]D"
1622        set _comp2mesh($cname) [list $mesh $v]
1623        set _comp2style($cname) [$_field get $cname.style]
1624        if {[$_field element $cname.flow] != ""} {
1625            set _comp2flowhints($cname) \
1626                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1627        }
1628        incr _counter
1629        foreach axis { x y z } {
1630            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1631        }
1632        set minmax [VectorLimits $v $_comp2size($cname)]
1633        lappend _comp2limits($cname) $cname $minmax
1634        lappend _comp2limits($cname) v $minmax
1635        return 1
1636    }
1637    error "unhandled case in field dim=$_dim element=$element"
1638}
1639
1640itcl::body Rappture::Field::AvsToVtk { cname contents } {
1641    package require vtk
1642
1643    set reader $this-datasetreader
1644    vtkAVSucdReader $reader
1645
1646    # Write the contents to a file just in case it's binary.
1647    set tmpfile $cname[pid].ucd
1648    set f [open "$tmpfile" "w"]
1649    fconfigure $f -translation binary -encoding binary
1650    puts $f $contents
1651    close $f
1652    $reader SetFileName $tmpfile
1653    $reader Update
1654    file delete $tmpfile
1655
1656    set tmpfile $cname[pid].vtk
1657    set writer $this-datasetwriter
1658    vtkDataSetWriter $writer
1659    $writer SetInputConnection [$reader GetOutputPort]
1660    $writer SetFileName $tmpfile
1661    $writer SetFileTypeToBinary
1662    $writer Write
1663    rename $reader ""
1664    rename $writer ""
1665
1666    set f [open "$tmpfile" "r"]
1667    fconfigure $f -translation binary -encoding binary
1668    set vtkdata [read $f]
1669    close $f
1670    file delete $tmpfile
1671    return $vtkdata
1672}
1673
1674itcl::body Rappture::Field::DicomToVtk { cname path } {
1675    package require vtk
1676
1677    if { ![file exists $path] } {
1678        puts stderr "path \"$path\" doesn't exist."
1679        return 0
1680    }
1681
1682    if { [file isdir $path] } {
1683        set files [glob -nocomplain $path/*.dcm]
1684        if { [llength $files] == 0 } {
1685            set files [glob -nocomplain $path/*]
1686            if { [llength $files] == 0 } {
1687                puts stderr "no dicom files found in \"$path\""
1688                return 0
1689            }
1690        }
1691
1692        #array set data [Rappture::DicomToVtk files $files]
1693        array set data [Rappture::DicomToVtk dir $path]
1694    } else {
1695        array set data [Rappture::DicomToVtk files [list $path]]
1696    }
1697
1698    if 0 {
1699        foreach key [array names data] {
1700            if {$key == "vtkdata"} {
1701                if 0 {
1702                    set f [open /tmp/$cname.vtk "w"]
1703                    fconfigure $f -translation binary -encoding binary
1704                    puts -nonewline $f $data(vtkdata)
1705                    close $f
1706                }
1707            } else {
1708                puts stderr "$key = \"$data($key)\""
1709            }
1710        }
1711    }
1712
1713    # Save viewer choice
1714    set viewer $_viewer
1715    ReadVtkDataSet $cname $data(vtkdata)
1716    # Restore viewer choice (ReadVtkDataSet wants to set it to contour/isosurface)
1717    set _viewer $viewer
1718    return $data(vtkdata)
1719}
1720
1721itcl::body Rappture::Field::GetTypeAndSize { cname } {
1722    array set type2components {
1723        "colorscalars"         4
1724        "normals"              3
1725        "scalars"              1
1726        "tcoords"              2
1727        "tensors"              9
1728        "vectors"              3
1729    }
1730    set type [$_field get $cname.elemtype]
1731    if { $type == "" } {
1732        set type "scalars"
1733    }
1734    if { ![info exists type2components($type)] } {
1735        error "unknown <elemtype> \"$type\" in field"
1736    }
1737    set size [$_field get $cname.elemsize]
1738    if { $size == "" } {
1739        set size $type2components($type)
1740    }
1741    set _comp2type($cname) $type
1742    set _comp2size($cname) $size
1743}
1744
1745itcl::body Rappture::Field::GetAssociation { cname } {
1746    set assoc [$_field get $cname.association]
1747    if { $assoc == "" } {
1748        set _comp2assoc($cname) "pointdata"
1749        return
1750    }
1751    switch -- $assoc {
1752        "pointdata" - "celldata" - "fielddata" {
1753            set _comp2assoc($cname) $assoc
1754            return
1755        }
1756        default {
1757            error "unknown field association \"$assoc\""
1758        }
1759    }
1760}
1761
1762#
1763# Compute the per-component limits or limits of vector magnitudes
1764#
1765itcl::body Rappture::Field::VectorLimits {vector vectorsize {comp -1}} {
1766    if {$vectorsize == 1} {
1767        set minmax [$vector limits]
1768    } else {
1769        set len [$vector length]
1770        if {[expr $len % $vectorsize] != 0} {
1771            error "Invalid vectorsize: $vectorsize"
1772        }
1773        if {$comp > $vectorsize-1} {
1774            error "Invalid vector component: $comp"
1775        }
1776        set numTuples [expr ($len/$vectorsize)]
1777        for {set i 0} {$i < $numTuples} {incr i} {
1778            if {$comp >= 0} {
1779                set idx [expr ($i * $vectorsize + $comp)]
1780                set val [$vector index $idx]
1781            } else {
1782                set idx [expr ($i * $vectorsize)]
1783                set mag 0
1784                for {set j 0} {$j < $vectorsize} {incr j} {
1785                    set val [$vector index $idx]
1786                    set mag [expr ($mag + $val * $val)]
1787                    incr idx
1788                }
1789                set val [expr (sqrt($mag))]
1790            }
1791            if (![info exists minmax]) {
1792                set minmax [list $val $val]
1793            } else {
1794                if {$val < [lindex $minmax 0]} {
1795                    lset minmax 0 $val
1796                }
1797                if {$val > [lindex $minmax 1]} {
1798                    lset minmax 1 $val
1799                }
1800            }
1801        }
1802    }
1803    return $minmax
1804}
Note: See TracBrowser for help on using the repository browser.