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

Last change on this file since 5178 was 5178, checked in by ldelgass, 5 years ago

merge r5159 from trunk

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