source: branches/1.3/gui/scripts/field.tcl @ 4665

Last change on this file since 4665 was 4665, checked in by gah, 10 years ago

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