source: trunk/vizservers/nanovis/nanovis.cpp @ 825

Last change on this file since 825 was 825, checked in by vrinside, 14 years ago

Image Loader initialization

  • Add BMP loader in Nv.cpp

Point Renderer code added..

  • Put a data container and manage the containters with std::vector
  • Renderer 1) scale factor part should be taken into account
File size: 138.0 KB
Line 
1/*
2 * ----------------------------------------------------------------------
3 * Nanovis: Visualization of Nanoelectronics Data
4 *
5 * ======================================================================
6 *  AUTHOR:  Wei Qiao <qiaow@purdue.edu>
7 *           Michael McLennan <mmclennan@purdue.edu>
8 *           Purdue Rendering and Perceptualization Lab (PURPL)
9 *
10 *  Copyright (c) 2004-2006  Purdue Research Foundation
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
17#include <getopt.h>
18#include <stdio.h>
19#include <math.h>
20#include <fstream>
21#include <iostream>
22#include <sstream>
23#include <string>
24#include <sys/time.h>
25#include <sys/types.h>
26#include <unistd.h>
27#include <fcntl.h>
28#include <signal.h>
29
30#include "Nv.h"
31#include "PointSetRenderer.h"
32#include "PointSet.h"
33
34#include "nanovis.h"
35#include "RpField1D.h"
36#include "RpFieldRect3D.h"
37#include "RpFieldPrism3D.h"
38#include "RpEncode.h"
39
40//transfer function headers
41#include "transfer-function/TransferFunctionMain.h"
42#include "transfer-function/ControlPoint.h"
43#include "transfer-function/TransferFunctionGLUTWindow.h"
44#include "transfer-function/ColorGradientGLUTWindow.h"
45#include "transfer-function/ColorPaletteWindow.h"
46#include "transfer-function/MainWindow.h"
47#include "ZincBlendeVolume.h"
48#include "NvLoadFile.h"
49#include "NvColorTableRenderer.h"
50#include "NvEventLog.h"
51#include "NvZincBlendeReconstructor.h"
52#include "HeightMap.h"
53#include "Grid.h"
54
55//#define  _LOCAL_ZINC_TEST_
56
57// R2 headers
58#include <R2/R2FilePath.h>
59#include <R2/R2Fonts.h>
60
61//render server
62
63extern VolumeRenderer* g_vol_render;
64extern PointSetRenderer* g_pointset_renderer;
65extern NvColorTableRenderer* g_color_table_renderer;
66extern Grid* g_grid;
67PlaneRenderer* plane_render;
68Camera* cam;
69
70//if true nanovis renders volumes in 3D, if not renders 2D plane
71bool volume_mode = true;
72
73// color table for built-in transfer function editor
74float color_table[256][4];     
75
76// default transfer function
77char *def_transfunc = "transfunc define default {\n\
78  0.0  1 1 1\n\
79  0.2  1 1 0\n\
80  0.4  0 1 0\n\
81  0.6  0 1 1\n\
82  0.8  0 0 1\n\
83  1.0  1 0 1\n\
84} {\n\
85  0.00  1.0\n\
86  0.05  0.0\n\
87  0.15  0.0\n\
88  0.20  1.0\n\
89  0.25  0.0\n\
90  0.35  0.0\n\
91  0.40  1.0\n\
92  0.45  0.0\n\
93  0.55  0.0\n\
94  0.60  1.0\n\
95  0.65  0.0\n\
96  0.75  0.0\n\
97  0.80  1.0\n\
98  0.85  0.0\n\
99  0.95  0.0\n\
100  1.00  1.0\n\
101}";
102
103/*
104#ifdef XINETD
105FILE* xinetd_log;
106#endif
107
108FILE* event_log;
109//log
110void init_event_log();
111void end_event_log();
112double cur_time;        //in seconds
113double get_time_interval();
114*/
115
116int render_window;              //the handle of the render window;
117bool axis_on = true;
118
119// forward declarations
120//void init_particles();
121void get_slice_vectors();
122Rappture::Outcome load_volume_stream(int index, std::iostream& fin);
123Rappture::Outcome load_volume_stream2(int index, std::iostream& fin);
124void load_volume(int index, int width, int height, int depth, int n_component, float* data, double vmin, double vmax,
125                double nzero_min);
126TransferFunction* get_transfunc(char *name);
127void resize_offscreen_buffer(int w, int h);
128void offscreen_buffer_capture();
129void bmp_header_add_int(unsigned char* header, int& pos, int data);
130void bmp_write(const char* cmd);
131void bmp_write_to_file();
132void display();
133void display_offscreen_buffer();
134void read_screen();
135
136//ParticleSystem* psys;
137//float psys_x=0.4, psys_y=0, psys_z=0;
138
139Lic* lic;
140
141//frame buffer for final rendering
142unsigned char* screen_buffer = NULL;
143NVISid final_fbo, final_color_tex, final_depth_rb;
144
145//bool advect=false;
146float vert[NMESH*NMESH*3];              //particle positions in main memory
147float slice_vector[NMESH*NMESH*4];      //per slice vectors in main memory
148
149int n_volumes = 0;
150// pointers to volumes, currently handle up to 10 volumes
151vector<Volume*> volume;
152
153vector<HeightMap*> g_heightMap;
154vector<PointSet*> g_pointSet;
155
156// maps transfunc name to TransferFunction object
157Tcl_HashTable tftable;
158
159// pointers to 2D planes, currently handle up 10
160Texture2D* plane[10];
161
162// value indicates "up" axis:  x=1, y=2, z=3, -x=-1, -y=-2, -z=-3
163int updir = 2;
164
165PerfQuery* perf;                        //perfromance counter
166
167//Nvidia CG shaders and their parameters
168//INSOO
169//CGcontext g_context;
170
171CGprogram m_passthru_fprog;
172CGparameter m_passthru_scale_param, m_passthru_bias_param;
173
174extern R2Fonts* g_fonts;
175
176/*
177CGprogram m_copy_texcoord_fprog;
178CGprogram m_one_volume_fprog;
179CGparameter m_vol_one_volume_param;
180CGparameter m_tf_one_volume_param;
181CGparameter m_mvi_one_volume_param;
182CGparameter m_mv_one_volume_param;
183CGparameter m_render_param_one_volume_param;
184*/
185
186/*
187CGprogram m_vert_std_vprog;
188CGparameter m_mvp_vert_std_param;
189CGparameter m_mvi_vert_std_param;
190*/
191
192using namespace std;
193
194#define RM_VOLUME 1
195#define RM_POINT 2
196
197int renderMode = RM_VOLUME;
198
199// Tcl interpreter for incoming messages
200static Tcl_Interp *interp;
201static Tcl_DString cmdbuffer;
202
203static int ScreenShotCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
204static int CameraCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
205static int CutplaneCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
206static int LegendCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
207static int ScreenCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
208static int TransfuncCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
209static int UpCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
210static int VolumeCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
211
212static int PlaneNewCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
213static int PlaneLinkCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
214static int PlaneEnableCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
215
216static int GridCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
217static int AxisCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]));
218
219static int GetVolumeIndices _ANSI_ARGS_((Tcl_Interp *interp, int argc, CONST84 char *argv[], vector<int>* vectorPtr));
220static int GetIndices(Tcl_Interp *interp, int argc, CONST84 char *argv[], vector<int>* vectorPtr);
221static int GetAxis _ANSI_ARGS_((Tcl_Interp *interp, char *str, int *valPtr));
222static int GetColor _ANSI_ARGS_((Tcl_Interp *interp, char *str, float *rgbPtr));
223
224/*
225 * ----------------------------------------------------------------------
226 * CLIENT COMMAND:
227 *   camera aim <x0> <y0> <z0>
228 *   camera angle <xAngle> <yAngle> <zAngle>
229 *   camera zoom <factor>
230 *
231 * Clients send these commands to manipulate the camera.  The "angle"
232 * operation controls the angle of the camera around the focal point.
233 * The "zoom" operation sets the zoom factor, moving the camera in
234 * and out.
235 * ----------------------------------------------------------------------
236 */
237static int
238CameraCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
239{
240        if (argc < 2) {
241                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
242                        " option arg arg...\"", (char*)NULL);
243                return TCL_ERROR;
244    }
245
246    char c = *argv[1];
247        if (c == 'a' && strcmp(argv[1],"angle") == 0) {
248        if (argc != 5) {
249                    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
250                            " angle xangle yangle zangle\"", (char*)NULL);
251                    return TCL_ERROR;
252        }
253
254        double xangle, yangle, zangle;
255            if (Tcl_GetDouble(interp, argv[2], &xangle) != TCL_OK) {
256                    return TCL_ERROR;
257            }
258            if (Tcl_GetDouble(interp, argv[3], &yangle) != TCL_OK) {
259                    return TCL_ERROR;
260            }
261            if (Tcl_GetDouble(interp, argv[4], &zangle) != TCL_OK) {
262                    return TCL_ERROR;
263            }
264            cam->rotate(xangle, yangle, zangle);
265
266            return TCL_OK;
267        }
268        else if (c == 'a' && strcmp(argv[1],"aim") == 0) {
269        if (argc != 5) {
270                    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
271                            " aim x y z\"", (char*)NULL);
272                    return TCL_ERROR;
273        }
274
275        double x0, y0, z0;
276            if (Tcl_GetDouble(interp, argv[2], &x0) != TCL_OK) {
277                    return TCL_ERROR;
278            }
279            if (Tcl_GetDouble(interp, argv[3], &y0) != TCL_OK) {
280                    return TCL_ERROR;
281            }
282            if (Tcl_GetDouble(interp, argv[4], &z0) != TCL_OK) {
283                    return TCL_ERROR;
284            }
285            cam->aim(x0, y0, z0);
286
287            return TCL_OK;
288        }
289        else if (c == 'z' && strcmp(argv[1],"zoom") == 0) {
290        if (argc != 3) {
291                    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
292                            " zoom factor\"", (char*)NULL);
293                    return TCL_ERROR;
294        }
295
296        double zoom;
297            if (Tcl_GetDouble(interp, argv[2], &zoom) != TCL_OK) {
298                    return TCL_ERROR;
299            }
300
301        live_obj_z = -2.5/zoom;
302                cam->move(live_obj_x, live_obj_y, live_obj_z);
303
304            return TCL_OK;
305    }
306
307        Tcl_AppendResult(interp, "bad option \"", argv[1],
308                "\": should be aim, angle, or zoom", (char*)NULL);
309        return TCL_ERROR;
310}
311
312static int
313ScreenShotCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
314{
315    int old_win_width = win_width;
316    int old_win_height = win_height;
317
318#ifdef XINETD
319    resize_offscreen_buffer(1024, 1024);
320    cam->set_screen_size(30, 90, 1024 - 60, 1024 - 120);
321    offscreen_buffer_capture();  //enable offscreen render
322    display();
323
324    // INSOO
325    // TBD
326    Volume* vol = volume[0];
327    TransferFunction* tf = g_vol_render->get_volume_shading(vol);
328    if (tf)
329    {
330        float data[512];
331        for (int i=0; i < 256; i++) {
332            data[i] = data[i+256] = (float)(i/255.0);
333        }
334        Texture2D* plane = new Texture2D(256, 2, GL_FLOAT, GL_LINEAR, 1, data);
335        g_color_table_renderer->render(1024, 1024, plane, tf, vol->range_min(), vol->range_max());
336        delete plane;
337    }
338
339    read_screen();
340    glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
341
342    bmp_write("nv>screenshot -bytes");
343   
344    resize_offscreen_buffer(old_win_width, old_win_height);
345#endif
346
347        return TCL_OK;
348}
349
350/*
351 * ----------------------------------------------------------------------
352 * CLIENT COMMAND:
353 *   cutplane state on|off <axis> ?<volume>...?
354 *   cutplane position <relvalue> <axis> ?<volume>...?
355 *
356 * Clients send these commands to manipulate the cutplanes in one or
357 * more data volumes.  The "state" command turns a cutplane on or
358 * off.  The "position" command changes the position to a relative
359 * value in the range 0-1.  The <axis> can be x, y, or z.  These
360 * operations are applied to the volumes represented by one or more
361 * <volume> indices.  If no volumes are specified, then all volumes
362 * are updated.
363 * ----------------------------------------------------------------------
364 */
365static int
366CutplaneCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
367{
368    if (argc < 2) {
369        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
370            " option ?arg arg...?\"", (char*)NULL);
371        return TCL_ERROR;
372    }
373
374    char c = *argv[1];
375    if (c == 's' && strcmp(argv[1],"state") == 0) {
376        if (argc < 4) {
377            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
378                " state on|off axis ?volume ...? \"", (char*)NULL);
379            return TCL_ERROR;
380        }
381
382        int state;
383        if (Tcl_GetBoolean(interp, argv[2], &state) != TCL_OK) {
384            return TCL_ERROR;
385        }
386
387        int axis;
388        if (GetAxis(interp, (char*) argv[3], &axis) != TCL_OK) {
389            return TCL_ERROR;
390        }
391
392        vector<int> ivol;
393        if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
394            return TCL_ERROR;
395        }
396
397        vector<int>::iterator iter = ivol.begin();
398        while (iter != ivol.end()) {
399            if (state) {
400                volume[*iter]->enable_cutplane(axis);
401            } else {
402                volume[*iter]->disable_cutplane(axis);
403            }
404            ++iter;
405        }
406        return TCL_OK;
407    }
408    else if (c == 'p' && strcmp(argv[1],"position") == 0) {
409        if (argc < 4) {
410            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
411                " position relval axis ?volume ...? \"", (char*)NULL);
412            return TCL_ERROR;
413        }
414
415        double relval;
416        if (Tcl_GetDouble(interp, argv[2], &relval) != TCL_OK) {
417            return TCL_ERROR;
418        }
419        // keep this just inside the volume so it doesn't disappear
420        if (relval < 0.01) { relval = 0.01; }
421        if (relval > 0.99) { relval = 0.99; }
422
423        int axis;
424        if (GetAxis(interp, (char*) argv[3], &axis) != TCL_OK) {
425            return TCL_ERROR;
426        }
427
428        vector<int> ivol;
429        if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
430            return TCL_ERROR;
431        }
432
433        vector<int>::iterator iter = ivol.begin();
434        while (iter != ivol.end()) {
435            volume[*iter]->move_cutplane(axis, (float)relval);
436            ++iter;
437        }
438        return TCL_OK;
439    }
440
441    Tcl_AppendResult(interp, "bad option \"", argv[1],
442        "\": should be position or state", (char*)NULL);
443    return TCL_ERROR;
444}
445
446/*
447 * ----------------------------------------------------------------------
448 * CLIENT COMMAND:
449 *   legend <volumeIndex> <width> <height>
450 *
451 * Clients use this to generate a legend image for the specified
452 * transfer function.  The legend image is a color gradient from 0
453 * to one, drawn in the given transfer function.  The resulting image
454 * is returned in the size <width> x <height>.
455 * ----------------------------------------------------------------------
456 */
457static int
458LegendCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
459{
460    if (argc != 4) {
461        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
462            " transfunc width height\"", (char*)NULL);
463        return TCL_ERROR;
464    }
465
466    TransferFunction *tf = NULL;
467    int ivol;
468    if (Tcl_GetInt(interp, argv[1], &ivol) != TCL_OK) {
469        return TCL_ERROR;
470    }
471    if (ivol < n_volumes) {
472        tf = g_vol_render->get_volume_shading(volume[ivol]);
473    }
474    if (tf == NULL) {
475        Tcl_AppendResult(interp, "transfer function not defined for volume ", argv[1], (char*)NULL);
476        return TCL_ERROR;
477    }
478
479    int old_width = win_width;
480    int old_height = win_height;
481
482    int width, height;
483    if (Tcl_GetInt(interp, argv[2], &width) != TCL_OK) {
484        return TCL_ERROR;
485    }
486    if (Tcl_GetInt(interp, argv[3], &height) != TCL_OK) {
487        return TCL_ERROR;
488    }
489
490    plane_render->set_screen_size(width, height);
491    resize_offscreen_buffer(width, height);
492
493    // generate data for the legend
494    float data[512];
495    for (int i=0; i < 256; i++) {
496        data[i] = data[i+256] = (float)(i/255.0);
497    }
498    plane[0] = new Texture2D(256, 2, GL_FLOAT, GL_LINEAR, 1, data);
499    int index = plane_render->add_plane(plane[0], tf);
500    plane_render->set_active_plane(index);
501
502    offscreen_buffer_capture();
503    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear screen
504    plane_render->render();
505
506// INSOO
507    glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, screen_buffer);
508    //glReadPixels(0, 0, width, height, GL_BGR, GL_UNSIGNED_BYTE, screen_buffer); // INSOO's
509
510    std::ostringstream result;
511    result << "nv>legend " << argv[1];
512    result << " " << volume[ivol]->range_min();
513    result << " " << volume[ivol]->range_max();
514    bmp_write(result.str().c_str());
515    write(0, "\n", 1);
516
517    plane_render->remove_plane(index);
518    resize_offscreen_buffer(old_width, old_height);
519
520    return TCL_OK;
521}
522
523
524
525/*
526 * ----------------------------------------------------------------------
527 * CLIENT COMMAND:
528 *   screen <width> <height>
529 *
530 * Clients send this command to set the size of the rendering area.
531 * Future images are generated at the specified width/height.
532 * ----------------------------------------------------------------------
533 */
534static int
535ScreenCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
536{
537    int w, h;
538
539    if (argc != 3) {
540        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
541            " width height\"", (char*)NULL);
542        return TCL_ERROR;
543    }
544    if (Tcl_GetInt(interp, argv[1], &w) != TCL_OK) {
545        return TCL_ERROR;
546    }
547    if (Tcl_GetInt(interp, argv[2], &h) != TCL_OK) {
548        return TCL_ERROR;
549    }
550    resize_offscreen_buffer(w, h);
551
552    return TCL_OK;
553}
554
555/*
556 * ----------------------------------------------------------------------
557 * CLIENT COMMAND:
558 *   transfunc define <name> <colormap> <alphamap>
559 *     where <colormap> = { <v> <r> <g> <b> ... }
560 *           <alphamap> = { <v> <w> ... }
561 *
562 * Clients send these commands to manipulate the transfer functions.
563 * ----------------------------------------------------------------------
564 */
565static int
566TransfuncCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
567{
568        if (argc < 2) {
569                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
570                        " option arg arg...\"", (char*)NULL);
571                return TCL_ERROR;
572    }
573
574    char c = *argv[1];
575        if (c == 'd' && strcmp(argv[1],"define") == 0) {
576        if (argc != 5) {
577                    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
578                            argv[1], " define name colormap alphamap\"", (char*)NULL);
579            return TCL_ERROR;
580        }
581
582        // decode the data and store in a series of fields
583        Rappture::Field1D rFunc, gFunc, bFunc, wFunc;
584        int cmapc, wmapc, i, j;
585        char **cmapv, **wmapv;
586
587        if (Tcl_SplitList(interp, argv[3], &cmapc, (const char***)&cmapv) != TCL_OK) {
588            return TCL_ERROR;
589        }
590        if (cmapc % 4 != 0) {
591            Tcl_Free((char*)cmapv);
592                    Tcl_AppendResult(interp, "bad colormap in transfunc: should be ",
593                "{ v r g b ... }", (char*)NULL);
594            return TCL_ERROR;
595        }
596
597        if (Tcl_SplitList(interp, argv[4], &wmapc, (const char***)&wmapv) != TCL_OK) {
598            return TCL_ERROR;
599        }
600        if (wmapc % 2 != 0) {
601            Tcl_Free((char*)cmapv);
602            Tcl_Free((char*)wmapv);
603                    Tcl_AppendResult(interp, "bad alphamap in transfunc: should be ",
604                "{ v w ... }", (char*)NULL);
605            return TCL_ERROR;
606        }
607
608        for (i=0; i < cmapc; i += 4) {
609            double vals[4];
610            for (j=0; j < 4; j++) {
611                if (Tcl_GetDouble(interp, cmapv[i+j], &vals[j]) != TCL_OK) {
612                    Tcl_Free((char*)cmapv);
613                    Tcl_Free((char*)wmapv);
614                    return TCL_ERROR;
615                }
616                if (vals[j] < 0 || vals[j] > 1) {
617                    Tcl_Free((char*)cmapv);
618                    Tcl_Free((char*)wmapv);
619                            Tcl_AppendResult(interp, "bad value \"", cmapv[i+j],
620                        "\": should be in the range 0-1", (char*)NULL);
621                    return TCL_ERROR;
622                }
623            }
624            rFunc.define(vals[0], vals[1]);
625            gFunc.define(vals[0], vals[2]);
626            bFunc.define(vals[0], vals[3]);
627        }
628
629        for (i=0; i < wmapc; i += 2) {
630            double vals[2];
631            for (j=0; j < 2; j++) {
632                if (Tcl_GetDouble(interp, wmapv[i+j], &vals[j]) != TCL_OK) {
633                    Tcl_Free((char*)cmapv);
634                    Tcl_Free((char*)wmapv);
635                    return TCL_ERROR;
636                }
637                if (vals[j] < 0 || vals[j] > 1) {
638                    Tcl_Free((char*)cmapv);
639                    Tcl_Free((char*)wmapv);
640                            Tcl_AppendResult(interp, "bad value \"", wmapv[i+j],
641                        "\": should be in the range 0-1", (char*)NULL);
642                    return TCL_ERROR;
643                }
644            }
645            wFunc.define(vals[0], vals[1]);
646        }
647        Tcl_Free((char*)cmapv);
648        Tcl_Free((char*)wmapv);
649
650        // sample the given function into discrete slots
651        const int nslots = 256;
652        float data[4*nslots];
653        for (i=0; i < nslots; i++) {
654            double xval = double(i)/(nslots-1);
655            data[4*i]   = rFunc.value(xval);
656            data[4*i+1] = gFunc.value(xval);
657            data[4*i+2] = bFunc.value(xval);
658            data[4*i+3] = wFunc.value(xval);
659        }
660
661        // find or create this transfer function
662        int newEntry;
663        Tcl_HashEntry *entryPtr;
664        TransferFunction *tf;
665
666        entryPtr = Tcl_CreateHashEntry(&tftable, argv[2], &newEntry);
667        if (newEntry) {
668            tf = new TransferFunction(nslots, data);
669            Tcl_SetHashValue(entryPtr, (ClientData)tf);
670        } else {
671            tf = (TransferFunction*)Tcl_GetHashValue(entryPtr);
672            tf->update(data);
673        }
674
675        return TCL_OK;
676    }
677
678
679    Tcl_AppendResult(interp, "bad option \"", argv[1],
680        "\": should be define", (char*)NULL);
681    return TCL_ERROR;
682}
683
684/*
685 * ----------------------------------------------------------------------
686 * CLIENT COMMAND:
687 *   up axis
688 *
689 * Clients use this to set the "up" direction for all volumes.  Volumes
690 * are oriented such that this direction points upward.
691 * ----------------------------------------------------------------------
692 */
693static int
694UpCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
695{
696    if (argc != 2) {
697        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
698            " x|y|z|-x|-y|-z\"", (char*)NULL);
699        return TCL_ERROR;
700    }
701
702    int sign = 1;
703    char *axisName = (char*)argv[1];
704    if (*axisName == '-') {
705        sign = -1;
706        axisName++;
707    }
708
709    int axis;
710    if (GetAxis(interp, axisName, &axis) != TCL_OK) {
711        return TCL_ERROR;
712    }
713
714    updir = (axis+1)*sign;
715
716    return TCL_OK;
717}
718
719/*
720 * ----------------------------------------------------------------------
721 * CLIENT COMMAND:
722 *   volume axis label x|y|z <value> ?<volumeId> ...?
723 *   volume data state on|off ?<volumeId> ...?
724 *   volume outline state on|off ?<volumeId> ...?
725 *   volume outline color on|off ?<volumeId> ...?
726 *   volume shading transfunc <name> ?<volumeId> ...?
727 *   volume shading diffuse <value> ?<volumeId> ...?
728 *   volume shading specular <value> ?<volumeId> ...?
729 *   volume shading opacity <value> ?<volumeId> ...?
730 *   volume state on|off ?<volumeId> ...?
731 *
732 * Clients send these commands to manipulate the volumes.
733 * ----------------------------------------------------------------------
734 */
735void NvLoadVolumeBinFile2(int index, char* filename);
736static int
737VolumeCmd(ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])
738{
739    if (argc < 2) {
740        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
741            " option arg arg...\"", (char*)NULL);
742        return TCL_ERROR;
743    }
744
745    char c = *argv[1];
746    if (c == 'a' && strcmp(argv[1],"axis") == 0) {
747        if (argc < 3) {
748            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
749                argv[1], " option ?arg arg...?\"", (char*)NULL);
750            return TCL_ERROR;
751        }
752        c = *argv[2];
753        if (c == 'l' && strcmp(argv[2],"label") == 0) {
754            if (argc < 4) {
755                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
756                    argv[1], " label x|y|z string ?volume ...?\"", (char*)NULL);
757                return TCL_ERROR;
758            }
759
760            int axis;
761            if (GetAxis(interp, (char*)argv[3], &axis) != TCL_OK) {
762                return TCL_ERROR;
763            }
764
765            vector<int> ivol;
766            if (GetVolumeIndices(interp, argc-5, argv+5, &ivol) != TCL_OK) {
767                return TCL_ERROR;
768            }
769
770            vector<int>::iterator iter = ivol.begin();
771            while (iter != ivol.end()) {
772                volume[*iter]->set_label(axis, (char*)argv[4]);
773                ++iter;
774            }
775            return TCL_OK;
776        }
777
778        Tcl_AppendResult(interp, "bad option \"", argv[2],
779            "\": should be label", (char*)NULL);
780        return TCL_ERROR;
781    }
782    else if (c == 'd' && strcmp(argv[1],"data") == 0) {
783        if (argc < 3) {
784            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
785                argv[1], " option ?arg arg...?\"", (char*)NULL);
786            return TCL_ERROR;
787        }
788        c = *argv[2];
789        if (c == 's' && strcmp(argv[2],"state") == 0) {
790            if (argc < 4) {
791                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
792                    argv[1], " state on|off ?volume ...?\"", (char*)NULL);
793                return TCL_ERROR;
794            }
795
796            int state;
797            if (Tcl_GetBoolean(interp, argv[3], &state) != TCL_OK) {
798                return TCL_ERROR;
799            }
800
801            vector<int> ivol;
802            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
803                return TCL_ERROR;
804            }
805
806            vector<int>::iterator iter = ivol.begin();
807            while (iter != ivol.end()) {
808                if (state) {
809                    volume[*iter]->enable_data();
810                } else {
811                    volume[*iter]->disable_data();
812                }
813                ++iter;
814            }
815            return TCL_OK;
816        }
817        else if (c == 'f' && strcmp(argv[2],"follows") == 0) {
818            printf("Data Loading\n");
819            //fflush(stdout);
820            //return TCL_OK;
821
822            int nbytes;
823            if (Tcl_GetInt(interp, argv[3], &nbytes) != TCL_OK) {
824                return TCL_ERROR;
825            }
826
827            Rappture::Outcome err;
828            Rappture::Buffer buf;
829
830            // DEBUG
831            int totalsize = nbytes;
832            char buffer[8096];
833            while (nbytes > 0)
834            {
835                int chunk = (sizeof(buffer) < nbytes) ? sizeof(buffer) : nbytes;
836                int status = fread(buffer, 1, chunk, stdin);
837                //printf("Begin Reading [%d Read : %d Left]\n", status, nbytes - status);
838                fflush(stdout);
839                if (status > 0) {
840                    buf.append(buffer,status);
841                    nbytes -= status;
842                } else {
843                    printf("data unpacking failed\n");
844                    Tcl_AppendResult(interp, "data unpacking failed: unexpected EOF",
845                        (char*)NULL);
846                    return TCL_ERROR;
847                }
848            }
849
850            err = Rappture::encoding::decode(buf,RPENC_Z|RPENC_B64|RPENC_HDR);
851            if (err) {
852                printf("ERROR -- DECODING\n");
853                fflush(stdout);
854                Tcl_AppendResult(interp, err.remark().c_str(), (char*)NULL);
855                return TCL_ERROR;
856            }
857
858            int n = n_volumes;
859            char header[6];
860            memcpy(header, buf.bytes(), sizeof(char) * 5);
861            header[5] = '\0';
862
863#ifdef _LOCAL_ZINC_TEST_
864            //FILE* fp = fopen("/home/nanohub/vrinside/nv/data/HOON/QDWL_100_100_50_strain_8000i.nd_zatom_12_1", "rb");
865            FILE* fp = fopen("/home/nanohub/vrinside/nv/data/HOON/GaAs_AlGaAs_2QD_B4.nd_zc_1_wf", "rb");
866            unsigned char* b = (unsigned char*) malloc(buf.size());
867            if (fp == 0)
868            {
869                printf("cannot open the file\n");
870                fflush(stdout);
871                return TCL_ERROR;
872            }
873            fread(b, buf.size(), 1, fp);
874            fclose(fp);
875#endif
876
877           
878            printf("Checking header[%s]\n", header);
879            fflush(stdout);
880            if (!strcmp(header, "<HDR>"))
881            {
882                Volume* vol = NULL;
883
884                printf("ZincBlende stream is in\n");
885                fflush(stdout);
886                //std::stringstream fdata(std::ios_base::out|std::ios_base::in|std::ios_base::binary);
887                //fdata.write(buf.bytes(),buf.size());
888                //vol = NvZincBlendeReconstructor::getInstance()->loadFromStream(fdata);
889               
890#ifdef _LOCAL_ZINC_TEST_
891                vol = NvZincBlendeReconstructor::getInstance()->loadFromMemory(b);
892#else
893                vol = NvZincBlendeReconstructor::getInstance()->loadFromMemory((void*) buf.bytes());
894#endif
895
896                printf("finish loading\n");
897                fflush(stdout);
898                if (vol)
899                {
900                    while (n_volumes <= n)
901                    {
902                        volume.push_back((Volume*) NULL);
903                        n_volumes++;
904                    }
905
906                    if (volume[n] != NULL)
907                    {
908                        delete volume[n];
909                        volume[n] = NULL;
910                    }
911
912                    float dx0 = -0.5;
913                    float dy0 = -0.5*vol->height/vol->width;
914                    float dz0 = -0.5*vol->depth/vol->width;
915                    vol->move(Vector3(dx0, dy0, dz0));
916
917                    volume[n] = vol;
918                }
919            }
920/*
921            else if (!strcmp(header, "<FET>"))
922            {
923                printf("FET loading...\n");
924                fflush(stdout);
925                std::stringstream fdata;
926                fdata.write(buf.bytes(),buf.size());
927                //err = load_volume_stream3(n, fdata);
928
929                if (err) {
930                    Tcl_AppendResult(interp, err.remark().c_str(), (char*)NULL);
931                    return TCL_ERROR;
932                }
933            }
934*/
935            else
936            {
937                printf("OpenDX loading...\n");
938                fflush(stdout);
939                std::stringstream fdata;
940                fdata.write(buf.bytes(),buf.size());
941                err = load_volume_stream(n, fdata);
942                //err = load_volume_stream2(n, fdata);
943                if (err) {
944                    Tcl_AppendResult(interp, err.remark().c_str(), (char*)NULL);
945                    return TCL_ERROR;
946                }
947            }
948           
949
950            //
951            // BE CAREFUL:  Set the number of slices to something
952            //   slightly different for each volume.  If we have
953            //   identical volumes at exactly the same position
954            //   with exactly the same number of slices, the second
955            //   volume will overwrite the first, so the first won't
956            //   appear at all.
957            //
958            if (volume[n])
959            {
960                volume[n]->set_n_slice(256-n);
961                volume[n]->disable_cutplane(0);
962                volume[n]->disable_cutplane(1);
963                volume[n]->disable_cutplane(2);
964
965                g_vol_render->add_volume(volume[n], get_transfunc("default"));
966            }
967
968            return TCL_OK;
969        }
970        Tcl_AppendResult(interp, "bad option \"", argv[2],
971            "\": should be follows or state", (char*)NULL);
972        return TCL_ERROR;
973    }
974    else if (c == 'o' && strcmp(argv[1],"outline") == 0) {
975        if (argc < 3) {
976            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
977                argv[1], " option ?arg arg...?\"", (char*)NULL);
978            return TCL_ERROR;
979        }
980        c = *argv[2];
981        if (c == 's' && strcmp(argv[2],"state") == 0) {
982            if (argc < 3) {
983                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
984                    argv[1], " state on|off ?volume ...? \"", (char*)NULL);
985                return TCL_ERROR;
986            }
987
988            int state;
989            if (Tcl_GetBoolean(interp, argv[3], &state) != TCL_OK) {
990                return TCL_ERROR;
991            }
992
993            vector<int> ivol;
994            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
995                return TCL_ERROR;
996            }
997
998            vector<int>::iterator iter = ivol.begin();
999            while (iter != ivol.end()) {
1000                if (state) {
1001                    volume[*iter]->enable_outline();
1002                } else {
1003                    volume[*iter]->disable_outline();
1004                }
1005                ++iter;
1006            }
1007            return TCL_OK;
1008        }
1009        else if (c == 'v' && strcmp(argv[2],"visible") == 0) {
1010            if (argv[3] == "false")
1011            {
1012                for (int i = 0; i < n_volumes; ++i)
1013                {
1014                    if (volume[i]) volume[i]->disable_outline();
1015                }
1016            }
1017            else if (argv[3] == "true")
1018            {
1019                for (int i = 0; i < n_volumes; ++i)
1020                {
1021                    if (volume[i]) volume[i]->enable_outline();
1022                }
1023            }
1024
1025            return TCL_OK;
1026        }
1027        else if (c == 'c' && strcmp(argv[2],"color") == 0) {
1028            if (argc < 3) {
1029                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1030                    argv[1], " color {R G B} ?volume ...? \"", (char*)NULL);
1031                return TCL_ERROR;
1032            }
1033
1034            float rgb[3];
1035            if (GetColor(interp, (char*) argv[3], rgb) != TCL_OK) {
1036                return TCL_ERROR;
1037            }
1038
1039            vector<int> ivol;
1040            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
1041                return TCL_ERROR;
1042            }
1043
1044            vector<int>::iterator iter = ivol.begin();
1045            while (iter != ivol.end()) {
1046                volume[*iter]->set_outline_color(rgb);
1047                ++iter;
1048            }
1049            return TCL_OK;
1050        }
1051
1052        Tcl_AppendResult(interp, "bad option \"", argv[2],
1053            "\": should be color or state", (char*)NULL);
1054        return TCL_ERROR;
1055    }
1056    else if (c == 's' && strcmp(argv[1],"shading") == 0) {
1057        if (argc < 3) {
1058            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1059                argv[1], " option ?arg arg...?\"", (char*)NULL);
1060            return TCL_ERROR;
1061        }
1062        c = *argv[2];
1063        if (c == 't' && strcmp(argv[2],"transfunc") == 0) {
1064            if (argc < 4) {
1065                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1066                    argv[1], " transfunc name ?volume ...?\"", (char*)NULL);
1067                return TCL_ERROR;
1068            }
1069
1070            TransferFunction *tf = get_transfunc((char*)argv[3]);
1071            if (tf == NULL) {
1072                Tcl_AppendResult(interp, "transfer function \"", argv[3],
1073                    "\" is not defined", (char*)NULL);
1074                return TCL_ERROR;
1075            }
1076
1077            vector<int> ivol;
1078            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
1079                return TCL_ERROR;
1080            }
1081
1082            vector<int>::iterator iter = ivol.begin();
1083            while (iter != ivol.end()) {
1084                g_vol_render->shade_volume(volume[*iter], tf);
1085                ++iter;
1086            }
1087            return TCL_OK;
1088        }
1089        else if (c == 'd' && strcmp(argv[2],"diffuse") == 0) {
1090            if (argc < 4) {
1091                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1092                    argv[1], " diffuse value ?volume ...?\"", (char*)NULL);
1093                return TCL_ERROR;
1094            }
1095
1096            double dval;
1097            if (Tcl_GetDouble(interp, argv[3], &dval) != TCL_OK) {
1098                return TCL_ERROR;
1099            }
1100
1101            vector<int> ivol;
1102            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
1103                return TCL_ERROR;
1104            }
1105
1106            vector<int>::iterator iter = ivol.begin();
1107            while (iter != ivol.end()) {
1108                volume[*iter]->set_diffuse((float)dval);
1109                ++iter;
1110            }
1111            return TCL_OK;
1112        }
1113        else if (c == 'o' && strcmp(argv[2],"opacity") == 0) {
1114            if (argc < 4) {
1115                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1116                    argv[1], " opacity value ?volume ...?\"", (char*)NULL);
1117                return TCL_ERROR;
1118            }
1119
1120            double dval;
1121            if (Tcl_GetDouble(interp, argv[3], &dval) != TCL_OK) {
1122                return TCL_ERROR;
1123            }
1124
1125            vector<int> ivol;
1126            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
1127                return TCL_ERROR;
1128            }
1129
1130            vector<int>::iterator iter = ivol.begin();
1131            while (iter != ivol.end()) {
1132                volume[*iter]->set_opacity_scale((float)dval);
1133                ++iter;
1134            }
1135            return TCL_OK;
1136        }
1137        else if (c == 's' && strcmp(argv[2],"specular") == 0) {
1138            if (argc < 4) {
1139                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1140                    argv[1], " specular value ?volume ...?\"", (char*)NULL);
1141                return TCL_ERROR;
1142            }
1143
1144            double dval;
1145            if (Tcl_GetDouble(interp, argv[3], &dval) != TCL_OK) {
1146                return TCL_ERROR;
1147            }
1148
1149            vector<int> ivol;
1150            if (GetVolumeIndices(interp, argc-4, argv+4, &ivol) != TCL_OK) {
1151                return TCL_ERROR;
1152            }
1153
1154            vector<int>::iterator iter = ivol.begin();
1155            while (iter != ivol.end()) {
1156                volume[*iter]->set_specular((float)dval);
1157                ++iter;
1158            }
1159            return TCL_OK;
1160        }
1161        Tcl_AppendResult(interp, "bad option \"", argv[2],
1162            "\": should be diffuse, opacity, specular, or transfunc", (char*)NULL);
1163        return TCL_ERROR;
1164    }
1165    else if (c == 's' && strcmp(argv[1],"state") == 0) {
1166        if (argc < 3) {
1167            Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1168                argv[1], " on|off ?volume...?\"", (char*)NULL);
1169            return TCL_ERROR;
1170        }
1171
1172        int state;
1173        if (Tcl_GetBoolean(interp, argv[2], &state) != TCL_OK) {
1174            return TCL_ERROR;
1175        }
1176
1177        vector<int> ivol;
1178        if (GetVolumeIndices(interp, argc-3, argv+3, &ivol) != TCL_OK) {
1179            return TCL_ERROR;
1180        }
1181
1182        vector<int>::iterator iter = ivol.begin();
1183        while (iter != ivol.end()) {
1184            if (state) {
1185                volume[*iter]->enable();
1186            } else {
1187                volume[*iter]->disable();
1188            }
1189            ++iter;
1190        }
1191        return TCL_OK;
1192    }
1193    else if (c == 't' && strcmp(argv[1],"test2") == 0) {
1194        volume[1]->disable_data();
1195        volume[1]->disable();
1196        return TCL_OK;
1197    }
1198
1199    Tcl_AppendResult(interp, "bad option \"", argv[1],
1200        "\": should be data, outline, shading, or state", (char*)NULL);
1201    return TCL_ERROR;
1202}
1203
1204int HeightMapCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]))
1205{
1206    if (argc < 2)
1207    {
1208        {   
1209        srand( (unsigned)time( NULL ) );
1210        int size = 20 * 20;
1211        float sigma = 5.0f;
1212        float mean = exp(0.0f) / (sigma * sqrt(2.0f));
1213        float* data = (float*) malloc(sizeof(float) * size);
1214
1215        float x;
1216        for (int i = 0; i < size; ++i)
1217        {
1218            x = - 10 + i%20;
1219            data[i] = exp(- (x * x)/(2 * sigma * sigma)) / (sigma * sqrt(2.0f)) / mean;
1220        }
1221
1222        HeightMap* heightMap = new HeightMap();
1223        heightMap->setHeight(0, 0, 1, 1, 20, 20, data);
1224        heightMap->setColorMap(get_transfunc("default"));
1225        heightMap->setVisible(true);
1226        heightMap->setLineContourVisible(true);
1227        g_heightMap.push_back(heightMap);
1228        }
1229
1230        return TCL_OK;
1231    }
1232
1233    char c = *argv[1];
1234    if (c == 'd' && strcmp(argv[1],"data") == 0)
1235    {
1236        //bytes
1237        vector<int> indices;
1238        if (strcmp(argv[2],"visible") == 0)
1239        {
1240            bool visible = !strcmp(argv[3], "true");
1241           
1242            if (GetIndices(interp, argc-4, argv+4, &indices) != TCL_OK)
1243            {
1244               return TCL_ERROR;
1245            }
1246
1247            for (int i = 0; i < indices.size(); ++i)
1248            {
1249                if ((indices[i] < g_heightMap.size()) && (g_heightMap[indices[i]] != NULL))
1250                {
1251                    g_heightMap[indices[i]]->setVisible(visible);
1252                }
1253            }
1254            return TCL_OK;
1255        }
1256        else if (c == 'f' && strcmp(argv[2],"follows") == 0) {
1257            int nbytes;
1258            if (Tcl_GetInt(interp, argv[3], &nbytes) != TCL_OK) {
1259                return TCL_ERROR;
1260            }
1261        }
1262    }
1263    else if (c == 'l' && (strcmp(argv[1], "linecontour") == 0))
1264    {
1265        //bytes
1266        vector<int> indices;
1267        if (strcmp(argv[2],"visible") == 0)
1268        {
1269           
1270            bool visible = !(strcmp("true", argv[3]));
1271            printf("heightmap linecontour visible %s\n", (visible)?"true":"false");
1272            if (GetIndices(interp, argc-4, argv+4, &indices) != TCL_OK)
1273            {
1274                return TCL_ERROR;
1275            }
1276
1277            for (int i = 0; i < indices.size(); ++i)
1278            {
1279                printf("heightmap index %d\n");
1280                if ((indices[i] < g_heightMap.size()) && (g_heightMap[indices[i]] != NULL))
1281                {
1282                    printf("heightmap index %d visible applied\n");
1283                    g_heightMap[indices[i]]->setLineContourVisible(visible);
1284                }
1285            }
1286            return TCL_OK;
1287        }
1288        else if (strcmp(argv[2],"color") == 0)
1289        {
1290            double r, g, b;
1291            if ((Tcl_GetDouble(interp, argv[3], &r) == TCL_OK) &&
1292                (Tcl_GetDouble(interp, argv[4], &g) == TCL_OK) &&
1293                (Tcl_GetDouble(interp, argv[5], &b) == TCL_OK)) {
1294                r = r / 255.0;
1295                g = g / 255.0;
1296                b = b / 255.0;
1297            }
1298            else
1299            {
1300                return TCL_ERROR;
1301            }
1302
1303            vector<int> indices;
1304            if (GetIndices(interp, argc-6, argv+6, &indices) != TCL_OK)
1305            {
1306                return TCL_ERROR;
1307            }
1308            for (int i = 0; i < indices.size(); ++i)
1309            {
1310                if ((indices[i] < g_heightMap.size()) && (g_heightMap[indices[i]] != NULL))
1311                {
1312                    g_heightMap[indices[i]]->setLineContourColor(r, g, b);
1313                }
1314            }
1315
1316            return TCL_OK;
1317        }
1318    }
1319    else if (c == 't' && (strcmp(argv[1], "transfunc") == 0))
1320    {
1321        TransferFunction *tf = get_transfunc((char*)argv[2]);
1322        if (tf == NULL) {
1323            Tcl_AppendResult(interp, "transfer function \"", argv[3],
1324                "\" is not defined", (char*)NULL);
1325            return TCL_ERROR;
1326        }
1327
1328        vector<int> indices;
1329        if (GetVolumeIndices(interp, argc - 3, argv + 3, &indices) != TCL_OK)
1330        {
1331            for (int i = 0; i < indices.size(); ++i)
1332            {
1333                if ((indices[i] < g_heightMap.size()) && (g_heightMap[indices[i]] != NULL))
1334                {
1335                    g_heightMap[indices[i]]->setColorMap(tf);
1336                }
1337            }
1338        }
1339        return TCL_OK;
1340    }
1341   
1342    Tcl_AppendResult(interp, "bad option \"", argv[1],
1343        "\": should be data, outline, shading, or state", (char*)NULL);
1344    return TCL_ERROR;
1345}
1346
1347int GridCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]))
1348{
1349    char c = *argv[1];
1350    if (c == 'v' && strcmp(argv[1],"visible") == 0)
1351    {
1352        if (strcmp(argv[2],"true") == 0)
1353        {
1354            g_grid->setVisible(true);
1355            return TCL_OK;
1356        }
1357        else if (strcmp(argv[2],"false") == 0)
1358        {
1359            g_grid->setVisible(false);
1360            return TCL_OK;
1361        }
1362    }
1363    else if (c == 'l' && strcmp(argv[1],"linecount") == 0)
1364    {
1365        int x, y, z;
1366
1367        if ((Tcl_GetInt(interp, argv[2], &x) == TCL_OK) &&
1368            (Tcl_GetInt(interp, argv[3], &y) == TCL_OK) &&
1369            (Tcl_GetInt(interp, argv[4], &z) == TCL_OK)) {
1370
1371            if (g_grid) g_grid->setGridLineCount(x, y, z);
1372
1373            return TCL_OK;
1374        }
1375    }
1376    else if (c == 'a' && strcmp(argv[1],"axiscolor") == 0)
1377    {
1378        int r, g, b;
1379        if ((Tcl_GetInt(interp, argv[2], &r) == TCL_OK) &&
1380            (Tcl_GetInt(interp, argv[3], &g) == TCL_OK) &&
1381            (Tcl_GetInt(interp, argv[4], &b) == TCL_OK)) {
1382
1383            if (g_grid) g_grid->setAxisColor(r / 255.0f, g / 255.0f, b / 255.0f);
1384            return TCL_OK;
1385        }
1386    }
1387    else if (c == 'l' && strcmp(argv[1],"linecolor") == 0)
1388    {
1389        int r, g, b;
1390        if ((Tcl_GetInt(interp, argv[2], &r) == TCL_OK) &&
1391            (Tcl_GetInt(interp, argv[3], &g) == TCL_OK) &&
1392            (Tcl_GetInt(interp, argv[4], &b) == TCL_OK)) {
1393
1394            if (g_grid) g_grid->setGridLineColor(r / 255.0f, g / 255.0f, b / 255.0f);
1395            return TCL_OK;
1396        }
1397    }
1398    else if (c == 'm')
1399    {
1400        if (strcmp(argv[1],"minmax") == 0)
1401        {
1402            double x1, y1, z1, x2, y2, z2;
1403            if ((Tcl_GetDouble(interp, argv[2], &x1) == TCL_OK) &&
1404                (Tcl_GetDouble(interp, argv[3], &y1) == TCL_OK) &&
1405                (Tcl_GetDouble(interp, argv[4], &z1) == TCL_OK) &&
1406                (Tcl_GetDouble(interp, argv[5], &x2) == TCL_OK) &&
1407                (Tcl_GetDouble(interp, argv[6], &y2) == TCL_OK) &&
1408                (Tcl_GetDouble(interp, argv[7], &z2) == TCL_OK)) {
1409
1410                if (g_grid) g_grid->setMinMax(Vector3(x1, y1, z1), Vector3(x2, y2, z2));
1411
1412                return TCL_OK;
1413            }
1414        }
1415    }
1416    else if (c == 'a' && strcmp(argv[1],"axisname") == 0)
1417    {
1418        int axisID = 0;
1419        if (!strcmp(argv[2], "x")) axisID = 0;
1420        if (!strcmp(argv[2], "y")) axisID = 1;
1421        if (!strcmp(argv[2], "z")) axisID = 2;
1422       
1423        if (g_grid) g_grid->setAxisName(axisID, argv[3]);
1424        return TCL_OK;
1425    }
1426
1427    Tcl_AppendResult(interp, "bad option \"", argv[1],
1428        "\": should be data, outline, shading, or state", (char*)NULL);
1429    return TCL_ERROR;
1430}
1431
1432int AxisCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]))
1433{
1434    if (argc < 2) {
1435        Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1436            " option arg arg...\"", (char*)NULL);
1437        return TCL_ERROR;
1438    }
1439
1440    char c = *argv[1];
1441    if (c == 'v' && strcmp(argv[1],"visible") == 0)
1442    {
1443        if (strcmp(argv[2],"true") == 0)
1444        {
1445            axis_on = true;
1446        }
1447        else if (strcmp(argv[2],"false") == 0)
1448        {
1449            axis_on = false;
1450        }
1451           
1452        return TCL_OK;
1453    }
1454
1455    Tcl_AppendResult(interp, "bad option \"", argv[1],
1456        "\": should be data, outline, shading, or state", (char*)NULL);
1457    return TCL_ERROR;
1458}
1459
1460int TestCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[]))
1461{
1462    std::ifstream fin;
1463    fin.open("/home/nanohub/vrinside/nv/data/fet/graph-finfet.txt", std::ios::binary);
1464                //std::stringstream fdata(std::ios_base::out|std::ios_base::in|std::ios_base::binary);
1465
1466    if (!fin)
1467    {
1468        printf("file not found\n");
1469        return TCL_OK;
1470    }
1471
1472    float temp1, temp2, temp3;
1473    float vmin, vmax;
1474    int count;
1475    char* start;
1476
1477    char line[128];
1478    char temp[24];
1479    float dx, dy, dz;
1480    float x0, y0, z0;
1481    do {
1482        fin.getline(line,sizeof(line)-1);
1483
1484        // skip leading blanks
1485        for (start=&line[0]; *start == ' ' || *start == '\t'; start++);
1486
1487        if (*start != '#')
1488        { 
1489            sscanf(start, "%s%f%f%f", temp, &temp1, &temp2, &temp3);
1490            if (!strcmp(temp, "size"))
1491            {
1492                dx = temp1;
1493                dy = temp2;
1494                dz = temp3;
1495            }
1496            else if (!strcmp(temp, "size"))
1497            {
1498                x0 = temp1;
1499                y0 = temp2;
1500                z0 = temp3;
1501            }
1502            else if (!strcmp(temp, "minmax"))
1503            {
1504                vmin = temp1;
1505                vmax = temp2;
1506            }
1507            else if (!strcmp(temp, "count"))
1508            {
1509                count = temp1;
1510            }
1511            else if (!strcmp(temp, "</FET>"))
1512            {
1513                break;
1514            }
1515        }
1516    } while (!fin.eof());
1517
1518    printf("data count %d \n", count);
1519    fflush(stdout);
1520
1521    Vector4* data = new Vector4[count];
1522    memset(data, 0, sizeof(Vector4) * count);
1523    float min = 1e21, max = 1e-21;
1524    for (int i = 0; i < count; ++i)
1525    {
1526        fin.getline(line,sizeof(line)-1);
1527        sscanf(line, "%f%f%f%f", &data[i].x, &data[i].y, &data[i].z, &data[i].w);
1528
1529        if (data[i].w < min) min = data[i].w;
1530        if (data[i].w > max) max = data[i].w;
1531    }
1532
1533    float d = max - min;
1534    for (int i = 0; i < count; ++i)
1535    {
1536        data[i].w = (data[i].w - min) / d;
1537    }
1538
1539    PointSet* newPointSet = new PointSet();
1540    newPointSet->initialize(data, count);
1541    newPointSet->setVisible(true);
1542    g_pointSet.push_back(newPointSet);
1543   
1544    return TCL_OK;
1545}
1546
1547/*
1548 * ----------------------------------------------------------------------
1549 * FUNCTION: GetVolumeIndices()
1550 *
1551 * Used internally to decode a series of volume index values and
1552 * store then in the specified vector.  If there are no volume index
1553 * arguments, this means "all volumes" to most commands, so all
1554 * active volume indices are stored in the vector.
1555 *
1556 * Updates pushes index values into the vector.  Returns TCL_OK or
1557 * TCL_ERROR to indicate an error.
1558 * ----------------------------------------------------------------------
1559 */
1560static int
1561GetVolumeIndices(Tcl_Interp *interp, int argc, CONST84 char *argv[],
1562    vector<int>* vectorPtr)
1563{
1564    if (argc == 0) {
1565        for (int n=0; n < volume.size(); n++) {
1566            if (volume[n] != NULL) {
1567                vectorPtr->push_back(n);
1568            }
1569        }
1570    } else {
1571        int ivol;
1572        for (int n=0; n < argc; n++) {
1573            if (Tcl_GetInt(interp, argv[n], &ivol) != TCL_OK) {
1574                return TCL_ERROR;
1575            }
1576            if (ivol < 0 || ivol >= volume.size()) {
1577                Tcl_AppendResult(interp, "bad volume index \"", argv[n],
1578                    "\"", (char*)NULL);
1579                return TCL_ERROR;
1580            }
1581            if (volume[ivol] != NULL) {
1582                vectorPtr->push_back(ivol);
1583            }
1584        }
1585    }
1586    return TCL_OK;
1587}
1588
1589static int
1590GetIndices(Tcl_Interp *interp, int argc, CONST84 char *argv[],
1591    vector<int>* vectorPtr)
1592{
1593    int ivol;
1594    for (int n=0; n < argc; n++)
1595    {
1596        if (Tcl_GetInt(interp, argv[n], &ivol) != TCL_OK) {
1597            return TCL_ERROR;
1598        }
1599        vectorPtr->push_back(ivol);
1600    }
1601    return TCL_OK;
1602}
1603
1604/*
1605 * ----------------------------------------------------------------------
1606 * FUNCTION: GetAxis()
1607 *
1608 * Used internally to decode an axis value from a string ("x", "y",
1609 * or "z") to its index (0, 1, or 2).  Returns TCL_OK if successful,
1610 * along with a value in valPtr.  Otherwise, it returns TCL_ERROR
1611 * and an error message in the interpreter.
1612 * ----------------------------------------------------------------------
1613 */
1614static int
1615GetAxis(Tcl_Interp *interp, char *str, int *valPtr)
1616{
1617    if (strcmp(str,"x") == 0) {
1618        *valPtr = 0;
1619        return TCL_OK;
1620    }
1621    else if (strcmp(str,"y") == 0) {
1622        *valPtr = 1;
1623        return TCL_OK;
1624    }
1625    else if (strcmp(str,"z") == 0) {
1626        *valPtr = 2;
1627        return TCL_OK;
1628    }
1629    Tcl_AppendResult(interp, "bad axis \"", str,
1630        "\": should be x, y, or z", (char*)NULL);
1631    return TCL_ERROR;
1632}
1633
1634/*
1635 * ----------------------------------------------------------------------
1636 * FUNCTION: GetColor()
1637 *
1638 * Used internally to decode a color value from a string ("R G B")
1639 * as a list of three numbers 0-1.  Returns TCL_OK if successful,
1640 * along with RGB values in valPtr.  Otherwise, it returns TCL_ERROR
1641 * and an error message in the interpreter.
1642 * ----------------------------------------------------------------------
1643 */
1644static int
1645GetColor(Tcl_Interp *interp, char *str, float *rgbPtr)
1646{
1647    int rgbc;
1648    char **rgbv;
1649    if (Tcl_SplitList(interp, str, &rgbc, (const char***)&rgbv) != TCL_OK) {
1650        return TCL_ERROR;
1651    }
1652    if (rgbc != 3) {
1653        Tcl_AppendResult(interp, "bad color \"", str,
1654            "\": should be {R G B} as double values 0-1", (char*)NULL);
1655        return TCL_ERROR;
1656    }
1657
1658    double rval, gval, bval;
1659    if (Tcl_GetDouble(interp, rgbv[0], &rval) != TCL_OK) {
1660        Tcl_Free((char*)rgbv);
1661        return TCL_ERROR;
1662    }
1663    if (Tcl_GetDouble(interp, rgbv[1], &gval) != TCL_OK) {
1664        Tcl_Free((char*)rgbv);
1665        return TCL_ERROR;
1666    }
1667    if (Tcl_GetDouble(interp, rgbv[2], &bval) != TCL_OK) {
1668        Tcl_Free((char*)rgbv);
1669        return TCL_ERROR;
1670    }
1671    Tcl_Free((char*)rgbv);
1672
1673    rgbPtr[0] = (float)rval;
1674    rgbPtr[1] = (float)gval;
1675    rgbPtr[2] = (float)bval;
1676
1677    return TCL_OK;
1678}
1679
1680/*
1681 * ----------------------------------------------------------------------
1682 * USAGE: debug("string", ...)
1683 *
1684 * Use this anywhere within the package to send a debug message
1685 * back to the client.  The string can have % fields as used by
1686 * the printf() package.  Any remaining arguments are treated as
1687 * field substitutions on that.
1688 * ----------------------------------------------------------------------
1689 */
1690void
1691debug(char *str)
1692{
1693    write(0, str, strlen(str));
1694}
1695
1696void
1697debug(char *str, double v1)
1698{
1699    char buffer[512];
1700    sprintf(buffer, str, v1);
1701    write(0, buffer, strlen(buffer));
1702}
1703
1704void
1705debug(char *str, double v1, double v2)
1706{
1707    char buffer[512];
1708    sprintf(buffer, str, v1, v2);
1709    write(0, buffer, strlen(buffer));
1710}
1711
1712void
1713debug(char *str, double v1, double v2, double v3)
1714{
1715    char buffer[512];
1716    sprintf(buffer, str, v1, v2, v3);
1717    write(0, buffer, strlen(buffer));
1718}
1719
1720
1721static int
1722PlaneNewCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])){
1723  fprintf(stderr, "load plane for 2D visualization command\n");
1724
1725  int index, w, h;
1726
1727  if (argc != 4) {
1728    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1729                " plane_index w h \"", (char*)NULL);
1730    return TCL_ERROR;
1731  }
1732  if (Tcl_GetInt(interp, argv[1], &index) != TCL_OK) {
1733        return TCL_ERROR;
1734  }
1735  if (Tcl_GetInt(interp, argv[2], &w) != TCL_OK) {
1736        return TCL_ERROR;
1737  }
1738  if (Tcl_GetInt(interp, argv[3], &h) != TCL_OK) {
1739        return TCL_ERROR;
1740  }
1741
1742  //Now read w*h*4 bytes. The server expects the plane to be a stream of floats
1743  char* tmp = new char[int(w*h*sizeof(float))];
1744  bzero(tmp, w*h*4);
1745  int status = read(0, tmp, w*h*sizeof(float));
1746  if (status <= 0){
1747    exit(0);
1748  }
1749 
1750  plane[index] = new Texture2D(w, h, GL_FLOAT, GL_LINEAR, 1, (float*)tmp);
1751 
1752  delete[] tmp;
1753  return TCL_OK;
1754}
1755
1756
1757static
1758int PlaneLinkCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])){
1759  fprintf(stderr, "link the plane to the 2D renderer command\n");
1760
1761  int plane_index, tf_index;
1762
1763  if (argc != 3) {
1764    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1765                " plane_index tf_index \"", (char*)NULL);
1766    return TCL_ERROR;
1767  }
1768  if (Tcl_GetInt(interp, argv[1], &plane_index) != TCL_OK) {
1769        return TCL_ERROR;
1770  }
1771  if (Tcl_GetInt(interp, argv[2], &tf_index) != TCL_OK) {
1772        return TCL_ERROR;
1773  }
1774
1775  //plane_render->add_plane(plane[plane_index], tf[tf_index]);
1776
1777  return TCL_OK;
1778}
1779
1780
1781//Enable a 2D plane for render
1782//The plane_index is the index mantained in the 2D plane renderer
1783static
1784int PlaneEnableCmd _ANSI_ARGS_((ClientData cdata, Tcl_Interp *interp, int argc, CONST84 char *argv[])){
1785  fprintf(stderr, "enable a plane so the 2D renderer can render it command\n");
1786
1787  int plane_index, mode;
1788
1789  if (argc != 3) {
1790    Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
1791                " plane_index mode \"", (char*)NULL);
1792    return TCL_ERROR;
1793  }
1794  if (Tcl_GetInt(interp, argv[1], &plane_index) != TCL_OK) {
1795        return TCL_ERROR;
1796  }
1797  if (Tcl_GetInt(interp, argv[2], &mode) != TCL_OK) {
1798        return TCL_ERROR;
1799  }
1800
1801  if(mode==0)
1802    plane_render->set_active_plane(-1);
1803  else
1804    plane_render->set_active_plane(plane_index);
1805
1806  return TCL_OK;
1807}
1808
1809//report errors related to CG shaders
1810void cgErrorCallback(void)
1811{
1812    CGerror lastError = cgGetError();
1813    if(lastError) {
1814        const char *listing = cgGetLastListing(g_context);
1815        printf("\n---------------------------------------------------\n");
1816        printf("%s\n\n", cgGetErrorString(lastError));
1817        printf("%s\n", listing);
1818        printf("-----------------------------------------------------\n");
1819        printf("Cg error, exiting...\n");
1820        cgDestroyContext(g_context);
1821        exit(-1);
1822    }
1823}
1824
1825/* Load a 3D vector field from a dx-format file
1826 */
1827void
1828load_vector_stream(int index, std::iostream& fin) {
1829    int dummy, nx, ny, nz, nxy, npts;
1830    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
1831    char line[128], type[128], *start;
1832
1833    do {
1834        fin.getline(line,sizeof(line)-1);
1835        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
1836            ;  // skip leading blanks
1837
1838        if (*start != '#') {  // skip comment lines
1839            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
1840                // found grid size
1841            }
1842            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
1843                // found origin
1844            }
1845            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
1846                // found one of the delta lines
1847                if (ddx != 0.0) { dx = ddx; }
1848                else if (ddy != 0.0) { dy = ddy; }
1849                else if (ddz != 0.0) { dz = ddz; }
1850            }
1851            else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
1852                if (npts != nx*ny*nz) {
1853                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
1854                    return;
1855                }
1856                break;
1857            }
1858            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
1859                if (npts != nx*ny*nz) {
1860                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
1861                    return;
1862                }
1863                break;
1864            }
1865        }
1866    } while (!fin.eof());
1867
1868    // read data points
1869    if (!fin.eof()) {
1870        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
1871        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
1872        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1873        Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
1874        Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
1875        Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
1876
1877        double vx, vy, vz;
1878        int nread = 0;
1879        for (int ix=0; ix < nx; ix++) {
1880            for (int iy=0; iy < ny; iy++) {
1881                for (int iz=0; iz < nz; iz++) {
1882                    if (fin.eof() || nread > npts) {
1883                        break;
1884                    }
1885                    fin.getline(line,sizeof(line)-1);
1886                    if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
1887                        int nindex = iz*nx*ny + iy*nx + ix;
1888                        xfield.define(nindex, vx);
1889                        yfield.define(nindex, vy);
1890                        zfield.define(nindex, vz);
1891                        nread++;
1892                    }
1893                }
1894            }
1895        }
1896
1897        // make sure that we read all of the expected points
1898        if (nread != nx*ny*nz) {
1899            std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << nread << " points" << std::endl;
1900            return;
1901        }
1902
1903        // figure out a good mesh spacing
1904        int nsample = 30;
1905        dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
1906        dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
1907        dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
1908        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1909
1910        nx = (int)ceil(dx/dmin);
1911        ny = (int)ceil(dy/dmin);
1912        nz = (int)ceil(dz/dmin);
1913
1914#ifndef NV40
1915        // must be an even power of 2 for older cards
1916            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1917            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1918            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1919#endif
1920
1921        float *data = new float[3*nx*ny*nz];
1922
1923        std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl;
1924
1925        // generate the uniformly sampled data that we need for a volume
1926        double vmin = 0.0;
1927        double vmax = 0.0;
1928        double nzero_min = 0.0;
1929        int ngen = 0;
1930        for (int iz=0; iz < nz; iz++) {
1931            double zval = z0 + iz*dmin;
1932            for (int iy=0; iy < ny; iy++) {
1933                double yval = y0 + iy*dmin;
1934                for (int ix=0; ix < nx; ix++) {
1935                    double xval = x0 + ix*dmin;
1936
1937                    vx = xfield.value(xval,yval,zval);
1938                    vy = yfield.value(xval,yval,zval);
1939                    vz = zfield.value(xval,yval,zval);
1940                    //vx = 1;
1941                    //vy = 1;
1942                    vz = 0;
1943                    double vm = sqrt(vx*vx + vy*vy + vz*vz);
1944
1945                    if (vm < vmin) { vmin = vm; }
1946                    if (vm > vmax) { vmax = vm; }
1947                    if (vm != 0.0f && vm < nzero_min)
1948                    {
1949                        nzero_min = vm;
1950                    }
1951
1952                    data[ngen++] = vx;
1953                    data[ngen++] = vy;
1954                    data[ngen++] = vz;
1955                }
1956            }
1957        }
1958
1959        ngen = 0;
1960        for (ngen=0; ngen < npts; ngen++) {
1961            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
1962        }
1963
1964        load_volume(index, nx, ny, nz, 3, data, vmin, vmax, nzero_min);
1965        delete [] data;
1966    } else {
1967        std::cerr << "WARNING: data not found in stream" << std::endl;
1968    }
1969}
1970
1971
1972/* Load a 3D volume from a dx-format file
1973 */
1974Rappture::Outcome
1975load_volume_stream2(int index, std::iostream& fin) {
1976    Rappture::Outcome result;
1977
1978    Rappture::MeshTri2D xymesh;
1979    int dummy, nx, ny, nz, nxy, npts;
1980    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
1981    char line[128], type[128], *start;
1982
1983    int isrect = 1;
1984
1985    float* voldata = 0;
1986    do {
1987        fin.getline(line,sizeof(line)-1);
1988        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
1989            ;  // skip leading blanks
1990
1991        if (*start != '#') {  // skip comment lines
1992            printf("%s\n", line);
1993            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
1994                // found grid size
1995                isrect = 1;
1996            }
1997            else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
1998                isrect = 0;
1999                double xx, yy, zz;
2000                for (int i=0; i < nxy; i++) {
2001                    fin.getline(line, sizeof(line));
2002                    fin.getline(line,sizeof(line)-1);
2003                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
2004                        xymesh.addNode( Rappture::Node2D(xx,yy) );
2005                    }
2006                }
2007                char mesg[256];
2008                sprintf(mesg,"test");
2009                result.error(mesg);
2010                return result;
2011
2012                char fpts[128];
2013                sprintf(fpts, "/tmp/tmppts%d", getpid());
2014                char fcells[128];
2015                sprintf(fcells, "/tmp/tmpcells%d", getpid());
2016
2017                std::ofstream ftmp(fpts);
2018                // save corners of bounding box first, to work around meshing
2019                // problems in voronoi utility
2020                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
2021                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
2022                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
2023                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
2024                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
2025                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
2026                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
2027                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
2028                for (int i=0; i < nxy; i++) {
2029                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
2030               
2031                }
2032                ftmp.close();
2033
2034                char cmdstr[512];
2035                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
2036                if (system(cmdstr) == 0) {
2037                    int cx, cy, cz;
2038                    std::ifstream ftri(fcells);
2039                    while (!ftri.eof()) {
2040                        ftri.getline(line,sizeof(line)-1);
2041                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
2042                            if (cx >= 4 && cy >= 4 && cz >= 4) {
2043                                // skip first 4 boundary points
2044                                xymesh.addCell(cx-4, cy-4, cz-4);
2045                            }
2046                        }
2047                    }
2048                    ftri.close();
2049                } else {
2050                    return result.error("triangularization failed");
2051                }
2052
2053                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
2054                system(cmdstr);
2055            }
2056            else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
2057                // found z-grid
2058            }
2059            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
2060                // found origin
2061            }
2062            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
2063                // found one of the delta lines
2064                if (ddx != 0.0) { dx = ddx; }
2065                else if (ddy != 0.0) { dy = ddy; }
2066                else if (ddz != 0.0) { dz = ddz; }
2067            }
2068            else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
2069                if (isrect && (npts != nx*ny*nz)) {
2070                    char mesg[256];
2071                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
2072                    return result.error(mesg);
2073                }
2074                else if (!isrect && (npts != nxy*nz)) {
2075                    char mesg[256];
2076                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
2077                    return result.error(mesg);
2078                }
2079                break;
2080            }
2081            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
2082                if (npts != nx*ny*nz) {
2083                    char mesg[256];
2084                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
2085                    return result.error(mesg);
2086                }
2087                break;
2088            }
2089        }
2090    } while (!fin.eof());
2091
2092    // read data points
2093    if (!fin.eof()) {
2094        if (isrect) {
2095            double dval[6];
2096            int nread = 0;
2097            int ix = 0;
2098            int iy = 0;
2099            int iz = 0;
2100            float* data = new float[nx *  ny *  nz * 4];
2101            memset(data, 0, nx*ny*nz*4);
2102            double vmin = 1e21;
2103            double nzero_min = 1e21;
2104            double vmax = -1e21;
2105
2106
2107            while (!fin.eof() && nread < npts) {
2108                fin.getline(line,sizeof(line)-1);
2109                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
2110
2111                for (int p=0; p < n; p++) {
2112                    int nindex = (iz*nx*ny + iy*nx + ix) * 4;
2113                    data[nindex] = dval[p];
2114
2115                    if (dval[p] < vmin) vmin = dval[p];
2116                    if (dval[p] > vmax) vmax = dval[p];
2117                    if (dval[p] != 0.0f && dval[p] < nzero_min)
2118                    {
2119                         nzero_min = dval[p];
2120                    }
2121
2122                    nread++;
2123                    if (++iz >= nz) {
2124                        iz = 0;
2125                        if (++iy >= ny) {
2126                            iy = 0;
2127                            ++ix;
2128                        }
2129                    }
2130                }
2131            }
2132
2133            // make sure that we read all of the expected points
2134            if (nread != nx*ny*nz) {
2135                char mesg[256];
2136                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
2137                result.error(mesg);
2138                return result;
2139            }
2140
2141            double dv = vmax - vmin;
2142            int count = nx*ny*nz;
2143            int ngen = 0;
2144            double v;
2145            printf("test2\n");
2146                        fflush(stdout);
2147            if (dv == 0.0) { dv = 1.0; }
2148            for (int i = 0; i < count; ++i)
2149            {
2150                v = data[ngen];
2151                // scale all values [0-1], -1 => out of bounds
2152                v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2153                data[ngen] = v;
2154                ngen += 4;
2155            }
2156
2157            // Compute the gradient of this data.  BE CAREFUL: center
2158            // calculation on each node to avoid skew in either direction.
2159            ngen = 0;
2160            for (int iz=0; iz < nz; iz++) {
2161                for (int iy=0; iy < ny; iy++) {
2162                    for (int ix=0; ix < nx; ix++) {
2163                        // gradient in x-direction
2164                        double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
2165                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
2166                        if (valm1 < 0 || valp1 < 0) {
2167                            data[ngen+1] = 0.0;
2168                        } else {
2169                            data[ngen+1] = valp1-valm1; // assume dx=1
2170                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2171                        }
2172
2173                        // gradient in y-direction
2174                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
2175                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
2176                        if (valm1 < 0 || valp1 < 0) {
2177                            data[ngen+2] = 0.0;
2178                        } else {
2179                            data[ngen+2] = valp1-valm1; // assume dx=1
2180                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
2181                        }
2182
2183                        // gradient in z-direction
2184                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
2185                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
2186                        if (valm1 < 0 || valp1 < 0) {
2187                            data[ngen+3] = 0.0;
2188                        } else {
2189                            data[ngen+3] = valp1-valm1; // assume dx=1
2190                            //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
2191                        }
2192
2193                        ngen += 4;
2194                    }
2195                }
2196            }
2197
2198            dx = nx;
2199            dy = ny;
2200            dz = nz;
2201            load_volume(index, nx, ny, nz, 4, data,
2202                vmin, vmax, nzero_min);
2203
2204            delete [] data;
2205
2206        } else {
2207            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
2208            Rappture::FieldPrism3D field(xymesh, zgrid);
2209
2210            double dval;
2211            int nread = 0;
2212            int ixy = 0;
2213            int iz = 0;
2214            while (!fin.eof() && nread < npts) {
2215                if (!(fin >> dval).fail()) {
2216                    int nid = nxy*iz + ixy;
2217                    field.define(nid, dval);
2218
2219                    nread++;
2220                    if (++iz >= nz) {
2221                        iz = 0;
2222                        ixy++;
2223                    }
2224                }
2225            }
2226
2227            // make sure that we read all of the expected points
2228            if (nread != nxy*nz) {
2229                char mesg[256];
2230                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
2231                return result.error(mesg);
2232            }
2233
2234            // figure out a good mesh spacing
2235            int nsample = 30;
2236            x0 = field.rangeMin(Rappture::xaxis);
2237            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
2238            y0 = field.rangeMin(Rappture::yaxis);
2239            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
2240            z0 = field.rangeMin(Rappture::zaxis);
2241            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
2242            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
2243
2244            nx = (int)ceil(dx/dmin);
2245            ny = (int)ceil(dy/dmin);
2246            nz = (int)ceil(dz/dmin);
2247#ifndef NV40
2248            // must be an even power of 2 for older cards
2249                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
2250                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
2251                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
2252#endif
2253            float *data = new float[4*nx*ny*nz];
2254
2255            double vmin = field.valueMin();
2256            double dv = field.valueMax() - field.valueMin();
2257            if (dv == 0.0) { dv = 1.0; }
2258
2259            // generate the uniformly sampled data that we need for a volume
2260            int ngen = 0;
2261            double nzero_min = 0.0;
2262            for (iz=0; iz < nz; iz++) {
2263                double zval = z0 + iz*dmin;
2264                for (int iy=0; iy < ny; iy++) {
2265                    double yval = y0 + iy*dmin;
2266                    for (int ix=0; ix < nx; ix++) {
2267                        double xval = x0 + ix*dmin;
2268                        double v = field.value(xval,yval,zval);
2269
2270                        if (v != 0.0f && v < nzero_min)
2271                        {
2272                            nzero_min = v;
2273                        }
2274                        // scale all values [0-1], -1 => out of bounds
2275                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2276                        data[ngen] = v;
2277
2278                        ngen += 4;
2279                    }
2280                }
2281            }
2282
2283            // Compute the gradient of this data.  BE CAREFUL: center
2284            // calculation on each node to avoid skew in either direction.
2285            ngen = 0;
2286            for (int iz=0; iz < nz; iz++) {
2287                for (int iy=0; iy < ny; iy++) {
2288                    for (int ix=0; ix < nx; ix++) {
2289                        // gradient in x-direction
2290                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
2291                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
2292                        if (valm1 < 0 || valp1 < 0) {
2293                            data[ngen+1] = 0.0;
2294                        } else {
2295                            //data[ngen+1] = valp1-valm1; // assume dx=1
2296                            data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
2297                        }
2298
2299                        // gradient in y-direction
2300                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
2301                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
2302                        if (valm1 < 0 || valp1 < 0) {
2303                            data[ngen+2] = 0.0;
2304                        } else {
2305                            //data[ngen+2] = valp1-valm1; // assume dy=1
2306                            data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
2307                        }
2308
2309                        // gradient in z-direction
2310                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
2311                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
2312                        if (valm1 < 0 || valp1 < 0) {
2313                            data[ngen+3] = 0.0;
2314                        } else {
2315                            //data[ngen+3] = valp1-valm1; // assume dz=1
2316                            data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
2317                        }
2318
2319                        ngen += 4;
2320                    }
2321                }
2322            }
2323
2324            load_volume(index, nx, ny, nz, 4, data,
2325                field.valueMin(), field.valueMax(), nzero_min);
2326
2327            delete [] data;
2328        }
2329    } else {
2330        return result.error("data not found in stream");
2331    }
2332
2333    //
2334    // Center this new volume on the origin.
2335    //
2336    float dx0 = -0.5;
2337    float dy0 = -0.5*dy/dx;
2338    float dz0 = -0.5*dz/dx;
2339    volume[index]->move(Vector3(dx0, dy0, dz0));
2340
2341    return result;
2342}
2343
2344Rappture::Outcome
2345load_volume_stream(int index, std::iostream& fin) {
2346    Rappture::Outcome result;
2347
2348    Rappture::MeshTri2D xymesh;
2349    int dummy, nx, ny, nz, nxy, npts;
2350    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
2351    char line[128], type[128], *start;
2352
2353    int isrect = 1;
2354
2355    do {
2356        fin.getline(line,sizeof(line)-1);
2357        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
2358            ;  // skip leading blanks
2359
2360        if (*start != '#') {  // skip comment lines
2361            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
2362                // found grid size
2363                isrect = 1;
2364            }
2365            else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
2366                isrect = 0;
2367                double xx, yy, zz;
2368                for (int i=0; i < nxy; i++) {
2369                    fin.getline(line,sizeof(line)-1);
2370                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
2371                        xymesh.addNode( Rappture::Node2D(xx,yy) );
2372                    }
2373                }
2374
2375                char fpts[128];
2376                sprintf(fpts, "/tmp/tmppts%d", getpid());
2377                char fcells[128];
2378                sprintf(fcells, "/tmp/tmpcells%d", getpid());
2379
2380                std::ofstream ftmp(fpts);
2381                // save corners of bounding box first, to work around meshing
2382                // problems in voronoi utility
2383                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
2384                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
2385                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
2386                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
2387                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
2388                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
2389                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
2390                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
2391                for (int i=0; i < nxy; i++) {
2392                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
2393               
2394                }
2395                ftmp.close();
2396
2397                char cmdstr[512];
2398                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
2399                if (system(cmdstr) == 0) {
2400                    int cx, cy, cz;
2401                    std::ifstream ftri(fcells);
2402                    while (!ftri.eof()) {
2403                        ftri.getline(line,sizeof(line)-1);
2404                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
2405                            if (cx >= 4 && cy >= 4 && cz >= 4) {
2406                                // skip first 4 boundary points
2407                                xymesh.addCell(cx-4, cy-4, cz-4);
2408                            }
2409                        }
2410                    }
2411                    ftri.close();
2412                } else {
2413                    return result.error("triangularization failed");
2414                }
2415
2416                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
2417                system(cmdstr);
2418            }
2419            else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
2420                // found z-grid
2421            }
2422            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
2423                // found origin
2424            }
2425            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
2426                // found one of the delta lines
2427                if (ddx != 0.0) { dx = ddx; }
2428                else if (ddy != 0.0) { dy = ddy; }
2429                else if (ddz != 0.0) { dz = ddz; }
2430            }
2431            else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
2432                if (isrect && (npts != nx*ny*nz)) {
2433                    char mesg[256];
2434                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
2435                    return result.error(mesg);
2436                }
2437                else if (!isrect && (npts != nxy*nz)) {
2438                    char mesg[256];
2439                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
2440                    return result.error(mesg);
2441                }
2442                break;
2443            }
2444            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
2445                if (npts != nx*ny*nz) {
2446                    char mesg[256];
2447                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
2448                    return result.error(mesg);
2449                }
2450                break;
2451            }
2452        }
2453    } while (!fin.eof());
2454
2455    // read data points
2456    if (!fin.eof()) {
2457        if (isrect) {
2458            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
2459            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
2460            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
2461            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
2462
2463            double dval[6];
2464            int nread = 0;
2465            int ix = 0;
2466            int iy = 0;
2467            int iz = 0;
2468            while (!fin.eof() && nread < npts) {
2469                fin.getline(line,sizeof(line)-1);
2470                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
2471
2472                for (int p=0; p < n; p++) {
2473                    int nindex = iz*nx*ny + iy*nx + ix;
2474                    field.define(nindex, dval[p]);
2475                    nread++;
2476                    if (++iz >= nz) {
2477                        iz = 0;
2478                        if (++iy >= ny) {
2479                            iy = 0;
2480                            ++ix;
2481                        }
2482                    }
2483                }
2484            }
2485
2486            // make sure that we read all of the expected points
2487            if (nread != nx*ny*nz) {
2488                char mesg[256];
2489                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
2490                result.error(mesg);
2491                return result;
2492            }
2493
2494            // figure out a good mesh spacing
2495            int nsample = 30;
2496            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
2497            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
2498            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
2499            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
2500
2501            nx = (int)ceil(dx/dmin);
2502            ny = (int)ceil(dy/dmin);
2503            nz = (int)ceil(dz/dmin);
2504
2505#ifndef NV40
2506            // must be an even power of 2 for older cards
2507                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
2508                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
2509                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
2510#endif
2511
2512            float *data = new float[4*nx*ny*nz];
2513
2514            double vmin = field.valueMin();
2515            double dv = field.valueMax() - field.valueMin();
2516            if (dv == 0.0) { dv = 1.0; }
2517
2518            // generate the uniformly sampled data that we need for a volume
2519            int ngen = 0;
2520            double nzero_min = 0.0;
2521            for (int iz=0; iz < nz; iz++) {
2522                double zval = z0 + iz*dmin;
2523                for (int iy=0; iy < ny; iy++) {
2524                    double yval = y0 + iy*dmin;
2525                    for (int ix=0; ix < nx; ix++) {
2526                        double xval = x0 + ix*dmin;
2527                        double v = field.value(xval,yval,zval);
2528
2529                        if (v != 0.0f && v < nzero_min)
2530                        {
2531                            nzero_min = v;
2532                        }
2533
2534                        // scale all values [0-1], -1 => out of bounds
2535                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2536
2537                        data[ngen] = v;
2538                        ngen += 4;
2539                    }
2540                }
2541            }
2542
2543            // Compute the gradient of this data.  BE CAREFUL: center
2544            // calculation on each node to avoid skew in either direction.
2545            ngen = 0;
2546            for (int iz=0; iz < nz; iz++) {
2547                for (int iy=0; iy < ny; iy++) {
2548                    for (int ix=0; ix < nx; ix++) {
2549                        // gradient in x-direction
2550                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
2551                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
2552                        if (valm1 < 0 || valp1 < 0) {
2553                            data[ngen+1] = 0.0;
2554                        } else {
2555                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2556                        }
2557
2558                        // gradient in y-direction
2559                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
2560                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
2561                        if (valm1 < 0 || valp1 < 0) {
2562                            data[ngen+2] = 0.0;
2563                        } else {
2564                            //data[ngen+2] = valp1-valm1; // assume dy=1
2565                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2566                        }
2567
2568                        // gradient in z-direction
2569                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
2570                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
2571                        if (valm1 < 0 || valp1 < 0) {
2572                            data[ngen+3] = 0.0;
2573                        } else {
2574                            //data[ngen+3] = valp1-valm1; // assume dz=1
2575                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
2576                        }
2577
2578                        ngen += 4;
2579                    }
2580                }
2581            }
2582
2583            load_volume(index, nx, ny, nz, 4, data,
2584                field.valueMin(), field.valueMax(), nzero_min);
2585
2586            delete [] data;
2587
2588        } else {
2589            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
2590            Rappture::FieldPrism3D field(xymesh, zgrid);
2591
2592            double dval;
2593            int nread = 0;
2594            int ixy = 0;
2595            int iz = 0;
2596            while (!fin.eof() && nread < npts) {
2597                if (!(fin >> dval).fail()) {
2598                    int nid = nxy*iz + ixy;
2599                    field.define(nid, dval);
2600
2601                    nread++;
2602                    if (++iz >= nz) {
2603                        iz = 0;
2604                        ixy++;
2605                    }
2606                }
2607            }
2608
2609            // make sure that we read all of the expected points
2610            if (nread != nxy*nz) {
2611                char mesg[256];
2612                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
2613                return result.error(mesg);
2614            }
2615
2616            // figure out a good mesh spacing
2617            int nsample = 30;
2618            x0 = field.rangeMin(Rappture::xaxis);
2619            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
2620            y0 = field.rangeMin(Rappture::yaxis);
2621            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
2622            z0 = field.rangeMin(Rappture::zaxis);
2623            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
2624            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
2625
2626            nx = (int)ceil(dx/dmin);
2627            ny = (int)ceil(dy/dmin);
2628            nz = (int)ceil(dz/dmin);
2629#ifndef NV40
2630            // must be an even power of 2 for older cards
2631                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
2632                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
2633                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
2634#endif
2635            float *data = new float[4*nx*ny*nz];
2636
2637            double vmin = field.valueMin();
2638            double dv = field.valueMax() - field.valueMin();
2639            if (dv == 0.0) { dv = 1.0; }
2640
2641            // generate the uniformly sampled data that we need for a volume
2642            int ngen = 0;
2643            double nzero_min = 0.0;
2644            for (iz=0; iz < nz; iz++) {
2645                double zval = z0 + iz*dmin;
2646                for (int iy=0; iy < ny; iy++) {
2647                    double yval = y0 + iy*dmin;
2648                    for (int ix=0; ix < nx; ix++) {
2649                        double xval = x0 + ix*dmin;
2650                        double v = field.value(xval,yval,zval);
2651
2652                        if (v != 0.0f && v < nzero_min)
2653                        {
2654                            nzero_min = v;
2655                        }
2656                        // scale all values [0-1], -1 => out of bounds
2657                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2658                        data[ngen] = v;
2659
2660                        ngen += 4;
2661                    }
2662                }
2663            }
2664
2665            // Compute the gradient of this data.  BE CAREFUL: center
2666            // calculation on each node to avoid skew in either direction.
2667            ngen = 0;
2668            for (int iz=0; iz < nz; iz++) {
2669                for (int iy=0; iy < ny; iy++) {
2670                    for (int ix=0; ix < nx; ix++) {
2671                        // gradient in x-direction
2672                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
2673                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
2674                        if (valm1 < 0 || valp1 < 0) {
2675                            data[ngen+1] = 0.0;
2676                        } else {
2677                            //data[ngen+1] = valp1-valm1; // assume dx=1
2678                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2679                        }
2680
2681                        // gradient in y-direction
2682                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
2683                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
2684                        if (valm1 < 0 || valp1 < 0) {
2685                            data[ngen+2] = 0.0;
2686                        } else {
2687                            //data[ngen+2] = valp1-valm1; // assume dy=1
2688                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
2689                        }
2690
2691                        // gradient in z-direction
2692                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
2693                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
2694                        if (valm1 < 0 || valp1 < 0) {
2695                            data[ngen+3] = 0.0;
2696                        } else {
2697                            //data[ngen+3] = valp1-valm1; // assume dz=1
2698                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
2699                        }
2700
2701                        ngen += 4;
2702                    }
2703                }
2704            }
2705
2706            load_volume(index, nx, ny, nz, 4, data,
2707                field.valueMin(), field.valueMax(), nzero_min);
2708
2709            delete [] data;
2710        }
2711    } else {
2712        return result.error("data not found in stream");
2713    }
2714
2715    //
2716    // Center this new volume on the origin.
2717    //
2718    float dx0 = -0.5;
2719    float dy0 = -0.5*dy/dx;
2720    float dz0 = -0.5*dz/dx;
2721    volume[index]->move(Vector3(dx0, dy0, dz0));
2722
2723    return result;
2724}
2725
2726Rappture::Outcome
2727load_volume_stream3(int index, std::iostream& fin)
2728{
2729    Rappture::Outcome result;
2730
2731    Rappture::MeshTri2D xymesh;
2732    int dummy, nx, ny, nz, nxy, npts;
2733    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
2734    char line[128], type[128], *start;
2735
2736    int isrect = 1;
2737
2738    float temp1, temp2, temp3;
2739    float vmin, vmax;
2740    int count;
2741
2742    char temp[24];
2743    do {
2744        fin.getline(line,sizeof(line)-1);
2745
2746        // skip leading blanks
2747        for (start=&line[0]; *start == ' ' || *start == '\t'; start++);
2748
2749        if (*start != '#')
2750        { 
2751            sscanf(start, "%s%f%f%f", temp, &temp1, &temp2, &temp3);
2752            if (!strcmp(temp, "size"))
2753            {
2754                dx = temp1;
2755                dy = temp2;
2756                dz = temp3;
2757            }
2758            else if (!strcmp(temp, "size"))
2759            {
2760                x0 = temp1;
2761                y0 = temp2;
2762                z0 = temp3;
2763            }
2764            else if (!strcmp(temp, "minmax"))
2765            {
2766                vmin = temp1;
2767                vmax = temp2;
2768            }
2769            else if (!strcmp(temp, "count"))
2770            {
2771                count = temp1;
2772            }
2773            else if (!strcmp(temp, "<FET>"))
2774            {
2775                break;
2776            }
2777        }
2778    } while (!fin.eof());
2779
2780    // read data points
2781    if (!fin.eof()) {
2782        if (isrect) {
2783            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
2784            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
2785            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
2786            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
2787
2788            double dval[6];
2789            int nread = 0;
2790            int ix = 0;
2791            int iy = 0;
2792            int iz = 0;
2793            while (!fin.eof() && nread < npts) {
2794                fin.getline(line,sizeof(line)-1);
2795                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
2796
2797                for (int p=0; p < n; p++) {
2798                    int nindex = iz*nx*ny + iy*nx + ix;
2799                    field.define(nindex, dval[p]);
2800                    nread++;
2801                    if (++iz >= nz) {
2802                        iz = 0;
2803                        if (++iy >= ny) {
2804                            iy = 0;
2805                            ++ix;
2806                        }
2807                    }
2808                }
2809            }
2810
2811            // make sure that we read all of the expected points
2812            if (nread != nx*ny*nz) {
2813                char mesg[256];
2814                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
2815                result.error(mesg);
2816                return result;
2817            }
2818
2819            // figure out a good mesh spacing
2820            int nsample = 30;
2821            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
2822            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
2823            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
2824            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
2825
2826            nx = (int)ceil(dx/dmin);
2827            ny = (int)ceil(dy/dmin);
2828            nz = (int)ceil(dz/dmin);
2829
2830#ifndef NV40
2831            // must be an even power of 2 for older cards
2832                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
2833                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
2834                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
2835#endif
2836
2837            float *data = new float[4*nx*ny*nz];
2838
2839            double vmin = field.valueMin();
2840            double dv = field.valueMax() - field.valueMin();
2841            if (dv == 0.0) { dv = 1.0; }
2842
2843            // generate the uniformly sampled data that we need for a volume
2844            int ngen = 0;
2845            double nzero_min = 0.0;
2846            for (int iz=0; iz < nz; iz++) {
2847                double zval = z0 + iz*dmin;
2848                for (int iy=0; iy < ny; iy++) {
2849                    double yval = y0 + iy*dmin;
2850                    for (int ix=0; ix < nx; ix++) {
2851                        double xval = x0 + ix*dmin;
2852                        double v = field.value(xval,yval,zval);
2853
2854                        if (v != 0.0f && v < nzero_min)
2855                        {
2856                            nzero_min = v;
2857                        }
2858
2859                        // scale all values [0-1], -1 => out of bounds
2860                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2861
2862                        data[ngen] = v;
2863                        ngen += 4;
2864                    }
2865                }
2866            }
2867
2868            // Compute the gradient of this data.  BE CAREFUL: center
2869            // calculation on each node to avoid skew in either direction.
2870            ngen = 0;
2871            for (int iz=0; iz < nz; iz++) {
2872                for (int iy=0; iy < ny; iy++) {
2873                    for (int ix=0; ix < nx; ix++) {
2874                        // gradient in x-direction
2875                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
2876                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
2877                        if (valm1 < 0 || valp1 < 0) {
2878                            data[ngen+1] = 0.0;
2879                        } else {
2880                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2881                        }
2882
2883                        // gradient in y-direction
2884                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
2885                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
2886                        if (valm1 < 0 || valp1 < 0) {
2887                            data[ngen+2] = 0.0;
2888                        } else {
2889                            //data[ngen+2] = valp1-valm1; // assume dy=1
2890                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
2891                        }
2892
2893                        // gradient in z-direction
2894                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
2895                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
2896                        if (valm1 < 0 || valp1 < 0) {
2897                            data[ngen+3] = 0.0;
2898                        } else {
2899                            //data[ngen+3] = valp1-valm1; // assume dz=1
2900                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
2901                        }
2902
2903                        ngen += 4;
2904                    }
2905                }
2906            }
2907
2908            load_volume(index, nx, ny, nz, 4, data,
2909                field.valueMin(), field.valueMax(), nzero_min);
2910
2911            delete [] data;
2912
2913        } else {
2914            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
2915            Rappture::FieldPrism3D field(xymesh, zgrid);
2916
2917            double dval;
2918            int nread = 0;
2919            int ixy = 0;
2920            int iz = 0;
2921            while (!fin.eof() && nread < npts) {
2922                if (!(fin >> dval).fail()) {
2923                    int nid = nxy*iz + ixy;
2924                    field.define(nid, dval);
2925
2926                    nread++;
2927                    if (++iz >= nz) {
2928                        iz = 0;
2929                        ixy++;
2930                    }
2931                }
2932            }
2933
2934            // make sure that we read all of the expected points
2935            if (nread != nxy*nz) {
2936                char mesg[256];
2937                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
2938                return result.error(mesg);
2939            }
2940
2941            // figure out a good mesh spacing
2942            int nsample = 30;
2943            x0 = field.rangeMin(Rappture::xaxis);
2944            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
2945            y0 = field.rangeMin(Rappture::yaxis);
2946            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
2947            z0 = field.rangeMin(Rappture::zaxis);
2948            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
2949            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
2950
2951            nx = (int)ceil(dx/dmin);
2952            ny = (int)ceil(dy/dmin);
2953            nz = (int)ceil(dz/dmin);
2954#ifndef NV40
2955            // must be an even power of 2 for older cards
2956                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
2957                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
2958                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
2959#endif
2960            float *data = new float[4*nx*ny*nz];
2961
2962            double vmin = field.valueMin();
2963            double dv = field.valueMax() - field.valueMin();
2964            if (dv == 0.0) { dv = 1.0; }
2965
2966            // generate the uniformly sampled data that we need for a volume
2967            int ngen = 0;
2968            double nzero_min = 0.0;
2969            for (iz=0; iz < nz; iz++) {
2970                double zval = z0 + iz*dmin;
2971                for (int iy=0; iy < ny; iy++) {
2972                    double yval = y0 + iy*dmin;
2973                    for (int ix=0; ix < nx; ix++) {
2974                        double xval = x0 + ix*dmin;
2975                        double v = field.value(xval,yval,zval);
2976
2977                        if (v != 0.0f && v < nzero_min)
2978                        {
2979                            nzero_min = v;
2980                        }
2981                        // scale all values [0-1], -1 => out of bounds
2982                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
2983                        data[ngen] = v;
2984
2985                        ngen += 4;
2986                    }
2987                }
2988            }
2989
2990            // Compute the gradient of this data.  BE CAREFUL: center
2991            // calculation on each node to avoid skew in either direction.
2992            ngen = 0;
2993            for (int iz=0; iz < nz; iz++) {
2994                for (int iy=0; iy < ny; iy++) {
2995                    for (int ix=0; ix < nx; ix++) {
2996                        // gradient in x-direction
2997                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
2998                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
2999                        if (valm1 < 0 || valp1 < 0) {
3000                            data[ngen+1] = 0.0;
3001                        } else {
3002                            //data[ngen+1] = valp1-valm1; // assume dx=1
3003                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
3004                        }
3005
3006                        // gradient in y-direction
3007                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
3008                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
3009                        if (valm1 < 0 || valp1 < 0) {
3010                            data[ngen+2] = 0.0;
3011                        } else {
3012                            //data[ngen+2] = valp1-valm1; // assume dy=1
3013                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
3014                        }
3015
3016                        // gradient in z-direction
3017                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
3018                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
3019                        if (valm1 < 0 || valp1 < 0) {
3020                            data[ngen+3] = 0.0;
3021                        } else {
3022                            //data[ngen+3] = valp1-valm1; // assume dz=1
3023                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
3024                        }
3025
3026                        ngen += 4;
3027                    }
3028                }
3029            }
3030
3031            load_volume(index, nx, ny, nz, 4, data,
3032                field.valueMin(), field.valueMax(), nzero_min);
3033
3034            delete [] data;
3035        }
3036    } else {
3037        return result.error("data not found in stream");
3038    }
3039
3040    //
3041    // Center this new volume on the origin.
3042    //
3043    float dx0 = -0.5;
3044    float dy0 = -0.5*dy/dx;
3045    float dz0 = -0.5*dz/dx;
3046    volume[index]->move(Vector3(dx0, dy0, dz0));
3047
3048    return result;
3049}
3050
3051/* Load a 3D volume
3052 * index: the index into the volume array, which stores pointers to 3D volume instances
3053 * data: pointer to an array of floats. 
3054 * n_component: the number of scalars for each space point.
3055 *              All component scalars for a point are placed consequtively in data array
3056 * width, height and depth: number of points in each dimension
3057 */
3058void load_volume(int index, int width, int height, int depth,
3059    int n_component, float* data, double vmin, double vmax, double nzero_min)
3060{
3061    while (n_volumes <= index) {
3062        volume.push_back(NULL);
3063        n_volumes++;
3064    }
3065
3066    if (volume[index] != NULL) {
3067        delete volume[index];
3068        volume[index] = NULL;
3069    }
3070
3071    volume[index] = new Volume(0.f, 0.f, 0.f, width, height, depth, 1.,
3072                                 n_component, data, vmin, vmax, nzero_min);
3073    assert(volume[index]!=0);
3074}
3075
3076// get a colormap 1D texture by name
3077TransferFunction*
3078get_transfunc(char *name) {
3079    Tcl_HashEntry *entryPtr;
3080    entryPtr = Tcl_FindHashEntry(&tftable, name);
3081    if (entryPtr) {
3082        return (TransferFunction*)Tcl_GetHashValue(entryPtr);
3083    }
3084    return NULL;
3085}
3086
3087
3088//Update the transfer function using local GUI in the non-server mode
3089void update_tf_texture()
3090{
3091  glutSetWindow(render_window);
3092
3093  //fprintf(stderr, "tf update\n");
3094  TransferFunction *tf = get_transfunc("default");
3095  if (tf == NULL) return;
3096
3097  float data[256*4];
3098  for(int i=0; i<256; i++){
3099    data[4*i+0] = color_table[i][0];
3100    data[4*i+1] = color_table[i][1];
3101    data[4*i+2] = color_table[i][2];
3102    data[4*i+3] = color_table[i][3];
3103    //fprintf(stderr, "(%f,%f,%f,%f) ", data[4*i+0], data[4*i+1], data[4*i+2], data[4*i+3]);
3104  }
3105
3106  tf->update(data);
3107
3108#ifdef EVENTLOG
3109  float param[3] = {0,0,0};
3110  Event* tmp = new Event(EVENT_ROTATE, param, NvGetTimeInterval());
3111  tmp->write(event_log);
3112  delete tmp;
3113#endif
3114}
3115
3116
3117//initialize frame buffer objects for offscreen rendering
3118void init_offscreen_buffer()
3119{
3120
3121  //initialize final fbo for final display
3122  glGenFramebuffersEXT(1, &final_fbo);
3123  glGenTextures(1, &final_color_tex);
3124  glGenRenderbuffersEXT(1, &final_depth_rb);
3125
3126  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
3127
3128  //initialize final color texture
3129  glBindTexture(GL_TEXTURE_2D, final_color_tex);
3130  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
3131  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
3132#ifdef NV40
3133  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, win_width, win_height, 0,
3134               GL_RGB, GL_INT, NULL);
3135#else
3136  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, win_width, win_height, 0,
3137               GL_RGB, GL_INT, NULL);
3138#endif
3139  glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,
3140                            GL_COLOR_ATTACHMENT0_EXT,
3141                            GL_TEXTURE_2D, final_color_tex, 0);
3142
3143  // initialize final depth renderbuffer
3144  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
3145  glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT,
3146                           GL_DEPTH_COMPONENT24, win_width, win_height);
3147  glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT,
3148                               GL_DEPTH_ATTACHMENT_EXT,
3149                               GL_RENDERBUFFER_EXT, final_depth_rb);
3150
3151  // Check framebuffer completeness at the end of initialization.
3152  CHECK_FRAMEBUFFER_STATUS();
3153  assert(glGetError()==0);
3154}
3155
3156
3157//resize the offscreen buffer
3158void resize_offscreen_buffer(int w, int h){
3159  win_width = w;
3160  win_height = h;
3161
3162    if (g_fonts)
3163    {
3164        g_fonts->resize(w, h);
3165    }
3166
3167  //fprintf(stderr, "screen_buffer size: %d\n", sizeof(screen_buffer));
3168  printf("screen_buffer size: %d %d\n", w, h);
3169
3170  if (screen_buffer) {
3171      delete[] screen_buffer;
3172      screen_buffer = NULL;
3173  }
3174
3175  screen_buffer = new unsigned char[4*win_width*win_height];
3176  assert(screen_buffer != NULL);
3177
3178  //delete the current render buffer resources
3179  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
3180  glDeleteTextures(1, &final_color_tex);
3181  glDeleteFramebuffersEXT(1, &final_fbo);
3182
3183  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
3184  glDeleteRenderbuffersEXT(1, &final_depth_rb);
3185fprintf(stdin,"  after glDeleteRenderbuffers\n");
3186
3187  //change the camera setting
3188  cam->set_screen_size(0, 0, win_width, win_height);
3189  plane_render->set_screen_size(win_width, win_height);
3190
3191  //Reinitialize final fbo for final display
3192  glGenFramebuffersEXT(1, &final_fbo);
3193  glGenTextures(1, &final_color_tex);
3194  glGenRenderbuffersEXT(1, &final_depth_rb);
3195fprintf(stdin,"  after glGenRenderbuffers\n");
3196
3197  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
3198
3199  //initialize final color texture
3200  glBindTexture(GL_TEXTURE_2D, final_color_tex);
3201  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
3202  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
3203#ifdef NV40
3204  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, win_width, win_height, 0,
3205               GL_RGB, GL_INT, NULL);
3206#else
3207  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, win_width, win_height, 0,
3208               GL_RGB, GL_INT, NULL);
3209#endif
3210  glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,
3211                            GL_COLOR_ATTACHMENT0_EXT,
3212                            GL_TEXTURE_2D, final_color_tex, 0);
3213fprintf(stdin,"  after glFramebufferTexture2DEXT\n");
3214       
3215  // initialize final depth renderbuffer
3216  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
3217  glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT,
3218                           GL_DEPTH_COMPONENT24, win_width, win_height);
3219  glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT,
3220                               GL_DEPTH_ATTACHMENT_EXT,
3221                               GL_RENDERBUFFER_EXT, final_depth_rb);
3222fprintf(stdin,"  after glFramebufferRenderEXT\n");
3223
3224  // Check framebuffer completeness at the end of initialization.
3225  CHECK_FRAMEBUFFER_STATUS();
3226  assert(glGetError()==0);
3227fprintf(stdin,"  after assert\n");
3228}
3229
3230
3231//init line integral convolution
3232void init_lic(){
3233  lic = new Lic(NMESH, win_width, win_height, 0.3, g_context, volume[1]->id,
3234                  volume[1]->aspect_ratio_width,
3235                  volume[1]->aspect_ratio_height,
3236                  volume[1]->aspect_ratio_depth);
3237}
3238
3239//init the particle system using vector field volume->[1]
3240/*
3241void init_particle_system()
3242{
3243    psys = new ParticleSystem(NMESH, NMESH, g_context, volume[1]->id,
3244        1./volume[1]->aspect_ratio_width,
3245        1./volume[1]->aspect_ratio_height,
3246        1./volume[1]->aspect_ratio_depth);
3247
3248    init_particles();   //fill initial particles
3249}
3250*/
3251
3252
3253void make_test_2D_data()
3254{
3255
3256  int w = 300;
3257  int h = 200;
3258  float* data = new float[w*h];
3259
3260  //procedurally make a gradient plane
3261  for(int j=0; j<h; j++){
3262    for(int i=0; i<w; i++){
3263      data[w*j+i] = float(i)/float(w);
3264    }
3265  }
3266
3267  plane[0] = new Texture2D(w, h, GL_FLOAT, GL_LINEAR, 1, data);
3268
3269  delete[] data;
3270}
3271
3272
3273/*----------------------------------------------------*/
3274void initGL(void)
3275{
3276   //buffer to store data read from the screen
3277   if (screen_buffer) {
3278       delete[] screen_buffer;
3279       screen_buffer = NULL;
3280   }
3281   screen_buffer = new unsigned char[4*win_width*win_height];
3282   assert(screen_buffer != NULL);
3283
3284   //create the camera with default setting
3285   cam = new Camera(0, 0, win_width, win_height,
3286                   live_obj_x, live_obj_y, live_obj_z,
3287                   0., 0., 100.,
3288                   (int)live_rot_x, (int)live_rot_y, (int)live_rot_z);
3289
3290   glEnable(GL_TEXTURE_2D);
3291   glShadeModel(GL_FLAT);
3292   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
3293
3294   glClearColor(0,0,0,1);
3295   glClear(GL_COLOR_BUFFER_BIT);
3296
3297   //initialize lighting
3298   GLfloat mat_specular[] = {1.0, 1.0, 1.0, 1.0};
3299   GLfloat mat_shininess[] = {30.0};
3300   GLfloat white_light[] = {1.0, 1.0, 1.0, 1.0};
3301   GLfloat green_light[] = {0.1, 0.5, 0.1, 1.0};
3302
3303   glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
3304   glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
3305   glLightfv(GL_LIGHT0, GL_DIFFUSE, white_light);
3306   glLightfv(GL_LIGHT0, GL_SPECULAR, white_light);     
3307   glLightfv(GL_LIGHT1, GL_DIFFUSE, green_light);
3308   glLightfv(GL_LIGHT1, GL_SPECULAR, white_light);     
3309
3310   // init table of transfer functions
3311   Tcl_InitHashTable(&tftable, TCL_STRING_KEYS);
3312
3313   //check if performance query is supported
3314   if(check_query_support()){
3315     //create queries to count number of rendered pixels
3316     perf = new PerfQuery();
3317   }
3318
3319   init_offscreen_buffer();    //frame buffer object for offscreen rendering
3320
3321   //create volume renderer and add volumes to it
3322   g_vol_render = new VolumeRenderer();
3323
3324   /*
3325   //I added this to debug : Wei
3326   float tmp_data[4*124];
3327   memset(tmp_data, 0, 4*4*124);
3328   TransferFunction* tmp_tf = new TransferFunction(124, tmp_data);
3329   g_vol_render->add_volume(volume[0], tmp_tf);
3330   volume[0]->get_cutplane(0)->enabled = false;
3331   volume[0]->get_cutplane(1)->enabled = false;
3332   volume[0]->get_cutplane(2)->enabled = false;
3333
3334   //volume[1]->move(Vector3(0.5, 0.6, 0.7));
3335   //g_vol_render->add_volume(volume[1], tmp_tf);
3336   */
3337
3338   //create an 2D plane renderer
3339   plane_render = new PlaneRenderer(g_context, win_width, win_height);
3340   make_test_2D_data();
3341
3342   plane_render->add_plane(plane[0], get_transfunc("default"));
3343
3344   assert(glGetError()==0);
3345
3346   //init_particle_system();
3347   //init_lic();
3348}
3349
3350
3351
3352void initTcl()
3353{
3354    interp = Tcl_CreateInterp();
3355    Tcl_MakeSafe(interp);
3356
3357    Tcl_DStringInit(&cmdbuffer);
3358
3359    // manipulate the viewpoint
3360    Tcl_CreateCommand(interp, "camera", CameraCmd,
3361        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3362
3363    // turn on/off cut planes in x/y/z directions
3364    Tcl_CreateCommand(interp, "cutplane", CutplaneCmd,
3365        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3366
3367    // request the legend for a plot (transfer function)
3368    Tcl_CreateCommand(interp, "legend", LegendCmd,
3369        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3370
3371    // change the size of the screen (size of picture generated)
3372    Tcl_CreateCommand(interp, "screen", ScreenCmd,
3373        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3374
3375    // manipulate transfer functions
3376    Tcl_CreateCommand(interp, "transfunc", TransfuncCmd,
3377        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3378
3379    // set the "up" direction for volumes
3380    Tcl_CreateCommand(interp, "up", UpCmd,
3381        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3382
3383    // manipulate volume data
3384    Tcl_CreateCommand(interp, "volume", VolumeCmd,
3385        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3386
3387    Tcl_CreateCommand(interp, "axis", AxisCmd,
3388        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3389
3390    Tcl_CreateCommand(interp, "grid", GridCmd,
3391        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3392
3393    Tcl_CreateCommand(interp, "heightmap", HeightMapCmd,
3394        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3395
3396    // get screenshot
3397    Tcl_CreateCommand(interp, "screenshot", ScreenShotCmd,
3398        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3399
3400    Tcl_CreateCommand(interp, "test", TestCmd,
3401        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
3402
3403    // create a default transfer function
3404    if (Tcl_Eval(interp, def_transfunc) != TCL_OK) {
3405        fprintf(stdin, "WARNING: bad default transfer function\n");
3406        fprintf(stdin, Tcl_GetStringResult(interp));
3407    }
3408}
3409
3410
3411void read_screen()
3412{
3413  glReadPixels(0, 0, win_width, win_height, GL_RGB, GL_UNSIGNED_BYTE, screen_buffer);
3414  assert(glGetError()==0);
3415}
3416
3417void display();
3418
3419
3420#if DO_RLE
3421char rle[512*512*3];
3422int rle_size;
3423
3424short offsets[512*512*3];
3425int offsets_size;
3426
3427void do_rle(){
3428  int len = win_width*win_height*3;
3429  rle_size = 0;
3430  offsets_size = 0;
3431
3432  int i=0;
3433  while(i<len){
3434    if (screen_buffer[i] == 0) {
3435      int pos = i+1;
3436      while ( (pos<len) && (screen_buffer[pos] == 0)) {
3437        pos++;
3438      }
3439      offsets[offsets_size++] = -(pos - i);
3440      i = pos;
3441    }
3442
3443    else {
3444      int pos;
3445      for (pos = i; (pos<len) && (screen_buffer[pos] != 0); pos++) {
3446        rle[rle_size++] = screen_buffer[pos];
3447      }
3448      offsets[offsets_size++] = (pos - i);
3449      i = pos;
3450    }
3451
3452  }
3453}
3454#endif
3455
3456// used internally to build up the BMP file header
3457// writes an integer value into the header data structure at pos
3458void
3459bmp_header_add_int(unsigned char* header, int& pos, int data)
3460{
3461    header[pos++] = data & 0xff;
3462    header[pos++] = (data >> 8) & 0xff;
3463    header[pos++] = (data >> 16) & 0xff;
3464    header[pos++] = (data >> 24) & 0xff;
3465}
3466
3467// INSOO
3468// FOR DEBUGGING
3469void
3470bmp_write_to_file()
3471{
3472    unsigned char header[54];
3473    int pos = 0;
3474    header[pos++] = 'B';
3475    header[pos++] = 'M';
3476
3477    // BE CAREFUL:  BMP files must have an even multiple of 4 bytes
3478    // on each scan line.  If need be, we add padding to each line.
3479    int pad = 0;
3480    if ((3*win_width) % 4 > 0) {
3481        pad = 4 - ((3*win_width) % 4);
3482    }
3483
3484    // file size in bytes
3485    int fsize = (3*win_width+pad)*win_height + sizeof(header);
3486    bmp_header_add_int(header, pos, fsize);
3487
3488    // reserved value (must be 0)
3489    bmp_header_add_int(header, pos, 0);
3490
3491    // offset in bytes to start of bitmap data
3492    bmp_header_add_int(header, pos, sizeof(header));
3493
3494    // size of the BITMAPINFOHEADER
3495    bmp_header_add_int(header, pos, 40);
3496
3497    // width of the image in pixels
3498    bmp_header_add_int(header, pos, win_width);
3499
3500    // height of the image in pixels
3501    bmp_header_add_int(header, pos, win_height);
3502
3503    // 1 plane + (24 bits/pixel << 16)
3504    bmp_header_add_int(header, pos, 1572865);
3505
3506    // no compression
3507    // size of image for compression
3508    bmp_header_add_int(header, pos, 0);
3509    bmp_header_add_int(header, pos, 0);
3510
3511    // x pixels per meter
3512    // y pixels per meter
3513    bmp_header_add_int(header, pos, 0);
3514    bmp_header_add_int(header, pos, 0);
3515
3516    // number of colors used (0 = compute from bits/pixel)
3517    // number of important colors (0 = all colors important)
3518    bmp_header_add_int(header, pos, 0);
3519    bmp_header_add_int(header, pos, 0);
3520
3521    // BE CAREFUL: BMP format wants BGR ordering for screen data
3522    unsigned char* scr = screen_buffer;
3523    for (int row=0; row < win_height; row++) {
3524        for (int col=0; col < win_width; col++) {
3525            unsigned char tmp = scr[2];
3526            scr[2] = scr[0];  // B
3527            scr[0] = tmp;     // R
3528            scr += 3;
3529        }
3530        scr += pad;  // skip over padding already in screen data
3531    }
3532
3533    FILE* f;
3534    f = fopen("/tmp/image.bmp", "wb");
3535    fwrite((void*) header, sizeof(header), 1, f);
3536    fwrite((void*) screen_buffer, (3*win_width+pad)*win_height, 1, f);
3537    fclose(f);
3538}
3539
3540void
3541bmp_write(const char* cmd)
3542{
3543    unsigned char header[54];
3544    int pos = 0;
3545    header[pos++] = 'B';
3546    header[pos++] = 'M';
3547
3548    // BE CAREFUL:  BMP files must have an even multiple of 4 bytes
3549    // on each scan line.  If need be, we add padding to each line.
3550    int pad = 0;
3551    if ((3*win_width) % 4 > 0) {
3552        pad = 4 - ((3*win_width) % 4);
3553    }
3554
3555    // file size in bytes
3556    int fsize = (3*win_width+pad)*win_height + sizeof(header);
3557    bmp_header_add_int(header, pos, fsize);
3558
3559    // reserved value (must be 0)
3560    bmp_header_add_int(header, pos, 0);
3561
3562    // offset in bytes to start of bitmap data
3563    bmp_header_add_int(header, pos, sizeof(header));
3564
3565    // size of the BITMAPINFOHEADER
3566    bmp_header_add_int(header, pos, 40);
3567
3568    // width of the image in pixels
3569    bmp_header_add_int(header, pos, win_width);
3570
3571    // height of the image in pixels
3572    bmp_header_add_int(header, pos, win_height);
3573
3574    // 1 plane + (24 bits/pixel << 16)
3575    bmp_header_add_int(header, pos, 1572865);
3576
3577    // no compression
3578    // size of image for compression
3579    bmp_header_add_int(header, pos, 0);
3580    bmp_header_add_int(header, pos, 0);
3581
3582    // x pixels per meter
3583    // y pixels per meter
3584    bmp_header_add_int(header, pos, 0);
3585    bmp_header_add_int(header, pos, 0);
3586
3587    // number of colors used (0 = compute from bits/pixel)
3588    // number of important colors (0 = all colors important)
3589    bmp_header_add_int(header, pos, 0);
3590    bmp_header_add_int(header, pos, 0);
3591
3592    // BE CAREFUL: BMP format wants BGR ordering for screen data
3593    unsigned char* scr = screen_buffer;
3594    for (int row=0; row < win_height; row++) {
3595        for (int col=0; col < win_width; col++) {
3596            unsigned char tmp = scr[2];
3597            scr[2] = scr[0];  // B
3598            scr[0] = tmp;     // R
3599            scr += 3;
3600        }
3601        scr += pad;  // skip over padding already in screen data
3602    }
3603
3604    std::ostringstream result;
3605    result << cmd << " " << fsize << "\n";
3606    write(0, result.str().c_str(), result.str().size());
3607
3608    write(0, header, sizeof(header));
3609    write(0, screen_buffer, (3*win_width+pad)*win_height);
3610
3611}
3612
3613
3614void xinetd_listen(){
3615
3616    int flags = fcntl(0, F_GETFL, 0);
3617    fcntl(0, F_SETFL, flags & ~O_NONBLOCK);
3618
3619    int status = TCL_OK;
3620    int npass = 0;
3621
3622    //
3623    //  Read and execute as many commands as we can from stdin...
3624    //
3625    while (status == TCL_OK) {
3626        //
3627        //  Read the next command from the buffer.  First time through
3628        //  we block here and wait if necessary until a command comes in.
3629        //
3630        //  BE CAREFUL:  Read only one command, up to a newline.
3631        //  The "volume data follows" command needs to be able to read
3632        //  the data immediately following the command, and we shouldn't
3633        //  consume it here.
3634        //
3635        while (1) {
3636            char c = getchar();
3637            if (c <= 0) {
3638                if (npass == 0) {
3639                    exit(0);  // EOF -- we're done!
3640                } else {
3641                    break;
3642                }
3643            }
3644            Tcl_DStringAppend(&cmdbuffer, &c, 1);
3645
3646            if (c=='\n' && Tcl_CommandComplete(Tcl_DStringValue(&cmdbuffer))) {
3647                break;
3648            }
3649        }
3650
3651        // no command? then we're done for now
3652        if (Tcl_DStringLength(&cmdbuffer) == 0) {
3653            break;
3654        }
3655
3656        // back to original flags during command evaluation...
3657        fcntl(0, F_SETFL, flags & ~O_NONBLOCK);
3658
3659        status = Tcl_Eval(interp, Tcl_DStringValue(&cmdbuffer));
3660        Tcl_DStringSetLength(&cmdbuffer, 0);
3661
3662        // non-blocking for next read -- we might not get anything
3663        fcntl(0, F_SETFL, flags | O_NONBLOCK);
3664        npass++;
3665    }
3666    fcntl(0, F_SETFL, flags);
3667
3668    if (status != TCL_OK) {
3669        std::ostringstream errmsg;
3670        errmsg << "ERROR: " << Tcl_GetStringResult(interp) << std::endl;
3671        write(0, errmsg.str().c_str(), errmsg.str().size());
3672        return;
3673    }
3674
3675    //
3676    //  Generate the latest frame and send it back to the client
3677    //
3678    // INSOO
3679    offscreen_buffer_capture();  //enable offscreen render
3680
3681    display();
3682
3683    // INSOO
3684#ifdef XINETD
3685   read_screen();
3686   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
3687#else
3688   display_offscreen_buffer(); //display the final rendering on screen
3689   read_screen();
3690   glutSwapBuffers();
3691#endif   
3692
3693#if DO_RLE
3694    do_rle();
3695    int sizes[2] = {  offsets_size*sizeof(offsets[0]), rle_size };
3696    fprintf(stderr, "Writing %d,%d\n", sizes[0], sizes[1]); fflush(stderr);
3697    write(0, &sizes, sizeof(sizes));
3698    write(0, offsets, offsets_size*sizeof(offsets[0]));
3699    write(0, rle, rle_size);    //unsigned byte
3700#else
3701    bmp_write("nv>image -bytes");
3702#endif
3703}
3704
3705
3706/*
3707//draw vectors
3708void draw_arrows(){
3709  glColor4f(0.,0.,1.,1.);
3710  for(int i=0; i<NMESH; i++){
3711    for(int j=0; j<NMESH; j++){
3712      Vector2 v = grid.get(i, j);
3713
3714      int x1 = i*DM;
3715      int y1 = j*DM;
3716
3717      int x2 = x1 + v.x;
3718      int y2 = y1 + v.y;
3719
3720      glBegin(GL_LINES);
3721        glVertex2d(x1, y1);
3722        glVertex2d(x2, y2);
3723      glEnd();
3724    }
3725  }
3726}
3727*/
3728
3729
3730/*----------------------------------------------------*/
3731void idle()
3732{
3733    glutSetWindow(render_window);
3734 
3735  /*
3736  struct timespec ts;
3737  ts.tv_sec = 0;
3738  ts.tv_nsec = 300000000;
3739  nanosleep(&ts, 0);
3740  */
3741
3742#ifdef XINETD
3743    xinetd_listen();
3744#else
3745    glutPostRedisplay();
3746#endif
3747}
3748
3749
3750void display_offscreen_buffer()
3751{
3752   glEnable(GL_TEXTURE_2D);
3753   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
3754   glBindTexture(GL_TEXTURE_2D, final_color_tex);
3755
3756   glViewport(0, 0, win_width, win_height);
3757   glMatrixMode(GL_PROJECTION);
3758   glLoadIdentity();
3759   gluOrtho2D(0, win_width, 0, win_height);
3760   glMatrixMode(GL_MODELVIEW);
3761   glLoadIdentity();
3762
3763   glColor3f(1.,1.,1.);         //MUST HAVE THIS LINE!!!
3764   glBegin(GL_QUADS);
3765   glTexCoord2f(0, 0); glVertex2f(0, 0);
3766   glTexCoord2f(1, 0); glVertex2f(win_width, 0);
3767   glTexCoord2f(1, 1); glVertex2f(win_width, win_height);
3768   glTexCoord2f(0, 1); glVertex2f(0, win_height);
3769   glEnd();
3770
3771}
3772
3773void offscreen_buffer_capture()
3774{
3775   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
3776}
3777
3778void draw_bounding_box(float x0, float y0, float z0,
3779                float x1, float y1, float z1,
3780                float r, float g, float b, float line_width)
3781{
3782        glDisable(GL_TEXTURE_2D);
3783
3784        glColor4d(r, g, b, 1.0);
3785        glLineWidth(line_width);
3786       
3787        glBegin(GL_LINE_LOOP);
3788
3789                glVertex3d(x0, y0, z0);
3790                glVertex3d(x1, y0, z0);
3791                glVertex3d(x1, y1, z0);
3792                glVertex3d(x0, y1, z0);
3793               
3794        glEnd();
3795
3796        glBegin(GL_LINE_LOOP);
3797
3798                glVertex3d(x0, y0, z1);
3799                glVertex3d(x1, y0, z1);
3800                glVertex3d(x1, y1, z1);
3801                glVertex3d(x0, y1, z1);
3802               
3803        glEnd();
3804
3805
3806        glBegin(GL_LINE_LOOP);
3807
3808                glVertex3d(x0, y0, z0);
3809                glVertex3d(x0, y0, z1);
3810                glVertex3d(x0, y1, z1);
3811                glVertex3d(x0, y1, z0);
3812               
3813        glEnd();
3814
3815        glBegin(GL_LINE_LOOP);
3816
3817                glVertex3d(x1, y0, z0);
3818                glVertex3d(x1, y0, z1);
3819                glVertex3d(x1, y1, z1);
3820                glVertex3d(x1, y1, z0);
3821               
3822        glEnd();
3823
3824        glEnable(GL_TEXTURE_2D);
3825}
3826
3827
3828
3829int particle_distance_sort(const void* a, const void* b){
3830  if((*((Particle*)a)).aux > (*((Particle*)b)).aux)
3831    return -1;
3832  else
3833    return 1;
3834}
3835
3836/*
3837void soft_read_verts()
3838{
3839  glReadPixels(0, 0, psys->psys_width, psys->psys_height, GL_RGB, GL_FLOAT, vert);
3840  //fprintf(stderr, "soft_read_vert");
3841
3842  //cpu sort the distance 
3843  Particle* p = (Particle*) malloc(sizeof(Particle)* psys->psys_width * psys->psys_height);
3844  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
3845    float x = vert[3*i];
3846    float y = vert[3*i+1];
3847    float z = vert[3*i+2];
3848
3849    float dis = (x-live_obj_x)*(x-live_obj_x) + (y-live_obj_y)*(y-live_obj_y) + (z-live_obj_z)*(z-live_obj_z);
3850    p[i].x = x;
3851    p[i].y = y;
3852    p[i].z = z;
3853    p[i].aux = dis;
3854  }
3855
3856  qsort(p, psys->psys_width * psys->psys_height, sizeof(Particle), particle_distance_sort);
3857
3858  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
3859    vert[3*i] = p[i].x;
3860    vert[3*i+1] = p[i].y;
3861    vert[3*i+2] = p[i].z;
3862  }
3863
3864  free(p);
3865}
3866*/
3867
3868/*
3869//display the content of a texture as a screen aligned quad
3870void display_texture(NVISid tex, int width, int height)
3871{
3872   glPushMatrix();
3873
3874   glEnable(GL_TEXTURE_2D);
3875   glBindTexture(GL_TEXTURE_RECTANGLE_NV, tex);
3876
3877   glViewport(0, 0, width, height);
3878   glMatrixMode(GL_PROJECTION);
3879   glLoadIdentity();
3880   gluOrtho2D(0, width, 0, height);
3881   glMatrixMode(GL_MODELVIEW);
3882   glLoadIdentity();
3883
3884   cgGLBindProgram(m_passthru_fprog);
3885   cgGLEnableProfile(CG_PROFILE_FP30);
3886
3887   cgGLSetParameter4f(m_passthru_scale_param, 1.0, 1.0, 1.0, 1.0);
3888   cgGLSetParameter4f(m_passthru_bias_param, 0.0, 0.0, 0.0, 0.0);
3889   
3890   draw_quad(width, height, width, height);
3891   cgGLDisableProfile(CG_PROFILE_FP30);
3892
3893   glPopMatrix();
3894
3895   assert(glGetError()==0);
3896}
3897*/
3898
3899
3900//draw vertices in the main memory
3901/*
3902void soft_display_verts()
3903{
3904  glDisable(GL_TEXTURE_2D);
3905  glEnable(GL_BLEND);
3906
3907  glPointSize(0.5);
3908  glColor4f(0,0.8,0.8,0.6);
3909  glBegin(GL_POINTS);
3910  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
3911    glVertex3f(vert[3*i], vert[3*i+1], vert[3*i+2]);
3912  }
3913  glEnd();
3914  //fprintf(stderr, "soft_display_vert");
3915}
3916*/
3917
3918#if 0
3919
3920//oddeven sort on GPU
3921void sortstep()
3922{
3923    // perform one step of the current sorting algorithm
3924
3925        /*
3926    // swap buffers
3927    int sourceBuffer = targetBuffer;
3928    targetBuffer = (targetBuffer+1)%2;   
3929    int pstage = (1<<stage);
3930    int ppass  = (1<<pass);
3931    int activeBitonicShader = 0;
3932
3933#ifdef _WIN32
3934    buffer->BindBuffer(wglTargets[sourceBuffer]);
3935#else
3936    buffer->BindBuffer(glTargets[sourceBuffer]);
3937#endif
3938    if (buffer->IsDoubleBuffered()) glDrawBuffer(glTargets[targetBuffer]);
3939    */
3940
3941    checkGLError("after db");
3942
3943    int pstage = (1<<stage);
3944    int ppass  = (1<<pass);
3945    //int activeBitonicShader = 0;
3946
3947    // switch on correct sorting shader
3948    oddevenMergeSort.bind();
3949    glUniform3fARB(oddevenMergeSort.getUniformLocation("Param1"), float(pstage+pstage),
3950                   float(ppass%pstage), float((pstage+pstage)-(ppass%pstage)-1));
3951    glUniform3fARB(oddevenMergeSort.getUniformLocation("Param2"), float(psys_width), float(psys_height), float(ppass));
3952    glUniform1iARB(oddevenMergeSort.getUniformLocation("Data"), 0);
3953    staticdebugmsg("sort","stage "<<pstage<<" pass "<<ppass);
3954
3955    // This clear is not necessary for sort to function. But if we are in interactive mode
3956    // unused parts of the texture that are visible will look bad.
3957    //if (!perfTest) glClear(GL_COLOR_BUFFER_BIT);
3958
3959    //buffer->Bind();
3960    //buffer->EnableTextureTarget();
3961
3962    // initiate the sorting step on the GPU
3963    // a full-screen quad
3964    glBegin(GL_QUADS);
3965    /*
3966    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,0.0f,0.0f,1.0f); glVertex2f(-1.0f,-1.0f);
3967    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),0.0f,1.0f,1.0f); glVertex2f(1.0f,-1.0f);
3968    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),float(psys_height),1.0f,0.0f); glVertex2f(1.0f,1.0f);       
3969    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,float(psys_height),0.0f,0.0f); glVertex2f(-1.0f,1.0f);   
3970    */
3971    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,0.0f,0.0f,1.0f); glVertex2f(0.,0.);       
3972    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),0.0f,1.0f,1.0f); glVertex2f(float(psys_width), 0.);
3973    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),float(psys_height),1.0f,0.0f); glVertex2f(float(psys_width), float(psys_height));   
3974    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,float(psys_height),0.0f,0.0f); glVertex2f(0., float(psys_height));       
3975    glEnd();
3976
3977
3978    // switch off sorting shader
3979    oddevenMergeSort.release();
3980
3981    //buffer->DisableTextureTarget();
3982
3983    assert(glGetError()==0);
3984}
3985
3986#endif
3987
3988
3989void draw_3d_axis()
3990{
3991    glDisable(GL_TEXTURE_2D);
3992    glEnable(GL_DEPTH_TEST);
3993
3994        //draw axes
3995        GLUquadric *obj;
3996        obj = gluNewQuadric();
3997       
3998        glDepthFunc(GL_LESS);
3999        glEnable(GL_DEPTH_TEST);
4000        glDisable(GL_BLEND);
4001
4002        int segments = 50;
4003
4004        glColor3f(0.8, 0.8, 0.8);
4005        glPushMatrix();
4006        glTranslatef(0.4, 0., 0.);
4007        glRotatef(90, 1, 0, 0);
4008        glRotatef(180, 0, 1, 0);
4009        glScalef(0.0005, 0.0005, 0.0005);
4010        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'x');
4011        glPopMatrix(); 
4012
4013        glPushMatrix();
4014        glTranslatef(0., 0.4, 0.);
4015        glRotatef(90, 1, 0, 0);
4016        glRotatef(180, 0, 1, 0);
4017        glScalef(0.0005, 0.0005, 0.0005);
4018        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'y');
4019        glPopMatrix(); 
4020
4021        glPushMatrix();
4022        glTranslatef(0., 0., 0.4);
4023        glRotatef(90, 1, 0, 0);
4024        glRotatef(180, 0, 1, 0);
4025        glScalef(0.0005, 0.0005, 0.0005);
4026        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'z');
4027        glPopMatrix(); 
4028
4029        glEnable(GL_LIGHTING);
4030        glEnable(GL_LIGHT0);
4031
4032        glColor3f(0.2, 0.2, 0.8);
4033        glPushMatrix();
4034        glutSolidSphere(0.02, segments, segments );
4035        glPopMatrix();
4036
4037        glPushMatrix();
4038        glRotatef(-90, 1, 0, 0);       
4039        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
4040        glPopMatrix(); 
4041
4042        glPushMatrix();
4043        glTranslatef(0., 0.3, 0.);
4044        glRotatef(-90, 1, 0, 0);       
4045        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
4046        glPopMatrix(); 
4047
4048        glPushMatrix();
4049        glRotatef(90, 0, 1, 0);
4050        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
4051        glPopMatrix(); 
4052
4053        glPushMatrix();
4054        glTranslatef(0.3, 0., 0.);
4055        glRotatef(90, 0, 1, 0);
4056        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
4057        glPopMatrix(); 
4058
4059        glPushMatrix();
4060        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
4061        glPopMatrix(); 
4062
4063        glPushMatrix();
4064        glTranslatef(0., 0., 0.3);
4065        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
4066        glPopMatrix(); 
4067
4068        glDisable(GL_LIGHTING);
4069        glDisable(GL_DEPTH_TEST);
4070        gluDeleteQuadric(obj);
4071
4072    glEnable(GL_TEXTURE_2D);
4073    glDisable(GL_DEPTH_TEST);
4074}
4075
4076/*
4077void draw_axis()
4078{
4079  glDisable(GL_TEXTURE_2D);
4080  glEnable(GL_DEPTH_TEST);
4081
4082  //red x
4083  glColor3f(1,0,0);
4084  glBegin(GL_LINES);
4085    glVertex3f(0,0,0);
4086    glVertex3f(1.5,0,0);
4087  glEnd();
4088
4089  //blue y
4090  glColor3f(0,0,1);
4091  glBegin(GL_LINES);
4092    glVertex3f(0,0,0);
4093    glVertex3f(0,1.5,0);
4094  glEnd();
4095
4096  //green z
4097  glColor3f(0,1,0);
4098  glBegin(GL_LINES);
4099    glVertex3f(0,0,0);
4100    glVertex3f(0,0,1.5);
4101  glEnd();
4102
4103  glEnable(GL_TEXTURE_2D);
4104  glDisable(GL_DEPTH_TEST);
4105}
4106*/
4107
4108
4109
4110/*----------------------------------------------------*/
4111void display()
4112{
4113    assert(glGetError()==0);
4114
4115    //lic->convolve(); //flow line integral convolution
4116    //psys->advect(); //advect particles
4117
4118    //start final rendering
4119    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear screen
4120
4121    if (volume_mode)
4122    {
4123        //3D rendering mode
4124        glEnable(GL_TEXTURE_2D);
4125        glEnable(GL_DEPTH_TEST);
4126
4127        //camera setting activated
4128        cam->activate();
4129
4130        //set up the orientation of items in the scene.
4131        glPushMatrix();
4132        switch (updir) {
4133        case 1:  // x
4134            glRotatef(90, 0, 0, 1);
4135            glRotatef(90, 1, 0, 0);
4136            break;
4137
4138        case 2:  // y
4139            // this is the default
4140            break;
4141
4142        case 3:  // z
4143            glRotatef(-90, 1, 0, 0);
4144            glRotatef(-90, 0, 0, 1);
4145            break;
4146
4147        case -1:  // -x
4148            glRotatef(-90, 0, 0, 1);
4149            break;
4150
4151        case -2:  // -y
4152            glRotatef(180, 0, 0, 1);
4153            glRotatef(-90, 0, 1, 0);
4154            break;
4155
4156        case -3:  // -z
4157            glRotatef(90, 1, 0, 0);
4158            break;
4159        }
4160
4161        glPushMatrix();
4162        //now render things in the scene
4163        if (axis_on)
4164        {
4165            draw_3d_axis();
4166        }
4167
4168        if (g_grid->isVisible())
4169        {
4170            g_grid->render();
4171        }
4172
4173        //lic->render();        //display the line integral convolution result
4174        //soft_display_verts();
4175        //perf->enable();
4176        //psys->render();
4177        //perf->disable();
4178        //fprintf(stderr, "particle pixels: %d\n", perf->get_pixel_count());
4179        //perf->reset();
4180
4181        perf->enable();
4182        g_vol_render->render_all();
4183
4184        perf->disable();
4185
4186        for (int ii = 0; ii < g_heightMap.size(); ++ii)
4187        {
4188            if (g_heightMap[ii]->isVisible())
4189                g_heightMap[ii]->render();
4190        }
4191
4192        glPopMatrix();
4193
4194        float mat[16];
4195        glGetFloatv(GL_MODELVIEW_MATRIX, mat);
4196
4197        for (int i = 0; i < g_pointSet.size(); ++i)
4198        {
4199            if (g_pointSet[i]->isVisible())
4200            {
4201                g_pointset_renderer->render(g_pointSet[i]->getCluster(), mat, g_pointSet[i]->getSortLevel());
4202            }
4203        }
4204
4205        glPopMatrix();
4206   }
4207   else {
4208        //2D rendering mode
4209        perf->enable();
4210        plane_render->render();
4211        perf->disable();
4212   }
4213
4214#ifdef XINETD
4215    float cost  = perf->get_pixel_count();
4216    write(3, &cost, sizeof(cost));
4217#endif
4218    perf->reset();
4219
4220}
4221
4222
4223void mouse(int button, int state, int x, int y){
4224  if(button==GLUT_LEFT_BUTTON){
4225
4226    if(state==GLUT_DOWN){
4227      left_last_x = x;
4228      left_last_y = y;
4229      left_down = true;
4230      right_down = false;
4231    }
4232    else{
4233      left_down = false;
4234      right_down = false;
4235    }
4236  }
4237  else{
4238    //fprintf(stderr, "right mouse\n");
4239
4240    if(state==GLUT_DOWN){
4241      //fprintf(stderr, "right mouse down\n");
4242      right_last_x = x;
4243      right_last_y = y;
4244      left_down = false;
4245      right_down = true;
4246    }
4247    else{
4248      //fprintf(stderr, "right mouse up\n");
4249      left_down = false;
4250      right_down = false;
4251    }
4252  }
4253}
4254
4255
4256void update_rot(int delta_x, int delta_y){
4257        live_rot_x += delta_x;
4258        live_rot_y += delta_y;
4259
4260        if(live_rot_x > 360.0)
4261                live_rot_x -= 360.0;   
4262        else if(live_rot_x < -360.0)
4263                live_rot_x += 360.0;
4264
4265        if(live_rot_y > 360.0)
4266                live_rot_y -= 360.0;   
4267        else if(live_rot_y < -360.0)
4268                live_rot_y += 360.0;
4269
4270        cam->rotate(live_rot_x, live_rot_y, live_rot_z);
4271}
4272
4273
4274void update_trans(int delta_x, int delta_y, int delta_z){
4275        live_obj_x += delta_x*0.03;
4276        live_obj_y += delta_y*0.03;
4277        live_obj_z += delta_z*0.03;
4278}
4279
4280void end_service();
4281
4282void keyboard(unsigned char key, int x, int y){
4283 
4284   bool log = false;
4285
4286   switch (key){
4287        case 'q':
4288#ifdef XINETD
4289        //end_service();
4290        NvExitService();
4291#endif
4292                exit(0);
4293                break;
4294        case '+':
4295                lic_slice_z+=0.05;
4296                lic->set_offset(lic_slice_z);
4297                break;
4298        case '-':
4299                lic_slice_z-=0.05;
4300                lic->set_offset(lic_slice_z);
4301                break;
4302        case ',':
4303                lic_slice_x+=0.05;
4304                //init_particles();
4305                break;
4306        case '.':
4307                lic_slice_x-=0.05;
4308                //init_particles();
4309                break;
4310        case '1':
4311                //advect = true;
4312                break;
4313        case '2':
4314                //psys_x+=0.05;
4315                break;
4316        case '3':
4317                //psys_x-=0.05;
4318                break;
4319        case 'w': //zoom out
4320                live_obj_z-=0.05;
4321                log = true;
4322                cam->move(live_obj_x, live_obj_y, live_obj_z);
4323                break;
4324        case 's': //zoom in
4325                live_obj_z+=0.05;
4326                log = true;
4327                cam->move(live_obj_x, live_obj_y, live_obj_z);
4328                break;
4329        case 'a': //left
4330                live_obj_x-=0.05;
4331                log = true;
4332                cam->move(live_obj_x, live_obj_y, live_obj_z);
4333                break;
4334        case 'd': //right
4335                live_obj_x+=0.05;
4336                log = true;
4337                cam->move(live_obj_x, live_obj_y, live_obj_z);
4338                break;
4339        case 'i':
4340                //init_particles();
4341                break;
4342        case 'v':
4343                g_vol_render->switch_volume_mode();
4344                break;
4345        case 'b':
4346                g_vol_render->switch_slice_mode();
4347                break;
4348        case 'n':
4349                resize_offscreen_buffer(win_width*2, win_height*2);
4350                break;
4351        case 'm':
4352                resize_offscreen_buffer(win_width/2, win_height/2);
4353                break;
4354
4355        default:
4356                break;
4357    }   
4358
4359#ifdef EVENTLOG
4360   if(log){
4361     float param[3] = {live_obj_x, live_obj_y, live_obj_z};
4362     Event* tmp = new Event(EVENT_MOVE, param, NvGetTimeInterval());
4363     tmp->write(event_log);
4364     delete tmp;
4365   }
4366#endif
4367}
4368
4369void motion(int x, int y)
4370{
4371
4372    int old_x, old_y;   
4373
4374    if(left_down){
4375      old_x = left_last_x;
4376      old_y = left_last_y;   
4377    }
4378    else if(right_down){
4379      old_x = right_last_x;
4380      old_y = right_last_y;   
4381    }
4382
4383    int delta_x = x - old_x;
4384    int delta_y = y - old_y;
4385
4386    //more coarse event handling
4387    //if(abs(delta_x)<10 && abs(delta_y)<10)
4388      //return;
4389
4390    if(left_down){
4391      left_last_x = x;
4392      left_last_y = y;
4393
4394      update_rot(-delta_y, -delta_x);
4395    }
4396    else if (right_down){
4397      //fprintf(stderr, "right mouse motion (%d,%d)\n", x, y);
4398
4399      right_last_x = x;
4400      right_last_y = y;
4401
4402      update_trans(0, 0, delta_x);
4403    }
4404
4405#ifdef EVENTLOG
4406    float param[3] = {live_rot_x, live_rot_y, live_rot_z};
4407    Event* tmp = new Event(EVENT_ROTATE, param, NvGetTimeInterval());
4408    tmp->write(event_log);
4409    delete tmp;
4410#endif
4411
4412    glutPostRedisplay();
4413}
4414
4415/*
4416#ifdef XINETD
4417void init_service(){
4418  //open log and map stderr to log file
4419  xinetd_log = fopen("/tmp/log.txt", "w");
4420  close(2);
4421  dup2(fileno(xinetd_log), 2);
4422  dup2(2,1);
4423
4424  //flush junk
4425  fflush(stdout);
4426  fflush(stderr);
4427}
4428
4429void end_service(){
4430  //close log file
4431  fclose(xinetd_log);
4432}
4433#endif
4434
4435void init_event_log(){
4436  event_log = fopen("event.txt", "w");
4437  assert(event_log!=0);
4438
4439  struct timeval time;
4440  gettimeofday(&time, NULL);
4441  cur_time = time.tv_sec*1000. + time.tv_usec/1000.;
4442}
4443
4444void end_event_log(){
4445  fclose(event_log);
4446}
4447
4448double get_time_interval(){
4449  struct timeval time;
4450  gettimeofday(&time, NULL);
4451  double new_time = time.tv_sec*1000. + time.tv_usec/1000.;
4452
4453  double interval = new_time - cur_time;
4454  cur_time = new_time;
4455  return interval;
4456}
4457*/
4458
4459void removeAllData()
4460{
4461    //
4462}
4463
4464
4465/*----------------------------------------------------*/
4466int main(int argc, char** argv)
4467{
4468
4469    char path [1000];
4470    while(1) {
4471        int c;
4472        int this_option_optind = optind ? optind : 1;
4473        int option_index = 0;
4474        struct option long_options[] = {
4475          // name, has_arg, flag, val
4476          { 0,0,0,0 },
4477        };
4478
4479        c = getopt_long(argc, argv, "+p:", long_options, &option_index);
4480        if (c == -1)
4481            break;
4482
4483        switch(c) {
4484            case 'p':
4485                strncpy(path, optarg, sizeof(path));
4486                break;
4487            default:
4488                fprintf(stderr,"Don't know what option '%c'.\n", c);
4489                return 1;
4490        }
4491    }
4492
4493    R2FilePath::getInstance()->setWorkingDirectory(argc, (const char**) argv);
4494
4495#ifdef XINETD
4496   signal(SIGPIPE,SIG_IGN);
4497   NvInitService();
4498#endif
4499
4500   glutInit(&argc, argv);
4501   glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
4502
4503   glutInitWindowSize(win_width, win_height);
4504
4505   glutInitWindowPosition(10, 10);
4506   render_window = glutCreateWindow(argv[0]);
4507   glutDisplayFunc(display);
4508
4509#ifndef XINETD
4510   glutMouseFunc(mouse);
4511   glutMotionFunc(motion);
4512   glutKeyboardFunc(keyboard);
4513#endif
4514
4515   glutIdleFunc(idle);
4516   glutReshapeFunc(resize_offscreen_buffer);
4517
4518   NvInit(path);
4519   initGL();
4520   initTcl();
4521
4522#ifdef EVENTLOG
4523   NvInitEventLog();
4524#endif
4525    //event loop
4526    glutMainLoop();
4527
4528    removeAllData();
4529
4530    NvExit();
4531
4532    return 0;
4533}
Note: See TracBrowser for help on using the repository browser.