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

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

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