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

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

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