Context Navigation

source: nanovis/branches/1.1/dxReader.cpp @ 4818

Last change on this file since 4818 was 4612, checked in by ldelgass, 10 years ago
merge r3597 from trunk
Property svn:eol-style set to `native`
File size: 17.7 KB

Line
1	/* -- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -- */
2	/*
3	* ----------------------------------------------------------------------
4	* Nanovis: Visualization of Nanoelectronics Data
5	*
6	* dxReader.cpp
7	* This module contains openDX readers for 2D and 3D volumes.
8	*
9	* ======================================================================
10	* AUTHOR: Wei Qiao <qiaow@purdue.edu>
11	* Michael McLennan <mmclennan@purdue.edu>
12	* Purdue Rendering and Perceptualization Lab (PURPL)
13	*
14	* Copyright (c) 2004-2013 HUBzero Foundation, LLC
15	*
16	* See the file "license.terms" for information on usage and
17	* redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
18	* ======================================================================
19	*/
20	#include <stdio.h>
21	#include <math.h>
22	#include <sys/types.h>
23	#include <unistd.h>
24	#include <float.h>
25
26	#include <fstream>
27	#include <iostream>
28	#include <sstream>
29	#include <string>
30
31	#include <RpOutcome.h>
32	#include <RpField1D.h>
33	#include <RpFieldRect3D.h>
34	#include <RpFieldPrism3D.h>
35
36	// common file/data reader functions
37	#include "ReaderCommon.h"
38
39	#include "config.h"
40	#include "nanovis.h"
41	#include "Unirect.h"
42	#include "Volume.h"
43	#include "ZincBlendeVolume.h"
44	#include "NvZincBlendeReconstructor.h"
45	#ifdef USE_POINTSET_RENDERER
46	#include "PointSet.h"
47	#endif
48
49	/**
50	* \brief Load a 3D volume from a dx-format file
51	*
52	* In DX format:
53	* rank 0 means scalars,
54	* rank 1 means vectors,
55	* rank 2 means matrices,
56	* rank 3 means tensors
57	*
58	* For rank 1, shape is a single number equal to the number of dimensions.
59	* e.g. rank 1 shape 3 means a 3-component vector field
60	*/
61	Volume *
62	load_dx_volume_stream(Rappture::Outcome& result, const char *tag,
63	std::iostream& fin)
64	{
65	TRACE("Enter tag:%s", tag);
66
67	Rappture::MeshTri2D xymesh;
68	int dummy, nx, ny, nz, nxy, npts;
69	// origin
70	double x0, y0, z0;
71	// deltas (cell dimensions)
72	double dx, dy, dz;
73	// length of volume edges
74	double lx, ly, lz;
75	// temps for delta values
76	double ddx, ddy, ddz;
77	char line[128], type[128], *start;
78
79	int isrect = 1;
80
81	dx = dy = dz = 0.0;
82	x0 = y0 = z0 = 0.0; // May not have an origin line.
83	nx = ny = nz = npts = nxy = 0;
84	while (!fin.eof()) {
85	fin.getline(line, sizeof(line) - 1);
86	if (fin.fail()) {
87	result.error("error in data stream");
88	return NULL;
89	}
90	for (start = line; start == ' ' \|\| start == '\t'; start++)
91	; // skip leading blanks
92
93	if (*start != '#') { // skip comment lines
94	if (sscanf(start, "object %d class gridpositions counts %d %d %d",
95	&dummy, &nx, &ny, &nz) == 4) {
96	// found grid size
97	isrect = 1;
98	} else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows",
99	&dummy, &nxy) == 2) {
100	isrect = 0;
101
102	double xx, yy, zz;
103	for (int i = 0; i < nxy; i++) {
104	fin.getline(line,sizeof(line)-1);
105	if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
106	xymesh.addNode(Rappture::Node2D(xx, yy));
107	}
108	}
109
110	char fpts[128];
111	sprintf(fpts, "/tmp/tmppts%d", getpid());
112	char fcells[128];
113	sprintf(fcells, "/tmp/tmpcells%d", getpid());
114
115	std::ofstream ftmp(fpts);
116	// save corners of bounding box first, to work around meshing
117	// problems in voronoi utility
118	int numBoundaryPoints = 4;
119
120	// Bump out the bounds by an epsilon to avoid problem when
121	// corner points are already nodes
122	double XEPS = (xymesh.rangeMax(Rappture::xaxis) -
123	xymesh.rangeMin(Rappture::xaxis)) / 10.0f;
124
125	double YEPS = (xymesh.rangeMax(Rappture::yaxis) -
126	xymesh.rangeMin(Rappture::yaxis)) / 10.0f;
127
128	ftmp << xymesh.rangeMin(Rappture::xaxis) - XEPS << " "
129	<< xymesh.rangeMin(Rappture::yaxis) - YEPS << std::endl;
130	ftmp << xymesh.rangeMax(Rappture::xaxis) + XEPS << " "
131	<< xymesh.rangeMin(Rappture::yaxis) - YEPS << std::endl;
132	ftmp << xymesh.rangeMax(Rappture::xaxis) + XEPS << " "
133	<< xymesh.rangeMax(Rappture::yaxis) + YEPS << std::endl;
134	ftmp << xymesh.rangeMin(Rappture::xaxis) - XEPS << " "
135	<< xymesh.rangeMax(Rappture::yaxis) + YEPS << std::endl;
136
137	for (int i = 0; i < nxy; i++) {
138	ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
139	}
140	ftmp.close();
141
142	char cmdstr[512];
143	sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
144	if (system(cmdstr) == 0) {
145	int cx, cy, cz;
146	std::ifstream ftri(fcells);
147	while (!ftri.eof()) {
148	ftri.getline(line,sizeof(line)-1);
149	if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
150	if (cx >= numBoundaryPoints &&
151	cy >= numBoundaryPoints &&
152	cz >= numBoundaryPoints) {
153	// skip boundary points we added
154	xymesh.addCell(cx - numBoundaryPoints,
155	cy - numBoundaryPoints,
156	cz - numBoundaryPoints);
157	}
158	}
159	}
160	ftri.close();
161	} else {
162	result.error("triangularization failed");
163	return NULL;
164	}
165	unlink(fpts);
166	unlink(fcells);
167	} else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
168	// found z-grid
169	} else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
170	// found origin
171	} else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
172	// found one of the delta lines
173	int count = 0;
174	if (ddx != 0.0) {
175	dx = ddx;
176	count++;
177	}
178	if (ddy != 0.0) {
179	dy = ddy;
180	count++;
181	}
182	if (ddz != 0.0) {
183	dz = ddz;
184	count++;
185	}
186	if (count > 1) {
187	result.addError("don't know how to handle multiple non-zero"
188	" delta values");
189	return NULL;
190	}
191	} else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows",
192	&dummy, type, &npts) == 3) {
193	if (isrect && (npts != nxnynz)) {
194	result.addError("inconsistent data: expected %d points"
195	" but found %d points", nxnynz, npts);
196	return NULL;
197	} else if (!isrect && (npts != nxy*nz)) {
198	result.addError("inconsistent data: expected %d points"
199	" but found %d points", nxy*nz, npts);
200	return NULL;
201	}
202	break;
203	} else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows",
204	&dummy, type, &npts) == 3) {
205	if (npts != nxnynz) {
206	result.addError("inconsistent data: expected %d points"
207	" but found %d points", nxnynz, npts);
208	return NULL;
209	}
210	break;
211	}
212	}
213	}
214
215	TRACE("found nx=%d ny=%d nxy=%d nz=%d\ndx=%f dy=%f dz=%f\nx0=%f y0=%f z0=%f",
216	nx, ny, nxy, nz, dx, dy, dz, x0, y0, z0);
217
218	lx = (nx - 1) * dx;
219	ly = (ny - 1) * dy;
220	lz = (nz - 1) * dz;
221
222	// read data points
223	if (fin.eof()) {
224	result.error("data not found in stream");
225	return NULL;
226	}
227	Volume *volume = NULL;
228	float *data = NULL;
229	double vmin = DBL_MAX;
230	double nzero_min = DBL_MAX;
231	double vmax = -DBL_MAX;
232	if (isrect) {
233	#ifdef DOWNSAMPLE_DATA
234	Rappture::Mesh1D xgrid(x0, x0 + lx, nx);
235	Rappture::Mesh1D ygrid(y0, y0 + ly, ny);
236	Rappture::Mesh1D zgrid(z0, z0 + lz, nz);
237	Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
238	#else // !DOWNSAMPLE_DATA
239	data = new float[nx * ny * nz * 4];
240	memset(data, 0, nxnynz*4);
241	#endif // DOWNSAMPLE_DATA
242	double dval[6];
243	int nread = 0;
244	int ix = 0;
245	int iy = 0;
246	int iz = 0;
247
248	while (!fin.eof() && nread < npts) {
249	fin.getline(line,sizeof(line)-1);
250	if (fin.fail()) {
251	result.addError("error reading data points");
252	return NULL;
253	}
254	int n = sscanf(line, "%lg %lg %lg %lg %lg %lg",
255	&dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
256	for (int p = 0; p < n; p++) {
257	#ifdef notdef
258	if (isnan(dval[p])) {
259	TRACE("Found NAN in input at %d,%d,%d", ix, iy, iz);
260	}
261	#endif
262	#ifdef DOWNSAMPLE_DATA
263	int nindex = iznxny + iy*nx + ix;
264	field.define(nindex, dval[p]);
265	#else // !DOWNSAMPLE_DATA
266	int nindex = (iznxny + iynx + ix) 4;
267	data[nindex] = dval[p];
268
269	if (dval[p] < vmin) {
270	vmin = dval[p];
271	} else if (dval[p] > vmax) {
272	vmax = dval[p];
273	}
274	if (dval[p] != 0.0 && dval[p] < nzero_min) {
275	nzero_min = dval[p];
276	}
277	#endif // DOWNSAMPLE_DATA
278	nread++;
279	if (++iz >= nz) {
280	iz = 0;
281	if (++iy >= ny) {
282	iy = 0;
283	++ix;
284	}
285	}
286	}
287	}
288
289	// make sure that we read all of the expected points
290	if (nread != npts) {
291	result.addError("inconsistent data: expected %d points"
292	" but found %d points", npts, nread);
293	return NULL;
294	}
295
296	#ifndef DOWNSAMPLE_DATA
297
298	// scale all values [0-1], -1 => out of bounds
299	normalizeScalar(data, npts, 4, vmin, vmax);
300	computeSimpleGradient(data, nx, ny, nz,
301	dx, dy, dz);
302	#else // DOWNSAMPLE_DATA
303	// figure out a good mesh spacing
304	int nsample = 64;
305	double dmin = pow((lxlylz)/((nsample-1)(nsample-1)(nsample-1)), 0.333);
306
307	nx = (int)ceil(lx/dmin) + 1;
308	ny = (int)ceil(ly/dmin) + 1;
309	nz = (int)ceil(lz/dmin) + 1;
310
311	#ifndef HAVE_NPOT_TEXTURES
312	// must be an even power of 2 for older cards
313	nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
314	ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
315	nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
316	#endif
317
318	dx = lx /(double)(nx - 1);
319	dy = ly /(double)(ny - 1);
320	dz = lz /(double)(nz - 1);
321
322	vmin = field.valueMin();
323	vmax = field.valueMax();
324	nzero_min = DBL_MAX;
325	double dv = vmax - vmin;
326	if (dv == 0.0) {
327	dv = 1.0;
328	}
329
330	int ngen = 0;
331	#ifdef FILTER_GRADIENTS
332	// Sobel filter expects a single float at
333	// each node
334	const int step = 1;
335	float cdata = new float[nxnynz step];
336	#else // !FILTER_GRADIENTS
337	// Leave 3 floats of space for gradients
338	// to be filled in by computeSimpleGradient()
339	const int step = 4;
340	data = new float[nxnynz * step];
341	#endif // FILTER_GRADIENTS
342
343	for (int iz = 0; iz < nz; iz++) {
344	double zval = z0 + iz*dz;
345	for (int iy = 0; iy < ny; iy++) {
346	double yval = y0 + iy*dy;
347	for (int ix = 0; ix < nx; ix++) {
348	double xval = x0 + ix*dx;
349	double v = field.value(xval, yval, zval);
350	#ifdef notdef
351	if (isnan(v)) {
352	TRACE("NAN at %d,%d,%d: (%g, %g, %g)", ix, iy, iz, xval, yval, zval);
353	}
354	#endif
355	if (v != 0.0 && v < nzero_min) {
356	nzero_min = v;
357	}
358	#ifdef FILTER_GRADIENTS
359	// NOT normalized, may have NaNs (FIXME: how does gradient filter handle NaN?)
360	cdata[ngen] = v;
361	#else // !FILTER_GRADIENTS
362	// scale all values [0-1], -1 => out of bounds
363	v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
364	data[ngen] = v;
365	#endif // FILTER_GRADIENTS
366	ngen += step;
367	}
368	}
369	}
370	#ifdef FILTER_GRADIENTS
371	// computeGradient returns a new array with gradients
372	// filled in, so data is now 4 floats per node
373	data = computeGradient(cdata, nx, ny, nz,
374	dx, dy, dz,
375	vmin, vmax);
376	delete [] cdata;
377	#else // !FILTER_GRADIENTS
378	computeSimpleGradient(data, nx, ny, nz,
379	dx, dy, dz);
380	#endif // FILTER_GRADIENTS
381
382	#endif // DOWNSAMPLE_DATA
383	} else {
384	Rappture::Mesh1D zgrid(z0, z0 + (nz-1)*dz, nz);
385	Rappture::FieldPrism3D field(xymesh, zgrid);
386
387	double dval;
388	int nread = 0;
389	int ixy = 0;
390	int iz = 0;
391	while (!fin.eof() && nread < npts) {
392	fin >> dval;
393	if (fin.fail()) {
394	result.addError("after %d of %d points: can't read number",
395	nread, npts);
396	return NULL;
397	} else {
398	int nid = nxy*iz + ixy;
399	field.define(nid, dval);
400
401	nread++;
402	if (++iz >= nz) {
403	iz = 0;
404	ixy++;
405	}
406	}
407	}
408
409	// make sure that we read all of the expected points
410	if (nread != npts) {
411	result.addError("inconsistent data: expected %d points"
412	" but found %d points", npts, nread);
413	return NULL;
414	}
415
416	x0 = field.rangeMin(Rappture::xaxis);
417	lx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
418	y0 = field.rangeMin(Rappture::yaxis);
419	ly = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
420	z0 = field.rangeMin(Rappture::zaxis);
421	lz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
422
423	// figure out a good mesh spacing
424	int nsample = 64;
425	double dmin = pow((lxlylz)/((nsample-1)(nsample-1)(nsample-1)), 0.333);
426
427	nx = (int)ceil(lx/dmin) + 1;
428	ny = (int)ceil(ly/dmin) + 1;
429	nz = (int)ceil(lz/dmin) + 1;
430	#ifndef HAVE_NPOT_TEXTURES
431	// must be an even power of 2 for older cards
432	nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
433	ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
434	nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
435	#endif
436
437	dx = lx /(double)(nx - 1);
438	dy = ly /(double)(ny - 1);
439	dz = lz /(double)(nz - 1);
440
441	vmin = field.valueMin();
442	vmax = field.valueMax();
443	nzero_min = DBL_MAX;
444	double dv = vmax - vmin;
445	if (dv == 0.0) {
446	dv = 1.0;
447	}
448
449	data = new float[4nxny*nz];
450	// generate the uniformly sampled data that we need for a volume
451	int ngen = 0;
452	for (int iz = 0; iz < nz; iz++) {
453	double zval = z0 + iz*dz;
454	for (int iy = 0; iy < ny; iy++) {
455	double yval = y0 + iy*dy;
456	for (int ix = 0; ix < nx; ix++) {
457	double xval = x0 + ix*dx;
458	double v = field.value(xval, yval, zval);
459
460	if (v != 0.0 && v < nzero_min) {
461	nzero_min = v;
462	}
463	// scale all values [0-1], -1 => out of bounds
464	v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
465	data[ngen] = v;
466
467	ngen += 4;
468	}
469	}
470	}
471
472	computeSimpleGradient(data, nx, ny, nz,
473	dx, dy, dz);
474	}
475
476	TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
477	TRACE("x0 = %lg y0 = %lg z0 = %lg", x0, y0, z0);
478	TRACE("lx = %lg ly = %lg lz = %lg", lx, ly, lz);
479	TRACE("dx = %lg dy = %lg dz = %lg", dx, dy, dz);
480	TRACE("dataMin = %lg dataMax = %lg nzero_min = %lg",
481	vmin, vmax, nzero_min);
482
483	volume = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
484	vmin, vmax, nzero_min);
485	volume->xAxis.setRange(x0, x0 + lx);
486	volume->yAxis.setRange(y0, y0 + ly);
487	volume->zAxis.setRange(z0, z0 + lz);
488	volume->updatePending = true;
489
490	// TBD..
491	#if 0 && defined(USE_POINTSET_RENDERER)
492	PointSet *pset = new PointSet();
493	pset->initialize(volume, (float*)data);
494	pset->setVisible(true);
495	NanoVis::pointSet.push_back(pset);
496	updateColor(pset);
497	volume->pointsetIndex = NanoVis::pointSet.size() - 1;
498	#endif
499	delete [] data;
500
501	//
502	// Center this new volume on the origin.
503	//
504	float dx0 = -0.5;
505	float dy0 = -0.5*ly/lx;
506	float dz0 = -0.5*lz/lx;
507	if (volume) {
508	volume->location(vrmath::Vector3f(dx0, dy0, dz0));
509	TRACE("Set volume location to %g %g %g", dx0, dy0, dz0);
510	}
511	return volume;
512	}

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

Download in other formats: