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

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

comment fixes
File size: 59.8 KB
Line 
1# -*- mode: tcl; indent-tabs-mode: nil -*-
2
3# ----------------------------------------------------------------------
4#  COMPONENT: field - extracts data from an XML description of a field
5#
6#  This object represents one field in an XML description of a device.
7#  It simplifies the process of extracting data vectors that represent
8#  the field.
9# ======================================================================
10#  AUTHOR:  Michael McLennan, Purdue University
11#  Copyright (c) 2004-2012  HUBzero Foundation, LLC
12#
13#  See the file "license.terms" for information on usage and
14#  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
15# ======================================================================
16
17# TODO:
18#    Vector field limits are wrong: need to compute magnitude limits and
19#    component-wise limits.
20
21#
22# Possible field dataset types:
23#
24# 2D Datasets
25#       vtk             (range of z-axis is zero).
26#       unirect2d       (deprecated)
27#       cloud           (x,y point coordinates) (deprecated)
28#       mesh
29# 3D Datasets
30#       vtk
31#       unirect3d       (deprecated)
32#       cloud           (x,y,z coordinates) (deprecated)
33#       mesh
34#       dx              (FIXME: make dx-to-vtk converter work)
35#       ucd avs
36#
37# Viewers:
38#       Format     Dim  Description                     Viewer          Server
39#       vtk         2   vtk file data.                  contour         vtkvis
40#       vtk         3   vtk file data.                  isosurface      vtkvis
41#       mesh        2   points-on-mesh                  heightmap       vtkvis
42#       mesh        3   points-on-mesh                  isosurface      vtkvis
43#       dx          3   DX                              volume          nanovis
44#       unirect2d   2   unirect3d + extents > 1 flow    flow            nanovis
45#       unirect3d   3   unirect2d + extents > 1 flow    flow            nanovis
46#       
47# With <views>, can specify which viewer for specific datasets.  So it's OK
48# for the same dataset to be viewed in more than one way.
49#  o Any 2D dataset can be viewed as a contour/heightmap.
50#  o Any 3D dataset can be viewed as a isosurface.
51#  o Any 2D dataset with vector data can be streamlines or flow. 
52#  o Any 3D uniform rectilinear dataset can be viewed as a volume.
53#  o Any 3D dataset with vector data can be streamlines or flow.
54#
55# Need <views> to properly do things like qdot: volume with a polydata
56# transparent shell.  The view will combine the two objects <field> and
57# <drawing> (??) into a single viewer.
58#
59package require Itcl
60package require BLT
61
62namespace eval Rappture {
63    # forward declaration
64}
65
66itcl::class Rappture::Field {
67    protected variable _dim 0;          # Dimension of the mesh
68    private variable _xmlobj "";        # ref to XML obj with field data
69    private variable _limits;           # maps axis name => {z0 z1} limits
70    private variable _field ""
71    private variable _comp2fldName ;    # cname => field names.
72    private variable _comp2type ;       # cname => type (e.g. "vectors")
73    private variable _comp2size ;       # cname => # of components in element
74    private variable _comp2assoc;       # cname => association (e.g. pointdata)
75    private variable _fld2Components;   # field name => number of components
76    private variable _fld2Label;        # field name => label
77    private variable _fld2Units;        # field name => units
78    private variable _hints
79    private variable _viewer "";        # Hints which viewer to use
80    private variable _xv "";            # For 1D meshes only.  Holds the points
81    private variable _isValid 0;        # Indicates if the field contains
82                                        # valid data.
83    private variable _isValidComponent; #  Array of valid components found
84    private variable _alwaysConvertDX 0;
85
86    constructor {xmlobj path} {
87        # defined below
88    }
89    destructor {
90        # defined below
91    }
92    public method blob { cname }
93    public method components {args}
94    public method controls {option args}
95    public method extents {{cname -overall}}
96    public method numComponents {cname}
97    public method fieldlimits {}
98    public method valueLimits { cname }
99    public method flowhints { cname }
100    public method hints {{key ""}}
101    public method isunirect2d {}
102    public method isunirect3d {}
103    public method limits {axis}
104    public method mesh {{cname -overall}}
105    public method style { cname }
106    public method type {}
107    public method values { cname }
108    public method vtkdata {cname}
109    public method xErrorValues { cname } {
110    }
111    public method yErrorValues { cname } {
112    }
113
114    public method fieldnames { cname } {
115        if { ![info exists _comp2fldName($cname)] } {
116            return ""
117        }
118        return $_comp2fldName($cname)
119    }
120    public method fieldinfo { fname } {
121        lappend out $_fld2Label($fname)
122        lappend out $_fld2Units($fname)
123        lappend out $_fld2Components($fname)
124        return $out
125    }
126    public method isvalid {} {
127        return $_isValid
128    }
129    public method viewer {} {
130        return $_viewer
131    }
132    protected method Build {}
133    protected method _getValue {expr}
134
135    private variable _path "";          # Path of this object in the XML
136    private variable _units ""   ;      # system of units for this field
137    private variable _zmax 0     ;# length of the device
138
139    private variable _comp2dims  ;# maps component name => dimensionality
140    private variable _comp2xy    ;# maps component name => x,y vectors
141    private variable _comp2vtk   ;# maps component name => vtk file data
142    private variable _comp2dx    ;# maps component name => OpenDX data
143    private variable _comp2unirect2d ;# maps component name => unirect2d obj
144    private variable _comp2unirect3d ;# maps component name => unirect3d obj
145    private variable _comp2style ;# maps component name => style settings
146    private variable _comp2cntls ;# maps component name => x,y control points
147    private variable _comp2extents
148    private variable _comp2limits;      #  Array of limits per component
149    private variable _type ""
150    private variable _comp2flowhints
151    private variable _comp2mesh
152    private common _counter 0    ;# counter for unique vector names
153
154    private method AvsToVtk { cname contents }
155    private method DicomToVtk { cname contents }
156    private method BuildPointsOnMesh { cname }
157    protected method GetAssociation { cname }
158    protected method GetTypeAndSize { cname }
159    protected method ReadVtkDataSet { cname contents }
160    private method InitHints {}
161
162    private method VerifyVtkDataSet { contents }
163    private method VectorLimits { vector vectorsize {comp -1} }
164    private variable _values ""
165}
166
167# ----------------------------------------------------------------------
168# CONSTRUCTOR
169# ----------------------------------------------------------------------
170itcl::body Rappture::Field::constructor {xmlobj path} {
171    package require vtk
172    if {![Rappture::library isvalid $xmlobj]} {
173        error "bad value \"$xmlobj\": should be Rappture::library"
174    }
175    set _xmlobj $xmlobj
176    set _path $path
177    set _field [$xmlobj element -as object $path]
178    set _units [$_field get units]
179
180    set xunits [$xmlobj get units]
181    if {"" == $xunits || "arbitrary" == $xunits} {
182        set xunits "um"
183    }
184
185    # determine the overall size of the device
186    set z0 [set z1 0]
187    foreach elem [$_xmlobj children components] {
188        switch -glob -- $elem {
189            box* {
190                if {![regexp {[0-9]$} $elem]} {
191                    set elem "${elem}0"
192                }
193                set z0 [$_xmlobj get components.$elem.corner0]
194                set z0 [Rappture::Units::convert $z0 \
195                    -context $xunits -to $xunits -units off]
196
197                set z1 [$_xmlobj get components.$elem.corner1]
198                set z1 [Rappture::Units::convert $z1 \
199                    -context $xunits -to $xunits -units off]
200
201                set _limits($elem) [list $z0 $z1]
202            }
203        }
204    }
205    set _zmax $z1
206
207    # build up vectors for various components of the field
208    Build
209    InitHints
210}
211
212# ----------------------------------------------------------------------
213# DESTRUCTOR
214# ----------------------------------------------------------------------
215itcl::body Rappture::Field::destructor {} {
216    itcl::delete object $_field
217    # don't destroy the _xmlobj! we don't own it!
218
219    foreach name [array names _comp2xy] {
220        eval blt::vector destroy $_comp2xy($name)
221    }
222    foreach name [array names _comp2unirect2d] {
223        itcl::delete object $_comp2unirect2d($name)
224    }
225    foreach name [array names _comp2unirect3d] {
226        itcl::delete object $_comp2unirect3d($name)
227    }
228    foreach name [array names _comp2flowhints] {
229        itcl::delete object $_comp2flowhints($name)
230    }
231    foreach name [array names _comp2mesh] {
232        # Data is in the form of a mesh and a vector.
233        foreach { mesh vector } $_comp2mesh($name) break
234        # Release the mesh (may be shared)
235        set class [$mesh info class]
236        ${class}::release $mesh
237        # Destroy the vector
238        blt::vector destroy $vector
239    }
240}
241
242# ----------------------------------------------------------------------
243# USAGE: components ?-name|-dimensions|-style? ?<pattern>?
244#
245# Returns a list of names or types for the various components of
246# this field.  If the optional glob-style <pattern> is specified,
247# then it returns only the components with names matching the pattern.
248# ----------------------------------------------------------------------
249itcl::body Rappture::Field::components {args} {
250    Rappture::getopts args params {
251        flag what -name default
252        flag what -dimensions
253        flag what -style
254        flag what -particles
255        flag what -flow
256        flag what -box
257    }
258
259    set pattern *
260    if {[llength $args] > 0} {
261        set pattern [lindex $args 0]
262        set args [lrange $args 1 end]
263    }
264    if {[llength $args] > 0} {
265        error "wrong # args: should be \"components ?switches? ?pattern?\""
266    }
267
268    # There's only one dimension of the field.  Components can't have
269    # different dimensions in the same field.  They would by definition be
270    # using different meshes and viewers.
271    if { $params(what) == "-dimensions" } {
272        return "${_dim}D"
273    }
274    # BE CAREFUL: return component names in proper order
275    set rlist ""
276    set components {}
277    # First compile a list of valid components that match the pattern
278    foreach cname [$_field children -type component] {
279        if { ![info exists _isValidComponent($cname)] } {
280            continue
281        }
282        if { [string match $pattern $cname] } {
283            lappend components $cname
284        }
285    }
286    # Now handle the tests.
287    switch -- $params(what) {
288        -name {
289            set rlist $components
290        }
291        -style {
292            foreach cname $components {
293                if { [info exists _comp2style($cname)] } {
294                    lappend rlist $_comp2style($cname)
295                }
296            }
297        }
298    }
299    return $rlist
300}
301
302# ----------------------------------------------------------------------
303# USAGE: mesh ?<name>?
304#
305# Returns a list {xvec yvec} for the specified field component <name>.
306# If the name is not specified, then it returns the vectors for the
307# overall field (sum of all components).
308# ----------------------------------------------------------------------
309itcl::body Rappture::Field::mesh {{cname -overall}} {
310    if {$cname == "-overall" || $cname == "component0"} {
311        set cname [lindex [components -name] 0]
312    }
313    if {[info exists _comp2xy($cname)]} {
314        return [lindex $_comp2xy($cname) 0]  ;# return xv
315    }
316    if { [info exists _comp2vtk($cname)] } {
317        # FIXME: extract mesh from VTK file data.
318        if { $_comp2dims($cname) == "1D" } {
319            return $_xv
320        }
321        error "method \"mesh\" is not implemented for VTK file data"
322    }
323    if {[info exists _comp2dx($cname)]} {
324        return ""  ;# no mesh -- it's embedded in the blob data
325    }
326    if {[info exists _comp2mesh($cname)]} {
327        return ""  ;# no mesh -- it's embedded in the value data
328    }
329    if {[info exists _comp2unirect2d($cname)]} {
330        set mobj [lindex $_comp2unirect2d($cname) 0]
331        return [$mobj mesh]
332    }
333    if {[info exists _comp2unirect3d($cname)]} {
334        set mobj [lindex $_comp2unirect3d($cname) 0]
335        return [$mobj mesh]
336    }
337    error "can't get field mesh: Unknown component \"$cname\": should be one of [join [lsort [array names _comp2dims]] {, }]"
338}
339
340# ----------------------------------------------------------------------
341# USAGE: values ?<name>?
342#
343# Returns a list {xvec yvec} for the specified field component <name>.
344# If the name is not specified, then it returns the vectors for the
345# overall field (sum of all components).
346# ----------------------------------------------------------------------
347itcl::body Rappture::Field::values {cname} {
348    if {$cname == "component0"} {
349        set cname "component"
350    }
351    if {[info exists _comp2xy($cname)]} {
352        return [lindex $_comp2xy($cname) 1]  ;# return yv
353    }
354    # VTK file data
355    if { [info exists _comp2vtk($cname)] } {
356        # FIXME: extract the values from the VTK file data
357        if { $_comp2dims($cname) == "1D" } {
358            return $_values
359        }
360        error "method \"values\" is not implemented for vtk file data"
361    }
362    # Points-on-mesh
363    if { [info exists _comp2mesh($cname)] } {
364        set vector [lindex $_comp2mesh($cname) 1]
365        return [$vector range 0 end]
366    }
367    if {[info exists _comp2dx($cname)]} {
368        error "method \"values\" is not implemented for dx file data"
369    }
370    if {[info exists _comp2unirect2d($cname)]} {
371        return [$_comp2unirect2d($cname) values]
372    }
373    if {[info exists _comp2unirect3d($cname)]} {
374        return [$_comp2unirect3d($cname) blob]
375    }
376    error "can't get field values. Unknown component \"$cname\": should be [join [lsort [array names _comp2dims]] {, }]"
377}
378
379# ----------------------------------------------------------------------
380# USAGE: blob ?<name>?
381#
382# Returns a string representing the blob of data for the mesh and values.
383# ----------------------------------------------------------------------
384itcl::body Rappture::Field::blob {cname} {
385    if {$cname == "component0"} {
386        set cname "component"
387    }
388    if {[info exists _comp2xy($cname)]} {
389        return ""
390    }
391    if { [info exists _comp2vtk($cname)] } {
392        error "blob not implemented for VTK file data"
393    }
394    if {[info exists _comp2dx($cname)]} {
395        return $_comp2dx($cname)  ;# return gzipped, base64-encoded DX data
396    }
397    if {[info exists _comp2unirect2d($cname)]} {
398        set blob [$_comp2unirect2d($cname) blob]
399        lappend blob "values" $_values
400        return $blob
401    }
402    if {[info exists _comp2unirect3d($cname)]} {
403        return [$_comp2unirect3d($cname) blob]
404    }
405    error "can't get field blob: Unknown component \"$cname\": should be one of [join [lsort [array names _comp2dims]] {, }]"
406}
407
408# ----------------------------------------------------------------------
409# USAGE: valueLimits <cname>
410#
411# Returns an array for the requested component with a list {min max}
412# representing the limits for each axis.
413# ----------------------------------------------------------------------
414itcl::body Rappture::Field::valueLimits { cname } {
415    if { [info exists _comp2limits($cname)] } {
416        return $_comp2limits($cname)
417    }
418    return ""
419}
420
421# ----------------------------------------------------------------------
422# USAGE: limits <axis>
423#
424# Returns a list {min max} representing the limits for the specified
425# axis.
426# ----------------------------------------------------------------------
427itcl::body Rappture::Field::limits {which} {
428    set min ""
429    set max ""
430    blt::vector tmp zero
431
432    foreach cname [array names _comp2dims] {
433        switch -- $_comp2dims($cname) {
434            1D {
435                switch -- $which {
436                    x - xlin {
437                        set pos 0; set log 0; set axis x
438                    }
439                    xlog {
440                        set pos 0; set log 1; set axis x
441                    }
442                    y - ylin - v - vlin {
443                        set pos 1; set log 0; set axis y
444                    }
445                    ylog - vlog {
446                        set pos 1; set log 1; set axis y
447                    }
448                    default {
449                        error "bad axis \"$which\": should be x, xlin, xlog, y, ylin, ylog, v, vlin, vlog"
450                    }
451                }
452
453                set vname [lindex $_comp2xy($cname) $pos]
454
455                if {$log} {
456                    # on a log scale, use abs value and ignore 0's
457                    $vname dup tmp
458                    $vname dup zero
459                    zero expr {tmp == 0}            ;# find the 0's
460                    tmp expr {abs(tmp)}             ;# get the abs value
461                    tmp expr {tmp + zero*max(tmp)}  ;# replace 0's with abs max
462                    set axisMin [blt::vector expr min(tmp)]
463                    set axisMax [blt::vector expr max(tmp)]
464                } 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            2D - 3D {
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    blt::vector destroy tmp zero
523    set val [$_field get "${axis}axis.min"]
524    if {"" != $val && "" != $min} {
525        if {$val > $min} {
526            # tool specified this min -- don't go any lower
527            set min $val
528        }
529    }
530    set val [$_field get "${axis}axis.max"]
531    if {"" != $val && "" != $max} {
532        if {$val < $max} {
533            # tool specified this max -- don't go any higher
534            set max $val
535        }
536    }
537    return [list $min $max]
538}
539
540# ----------------------------------------------------------------------
541# 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            set vector ""
922            if 0 {
923                set f [open /tmp/$_path.$cname.dx "w"]
924                puts -nonewline $f $contents
925                close $f
926            }
927            # This is temporary.  I put a check for this in the DxToVtk
928            # parser. 
929            if { [string range $contents  0 3] == "<DX>" } {
930                set contents [string range $contents 4 end]
931            }
932            if { [catch { Rappture::DxToVtk $contents } vtkdata] == 0 } {
933                ReadVtkDataSet $cname $vtkdata
934            } else {
935                puts stderr "Can't parse dx data: $vtkdata"
936            }
937            if 0 {
938                set f [open /tmp/$_path.$cname.vtk "w"]
939                puts -nonewline $f $vtkdata
940                close $f
941            }
942            if { $_alwaysConvertDX ||
943                 ($_viewer != "nanovis" && $_viewer != "flowvis") } {
944                set _type "vtk"
945                set _comp2vtk($cname) $vtkdata
946            } else {
947                set _type "dx"
948                set _comp2dx($cname) $data
949            }
950            set _comp2style($cname) [$_field get $cname.style]
951            if {[$_field element $cname.flow] != ""} {
952                set _comp2flowhints($cname) \
953                    [Rappture::FlowHints ::\#auto $_field $cname $_units]
954            }
955            set _dim 3
956            incr _counter
957        } elseif { $type == "dicom"} {
958            set contents [$_field get $cname.dicom]
959            if { $contents == "" } {
960                continue;               # Ignore this component
961            }
962            set viewer [$_field get "about.view"]
963            if { $viewer != "" } {
964                set _viewer $viewer
965            }
966            set vtkdata [DicomToVtk $cname $contents]
967            if { $_viewer == "" } {
968                set _viewer [expr {($_dim == 3) ? "vtkvolume" : "vtkimage"}]
969            }
970            set _comp2vtk($cname) $vtkdata
971            set _comp2style($cname) [$_field get $cname.style]
972            incr _counter
973        } elseif { $type == "ucd"} {
974            set contents [$_field get $cname.ucd]
975            if { $contents == "" } {
976                continue;               # Ignore this component
977            }
978            set vtkdata [AvsToVtk $cname $contents]
979            ReadVtkDataSet $cname $vtkdata
980            set _comp2vtk($cname) $vtkdata
981            set _comp2style($cname) [$_field get $cname.style]
982            incr _counter
983        }
984        set _isValidComponent($cname) 1
985        #puts stderr "Field $cname type is: $_type"
986    }
987    if { [array size _isValidComponent] == 0 } {
988        puts stderr "ERROR: All components of field \"$_path\" are invalid."
989        return 0
990    }
991    # Sanity check.  Verify that all components of the field have the same
992    # dimension.
993    set dim ""
994    foreach cname [array names _comp2dims] {
995        if { $dim == "" } {
996            set dim $_comp2dims($cname)
997            continue
998        }
999        if { $dim != $_comp2dims($cname) } {
1000            puts stderr "WARNING: A field can't have components of different dimensions: [join [array get _comp2dims] ,]"
1001            return 0
1002        }
1003    }
1004
1005    # FIXME: about.scalars and about.vectors are temporary.  With views
1006    #        the label and units for each field will be specified there.
1007    #
1008    # FIXME: Test that every <field><component> has the same field names,
1009    #        units, components.
1010    #
1011    # Override what we found in the VTK file with names that the user
1012    # selected.  We override the field label and units.
1013    foreach { fname label units } [$_field get about.scalars] {
1014        if { ![info exists _fld2Name($fname)] } {
1015            set _fld2Name($fname) $fname
1016            set _fld2Components($fname) 1
1017        }
1018        set _fld2Label($fname) $label
1019        set _fld2Units($fname) $units
1020    }
1021    foreach { fname label units } [$_field get about.vectors] {
1022        if { ![info exists _fld2Name($fname)] } {
1023            set _fld2Name($fname) $fname
1024            # We're just marking the field as vector (> 1) for now.
1025            set _fld2Components($fname) 3
1026        }
1027        set _fld2Label($fname) $label
1028        set _fld2Units($fname) $units
1029    }
1030    set _isValid 1
1031    return 1
1032}
1033
1034# ----------------------------------------------------------------------
1035# USAGE: _getValue <expr>
1036#
1037# Used internally to get the value for an expression <expr>.  Returns
1038# a list of the form {val parameterPath}, where val is the numeric
1039# value of the expression, and parameterPath is the XML path to the
1040# parameter representing the value, or "" if the <expr> does not
1041# depend on any parameters.
1042# ----------------------------------------------------------------------
1043itcl::body Rappture::Field::_getValue {expr} {
1044    #
1045    # First, look for the expression among the <parameter>'s
1046    # associated with the device.
1047    #
1048    set found 0
1049    foreach pcomp [$_xmlobj children parameters] {
1050        set id [$_xmlobj element -as id parameters.$pcomp]
1051        if {[string equal $id $expr]} {
1052            set val [$_xmlobj get parameters.$pcomp.current]
1053            if {"" == $val} {
1054                set val [$_xmlobj get parameters.$pcomp.default]
1055            }
1056            if {"" != $val} {
1057                set expr $val
1058                set found 1
1059                break
1060            }
1061        }
1062    }
1063    if {$found} {
1064        set pcomp "parameters.$pcomp"
1065    } else {
1066        set pcomp ""
1067    }
1068
1069    if {$_units != ""} {
1070        set expr [Rappture::Units::convert $expr \
1071            -context $_units -to $_units -units off]
1072    }
1073
1074    return [list $expr $pcomp]
1075}
1076
1077#
1078# isunirect2d  --
1079#
1080# Returns if the field is a unirect2d object. 
1081#
1082itcl::body Rappture::Field::isunirect2d { } {
1083    return [expr [array size _comp2unirect2d] > 0]
1084}
1085
1086#
1087# isunirect3d  --
1088#
1089# Returns if the field is a unirect3d object. 
1090#
1091itcl::body Rappture::Field::isunirect3d { } {
1092    return [expr [array size _comp2unirect3d] > 0]
1093}
1094
1095#
1096# flowhints  --
1097#
1098# Returns the hints associated with a flow vector field. 
1099#
1100itcl::body Rappture::Field::flowhints { cname } {
1101    if { [info exists _comp2flowhints($cname)] } {
1102        return $_comp2flowhints($cname)
1103    }
1104    return ""
1105}
1106
1107#
1108# style  --
1109#
1110# Returns the style associated with a component of the field. 
1111#
1112itcl::body Rappture::Field::style { cname } {
1113    if { [info exists _comp2style($cname)] } {
1114        return $_comp2style($cname)
1115    }
1116    return ""
1117}
1118
1119#
1120# type  --
1121#
1122# Returns the data storage type of the field.
1123#
1124# FIXME: What are the valid types?
1125#
1126itcl::body Rappture::Field::type {} {
1127    return $_type
1128}
1129
1130#
1131# numComponents --
1132#
1133# Returns the number of components in the field component.
1134#
1135itcl::body Rappture::Field::numComponents {cname} {
1136    set name $cname
1137    regsub -all { } $name {_} name
1138    if {[info exists _fld2Components($name)] } {
1139        return $_fld2Components($name)
1140    }
1141    return 1;                           # Default to scalar.
1142}
1143
1144#
1145# extents --
1146#
1147# Returns if the field is a unirect2d object. 
1148#
1149itcl::body Rappture::Field::extents {{cname -overall}} {
1150    if {$cname == "-overall" } {
1151        set max 0
1152        foreach cname [$_field children -type component] {
1153            if { ![info exists _comp2unirect3d($cname)] &&
1154                 ![info exists _comp2extents($cname)] } {
1155                continue
1156            }
1157            set value $_comp2extents($cname)
1158            if { $max < $value } {
1159                set max $value
1160            }
1161        }
1162        return $max
1163    }
1164    if { $cname == "component0"} {
1165        set cname [lindex [components -name] 0]
1166    }
1167    return $_comp2extents($cname)
1168}
1169
1170itcl::body Rappture::Field::VerifyVtkDataSet { contents } {
1171    package require vtk
1172
1173    set reader $this-datasetreader
1174    vtkDataSetReader $reader
1175
1176    # Write the contents to a file just in case it's binary.
1177    set tmpfile file[pid].vtk
1178    set f [open "$tmpfile" "w"]
1179    fconfigure $f -translation binary -encoding binary
1180    puts $f $contents
1181    close $f
1182
1183    $reader SetFileName $tmpfile
1184    $reader ReadAllNormalsOn
1185    $reader ReadAllTCoordsOn
1186    $reader ReadAllScalarsOn
1187    $reader ReadAllColorScalarsOn
1188    $reader ReadAllVectorsOn
1189    $reader ReadAllTensorsOn
1190    $reader ReadAllFieldsOn
1191    $reader Update
1192    file delete $tmpfile
1193
1194    set dataset [$reader GetOutput]
1195    set dataAttrs [$dataset GetPointData]
1196    if { $dataAttrs == ""} {
1197        puts stderr "WARNING: No point data found in \"$_path\""
1198        rename $reader ""
1199        return 0
1200    }
1201    rename $reader ""
1202}
1203
1204itcl::body Rappture::Field::ReadVtkDataSet { cname contents } {
1205    package require vtk
1206
1207    set reader $this-datasetreader
1208    vtkDataSetReader $reader
1209
1210    # Write the contents to a file just in case it's binary.
1211    set tmpfile file[pid].vtk
1212    set f [open "$tmpfile" "w"]
1213    fconfigure $f -translation binary -encoding binary
1214    puts $f $contents
1215    close $f
1216
1217    $reader SetFileName $tmpfile
1218    $reader ReadAllNormalsOn
1219    $reader ReadAllTCoordsOn
1220    $reader ReadAllScalarsOn
1221    $reader ReadAllColorScalarsOn
1222    $reader ReadAllVectorsOn
1223    $reader ReadAllTensorsOn
1224    $reader ReadAllFieldsOn
1225    $reader Update
1226    file delete $tmpfile
1227
1228    set dataset [$reader GetOutput]
1229    set limits {}
1230    foreach {xmin xmax ymin ymax zmin zmax} [$dataset GetBounds] break
1231    # Figure out the dimension of the mesh from the bounds.
1232    set _dim 0
1233    if { $xmax > $xmin } {
1234        incr _dim
1235    }
1236    if { $ymax > $ymin } {
1237        incr _dim
1238    }
1239    if { $zmax > $zmin } {
1240        incr _dim
1241    }
1242    if { $_viewer == "" } {
1243        if { $_dim == 2 } {
1244            set _viewer contour
1245        } else {
1246            set _viewer isosurface
1247        }
1248    }
1249    set _comp2dims($cname) ${_dim}D
1250    if { $_dim < 2 } {
1251        set numPoints [$dataset GetNumberOfPoints]
1252        set xv [blt::vector create \#auto]
1253        for { set i 0 } { $i < $numPoints } { incr i } {
1254            set point [$dataset GetPoint $i]
1255            $xv append [lindex $point 0]
1256        }
1257        set yv [blt::vector create \#auto]
1258        set dataAttrs [$dataset GetPointData]
1259        if { $dataAttrs == ""} {
1260            puts stderr "WARNING: No point data found in \"$_path\""
1261            rename $reader ""
1262            return 0
1263        }
1264        set array [$dataAttrs GetScalars]
1265        if { $array == ""} {
1266            puts stderr "WARNING: No scalar point data found in \"$_path\""
1267            rename $reader ""
1268            return 0
1269        }
1270        set numTuples [$array GetNumberOfTuples]
1271        for { set i 0 } { $i < $numTuples } { incr i } {
1272            $yv append [$array GetComponent $i 0]
1273        }
1274        $xv sort $yv
1275        set _comp2xy($cname) [list $xv $yv]
1276    }
1277    lappend limits x [list $xmin $xmax]
1278    lappend limits y [list $ymin $ymax]
1279    lappend limits z [list $zmin $zmax]
1280    set dataAttrs [$dataset GetPointData]
1281    if { $dataAttrs == ""} {
1282        puts stderr "WARNING: No point data found in \"$_path\""
1283        rename $reader ""
1284        return 0
1285    }
1286    set vmin 0
1287    set vmax 1
1288    set numArrays [$dataAttrs GetNumberOfArrays]
1289    if { $numArrays > 0 } {
1290        for {set i 0} {$i < [$dataAttrs GetNumberOfArrays] } {incr i} {
1291            set array [$dataAttrs GetArray $i]
1292            set fname  [$dataAttrs GetArrayName $i]
1293            foreach {min max} [$array GetRange -1] break
1294            if {$i == 0} {
1295                set vmin $min
1296                set vmax $max
1297            }
1298            lappend limits $fname [list $min $max]
1299            set _fld2Units($fname) ""
1300            set _fld2Label($fname) $fname
1301            # Let the VTK file override the <type> designated.
1302            set _fld2Components($fname) [$array GetNumberOfComponents]
1303            lappend _comp2fldName($cname) $fname
1304        }
1305    }
1306   
1307    lappend limits v [list $vmin $vmax]
1308    set _comp2limits($cname) $limits
1309
1310    rename $reader ""
1311}
1312
1313#
1314# vtkdata --
1315#
1316#       Returns a string representing the mesh and field data for a specific
1317#       component in the legacy VTK file format.
1318#
1319itcl::body Rappture::Field::vtkdata {cname} {
1320    if {$cname == "component0"} {
1321        set cname "component"
1322    }
1323    # DX: Convert DX to VTK
1324    if {[info exists _comp2dx($cname)]} {
1325        set data $_comp2dx($cname)
1326        set data [Rappture::encoding::decode $data]
1327        return [Rappture::DxToVtk $data]
1328    }
1329    # Unirect3d: isosurface
1330    if {[info exists _comp2unirect3d($cname)]} {
1331        return [$_comp2unirect3d($cname) vtkdata]
1332    }
1333    # VTK file data:
1334    if { [info exists _comp2vtk($cname)] } {
1335        return $_comp2vtk($cname)
1336    }
1337    # Points on mesh:  Construct VTK file output.
1338    if { [info exists _comp2mesh($cname)] } {
1339        # Data is in the form mesh and vector
1340        foreach {mesh vector} $_comp2mesh($cname) break
1341        set label $cname
1342        regsub -all { } $label {_} label
1343        append out "# vtk DataFile Version 3.0\n"
1344        append out "[hints label]\n"
1345        append out "ASCII\n"
1346        append out [$mesh vtkdata]
1347
1348        if { $_comp2assoc($cname) == "pointdata" } {
1349            set vtkassoc "POINT_DATA"
1350        } elseif { $_comp2assoc($cname) == "celldata" } {
1351            set vtkassoc "CELL_DATA"
1352        } elseif { $_comp2assoc($cname) == "fielddata" } {
1353            set vtkassoc "FIELD"
1354        } else {
1355            error "unknown association \"$_comp2assoc($cname)\""
1356        }
1357        set elemSize [numComponents $cname]
1358        set numValues [expr [$vector length] / $elemSize]
1359        if { $_comp2assoc($cname) == "fielddata" } {
1360            append out "$vtkassoc FieldData 1\n"
1361            append out "$label $elemSize $numValues double\n"
1362        } else {
1363            append out "$vtkassoc $numValues\n"
1364            if { $_comp2type($cname) == "colorscalars" } {
1365                # Must be float for ASCII, unsigned char for BINARY
1366                append out "COLOR_SCALARS $label $elemSize\n"
1367            } elseif { $_comp2type($cname) == "normals" } {
1368                # elemSize must equal 3
1369                append out "NORMALS $label double\n"
1370            } elseif { $_comp2type($cname) == "scalars" } {
1371                # elemSize can be 1, 2, 3 or 4
1372                append out "SCALARS $label double $elemSize\n"
1373                append out "LOOKUP_TABLE default\n"
1374            } elseif { $_comp2type($cname) == "tcoords" } {
1375                # elemSize must be 1, 2, or 3
1376                append out "TEXTURE_COORDINATES $label $elemSize double\n"
1377            } elseif { $_comp2type($cname) == "tensors" } {
1378                # elemSize must equal 9
1379                append out "TENSORS $label double\n"
1380            } elseif { $_comp2type($cname) == "vectors" } {
1381                # elemSize must equal 3
1382                append out "VECTORS $label double\n"
1383            } else {
1384                error "unknown element type \"$_comp2type($cname)\""
1385            }
1386        }
1387        append out [$vector range 0 end]
1388        append out "\n"
1389        if 0 {
1390            VerifyVtkDataSet $out
1391        }
1392        return $out
1393    }
1394    error "can't find vtkdata for $cname. This method should only be called by the vtkheightmap widget"
1395}
1396
1397#
1398# BuildPointsOnMesh --
1399#
1400#       Parses the field XML description to build a mesh and values vector
1401#       representing the field.  Right now we handle the deprecated types
1402#       of "cloud", "unirect2d", and "unirect3d" (mostly for flows).
1403#
1404itcl::body Rappture::Field::BuildPointsOnMesh {cname} {
1405    #
1406    # More complex 2D/3D data is represented by a mesh
1407    # object and an associated vector for field values.
1408    #
1409    set path [$_field get $cname.mesh]
1410    if {[$_xmlobj element $path] == ""} {
1411        # Unknown mesh designated.
1412        return 0
1413    }
1414    set viewer [$_field get "about.view"]
1415    if { $viewer != "" } {
1416        set _viewer $viewer
1417    }
1418    set element [$_xmlobj element -as type $path]
1419    set name $cname
1420    regsub -all { } $name {_} name
1421    set _fld2Label($name) $name
1422    set label [hints zlabel]
1423    if { $label != "" } {
1424        set _fld2Label($name) $label
1425    }
1426    set _fld2Units($name) [hints zunits]
1427    set _fld2Components($name) $_comp2size($cname)
1428    lappend _comp2fldName($cname) $name
1429
1430    # Handle bizarre cases that hopefully will be deprecated.
1431    if { $element == "unirect3d" } {
1432        # Special case: unirect3d (should be deprecated) + flow.
1433        if { [$_field element $cname.extents] != "" } {
1434            set vectorsize [$_field get $cname.extents]
1435        } else {
1436            set vectorsize 1
1437        }
1438        set _type unirect3d
1439        set _dim 3
1440        if { $_viewer == "" } {
1441            set _viewer flowvis
1442        }
1443        set _comp2dims($cname) "3D"
1444        set _comp2unirect3d($cname) \
1445            [Rappture::Unirect3d \#auto $_xmlobj $_field $cname $vectorsize]
1446        set _comp2style($cname) [$_field get $cname.style]
1447        set limits {}
1448        foreach axis { x y z } {
1449            lappend limits $axis [$_comp2unirect3d($cname) limits $axis]
1450        }
1451        # Get the data limits
1452        set vector [$_comp2unirect3d($cname) valuesObj]
1453        set minmax [VectorLimits $vector $vectorsize]
1454        lappend limits $cname $minmax
1455        lappend limits v      $minmax
1456        set _comp2limits($cname) $limits
1457        if {[$_field element $cname.flow] != ""} {
1458            set _comp2flowhints($cname) \
1459                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1460        }
1461        incr _counter
1462        return 1
1463    }
1464    if { $element == "unirect2d" && [$_field element $cname.flow] != "" } {
1465        # Special case: unirect2d (normally deprecated) + flow.
1466        if { [$_field element $cname.extents] != "" } {
1467            set vectorsize [$_field get $cname.extents]
1468        } else {
1469            set vectorsize 1
1470        }
1471        set _type unirect2d
1472        set _dim 2
1473        if { $_viewer == "" } {
1474            set _viewer "flowvis"
1475        }
1476        set _comp2dims($cname) "2D"
1477        set _comp2unirect2d($cname) \
1478            [Rappture::Unirect2d \#auto $_xmlobj $path]
1479        set _comp2style($cname) [$_field get $cname.style]
1480        set _comp2flowhints($cname) \
1481            [Rappture::FlowHints ::\#auto $_field $cname $_units]
1482        set _values [$_field get $cname.values]
1483        set limits {}
1484        foreach axis { x y z } {
1485            lappend limits $axis [$_comp2unirect2d($cname) limits $axis]
1486        }
1487        set xv [blt::vector create \#auto]
1488        $xv set $_values
1489        set minmax [VectorLimits $xv $vectorsize]
1490        lappend limits $cname $minmax
1491        lappend limits v $minmax
1492        blt::vector destroy $xv
1493        set _comp2limits($cname) $limits
1494        incr _counter
1495        return 1
1496    }
1497    switch -- $element {
1498        "cloud" {
1499            set mesh [Rappture::Cloud::fetch $_xmlobj $path]
1500            set _type cloud
1501        }
1502        "mesh" {
1503            set mesh [Rappture::Mesh::fetch $_xmlobj $path]
1504            set _type mesh
1505        }           
1506        "unirect2d" {
1507            if { $_viewer == "" } {
1508                set _viewer "heightmap"
1509            }
1510            set mesh [Rappture::Unirect2d::fetch $_xmlobj $path]
1511            set _type unirect2d
1512        }
1513    }
1514    if { ![$mesh isvalid] } {
1515        puts stderr "Mesh is invalid"
1516        return 0
1517    }
1518    set _dim [$mesh dimensions]
1519    if { $_dim == 3 } {
1520        set dim 0
1521        foreach axis {x y z} {
1522            foreach {min max} [$mesh limits $axis] {
1523                if { $min < $max } {
1524                    incr dim
1525                }
1526            }
1527        }
1528        if { $dim != 3 } {
1529            set _dim $dim
1530        }
1531    }
1532
1533    if {$_dim == 1} {
1534        # 1D data: Create vectors for graph widget.
1535        # Is this used anywhere?
1536        #
1537        # OOPS!  This is 1D data
1538        # Forget the cloud/field -- store BLT vectors
1539        #
1540        # Is there a natural growth path in generating output from 1D to
1541        # higher dimensions?  If there isn't, let's kill this in favor
1542        # or explicitly using a <curve> instead.  Otherwise, the features
1543        # (methods such as xmarkers) or the <curve> need to be added
1544        # to the <field>.
1545        #
1546        set xv [blt::vector create x$_counter]
1547        set yv [blt::vector create y$_counter]
1548       
1549        $yv set [$mesh points]
1550        $xv seq 0 1 [$yv length]
1551        # sort x-coords in increasing order
1552        $xv sort $yv
1553        set _comp2dims($cname) "1D"
1554        set _comp2xy($cname) [list $xv $yv]
1555        incr _counter
1556        return 1
1557    }
1558    if {$_dim == 2} {
1559        # 2D data: By default surface or contour plot using heightmap widget.
1560        set v [blt::vector create \#auto]
1561        $v set [$_field get $cname.values]
1562        if { [$v length] == 0 } {
1563            return 0
1564        }
1565        if { $_viewer == "" } {
1566            set _viewer "contour"
1567        }
1568        set numFieldValues [$v length]
1569        set numComponentsPerTuple [numComponents $cname]
1570        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1571            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1572            return 0
1573        }
1574        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1575        if { $_comp2assoc($cname) == "pointdata" } {
1576            set numPoints [$mesh numpoints]
1577            if { $numPoints != $numFieldTuples } {
1578                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1579                return 0
1580            }
1581        } elseif { $_comp2assoc($cname) == "celldata" } {
1582            set numCells [$mesh numcells]
1583            if { $numCells != $numFieldTuples } {
1584                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1585                return 0
1586            }
1587        }
1588        set _comp2dims($cname) "[$mesh dimensions]D"
1589        set _comp2mesh($cname) [list $mesh $v]
1590        set _comp2style($cname) [$_field get $cname.style]
1591        if {[$_field element $cname.flow] != ""} {
1592            set _comp2flowhints($cname) \
1593                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1594        }
1595        incr _counter
1596        array unset _comp2limits $cname
1597        foreach axis { x y z } {
1598            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1599        }
1600        set minmax [VectorLimits $v $_comp2size($cname)]
1601        lappend _comp2limits($cname) $cname $minmax
1602        lappend _comp2limits($cname) v $minmax
1603        return 1
1604    }
1605    if {$_dim == 3} {
1606        # 3D data: By default isosurfaces plot using isosurface widget.
1607        if { $_viewer == "" } {
1608            set _viewer "isosurface"
1609        }
1610        set v [blt::vector create \#auto]
1611        $v set [$_field get $cname.values]
1612        if { [$v length] == 0 } {
1613            return 0
1614        }
1615        set numFieldValues [$v length]
1616        set numComponentsPerTuple [numComponents $cname]
1617        if { [expr $numFieldValues % $numComponentsPerTuple] != 0 } {
1618            puts stderr "ERROR: Number of field values ($numFieldValues) not divisble by elemsize ($numComponentsPerTuple)"
1619            return 0
1620        }
1621        set numFieldTuples [expr $numFieldValues / $numComponentsPerTuple]
1622        if { $_comp2assoc($cname) == "pointdata" } {
1623            set numPoints [$mesh numpoints]
1624            if { $numPoints != $numFieldTuples } {
1625                puts stderr "ERROR: Number of points in mesh ($numPoints) and number of field tuples ($numFieldTuples) don't agree"
1626                return 0
1627            }
1628        } elseif { $_comp2assoc($cname) == "celldata" } {
1629            set numCells [$mesh numcells]
1630            if { $numCells != $numFieldTuples } {
1631                puts stderr "ERROR: Number of cells in mesh ($numCells) and number of field tuples ($numFieldTuples) don't agree"
1632                return 0
1633            }
1634        }
1635        set _comp2dims($cname) "[$mesh dimensions]D"
1636        set _comp2mesh($cname) [list $mesh $v]
1637        set _comp2style($cname) [$_field get $cname.style]
1638        if {[$_field element $cname.flow] != ""} {
1639            set _comp2flowhints($cname) \
1640                [Rappture::FlowHints ::\#auto $_field $cname $_units]
1641        }
1642        incr _counter
1643        foreach axis { x y z } {
1644            lappend _comp2limits($cname) $axis [$mesh limits $axis]
1645        }
1646        set minmax [VectorLimits $v $_comp2size($cname)]
1647        lappend _comp2limits($cname) $cname $minmax
1648        lappend _comp2limits($cname) v $minmax
1649        return 1
1650    }
1651    error "unhandled case in field dim=$_dim element=$element"
1652}
1653
1654itcl::body Rappture::Field::AvsToVtk { cname contents } {
1655    package require vtk
1656
1657    set reader $this-datasetreader
1658    vtkAVSucdReader $reader
1659
1660    # Write the contents to a file just in case it's binary.
1661    set tmpfile $cname[pid].ucd
1662    set f [open "$tmpfile" "w"]
1663    fconfigure $f -translation binary -encoding binary
1664    puts $f $contents
1665    close $f
1666    $reader SetFileName $tmpfile
1667    $reader Update
1668    file delete $tmpfile
1669
1670    set tmpfile $cname[pid].vtk
1671    set writer $this-datasetwriter
1672    vtkDataSetWriter $writer
1673    $writer SetInputConnection [$reader GetOutputPort]
1674    $writer SetFileName $tmpfile
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    foreach key [array names data] {
1712        if {$key == "vtkdata"} {
1713            if {1} {
1714                set f [open /tmp/$cname.vtk "w"]
1715                fconfigure $f -translation binary -encoding binary
1716                puts -nonewline $f $data(vtkdata)
1717                close $f
1718            }
1719        } else {
1720            puts stderr "$key = \"$data($key)\""
1721        }
1722    }
1723
1724    # Save viewer choice
1725    set viewer $_viewer
1726    ReadVtkDataSet $cname $data(vtkdata)
1727    # Restore viewer choice (ReadVtkDataSet wants to set it to contour/isosurface)
1728    set _viewer $viewer
1729    return $data(vtkdata)
1730}
1731
1732itcl::body Rappture::Field::GetTypeAndSize { cname } {
1733    array set type2components {
1734        "colorscalars"         4
1735        "normals"              3
1736        "scalars"              1
1737        "tcoords"              2
1738        "tensors"              9
1739        "vectors"              3
1740    }
1741    set type [$_field get $cname.elemtype]
1742    if { $type == "" } {
1743        set type "scalars"
1744    }
1745    if { ![info exists type2components($type)] } {
1746        error "unknown <elemtype> \"$type\" in field"
1747    }
1748    set size [$_field get $cname.elemsize]
1749    if { $size == "" } {
1750        set size $type2components($type)
1751    }
1752    set _comp2type($cname) $type
1753    set _comp2size($cname) $size
1754}
1755
1756itcl::body Rappture::Field::GetAssociation { cname } {
1757    set assoc [$_field get $cname.association]
1758    if { $assoc == "" } {
1759        set _comp2assoc($cname) "pointdata"
1760        return
1761    }
1762    switch -- $assoc {
1763        "pointdata" - "celldata" - "fielddata" {
1764            set _comp2assoc($cname) $assoc
1765            return
1766        }
1767        default {
1768            error "unknown field association \"$assoc\""
1769        }
1770    }
1771}
1772
1773#
1774# Compute the per-component limits or limits of vector magnitudes
1775#
1776itcl::body Rappture::Field::VectorLimits {vector vectorsize {comp -1}} {
1777    if {$vectorsize == 1} {
1778        set minmax [$vector limits]
1779    } else {
1780        set len [$vector length]
1781        if {[expr $len % $vectorsize] != 0} {
1782            error "Invalid vectorsize: $vectorsize"
1783        }
1784        if {$comp > $vectorsize-1} {
1785            error "Invalid vector component: $comp"
1786        }
1787        set numTuples [expr ($len/$vectorsize)]
1788        for {set i 0} {$i < $numTuples} {incr i} {
1789            if {$comp >= 0} {
1790                set idx [expr ($i * $vectorsize + $comp)]
1791                set val [$vector index $idx]
1792            } else {
1793                set idx [expr ($i * $vectorsize)]
1794                set mag 0
1795                for {set j 0} {$j < $vectorsize} {incr j} {
1796                    set val [$vector index $idx]
1797                    set mag [expr ($mag + $val * $val)]
1798                    incr idx
1799                }
1800                set val [expr (sqrt($mag))]
1801            }
1802            if (![info exists minmax]) {
1803                set minmax [list $val $val]
1804            } else {
1805                if {$val < [lindex $minmax 0]} {
1806                    lset minmax 0 $val
1807                }
1808                if {$val > [lindex $minmax 1]} {
1809                    lset minmax 1 $val
1810                }
1811            }
1812        }
1813    }
1814    return $minmax
1815}
Note: See TracBrowser for help on using the repository browser.