Context Navigation

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

Last change on this file since 5330 was 5330, checked in by ldelgass, 9 years ago
Remove old client side VTK molecule viewer. It is incomplete and doesn't work with VTK 6.x anyway. This wasn't being used as a fallback in any case because there is always a default pymol server (localhost).
File size: 61.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: