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

Last change on this file since 827 was 827, checked in by vrinside, 16 years ago

Added color blending with point sprite texture

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