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

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

Add VTK download from nano/flowvisviewer

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