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

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

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