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

Last change on this file since 829 was 829, checked in by vrinside, 13 years ago

Moved all the functions related to TCL to Command.cpp

File size: 75.2 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
58// R2 headers
59#include <R2/R2FilePath.h>
60#include <R2/R2Fonts.h>
61
62//render server
63
64extern VolumeRenderer* g_vol_render;
65extern PointSetRenderer* g_pointset_renderer;
66extern NvColorTableRenderer* g_color_table_renderer;
67extern Grid* g_grid;
68PlaneRenderer* plane_render;
69Camera* cam;
70
71//if true nanovis renders volumes in 3D, if not renders 2D plane
72bool volume_mode = true;
73
74// color table for built-in transfer function editor
75float color_table[256][4];     
76
77// default transfer function
78char *def_transfunc = "transfunc define default {\n\
79  0.0  1 1 1\n\
80  0.2  1 1 0\n\
81  0.4  0 1 0\n\
82  0.6  0 1 1\n\
83  0.8  0 0 1\n\
84  1.0  1 0 1\n\
85} {\n\
86  0.00  1.0\n\
87  0.05  0.0\n\
88  0.15  0.0\n\
89  0.20  1.0\n\
90  0.25  0.0\n\
91  0.35  0.0\n\
92  0.40  1.0\n\
93  0.45  0.0\n\
94  0.55  0.0\n\
95  0.60  1.0\n\
96  0.65  0.0\n\
97  0.75  0.0\n\
98  0.80  1.0\n\
99  0.85  0.0\n\
100  0.95  0.0\n\
101  1.00  1.0\n\
102}";
103
104/*
105#ifdef XINETD
106FILE* xinetd_log;
107#endif
108
109FILE* event_log;
110//log
111void init_event_log();
112void end_event_log();
113double cur_time;        //in seconds
114double get_time_interval();
115*/
116
117int render_window;              //the handle of the render window;
118bool axis_on = true;
119
120// in Command.cpp
121extern void xinetd_listen();
122extern void initTcl();
123
124
125// forward declarations
126//void init_particles();
127void get_slice_vectors();
128Rappture::Outcome load_volume_stream(int index, std::iostream& fin);
129Rappture::Outcome load_volume_stream2(int index, std::iostream& fin);
130void load_volume(int index, int width, int height, int depth, int n_component, float* data, double vmin, double vmax,
131                double nzero_min);
132TransferFunction* get_transfunc(char *name);
133void resize_offscreen_buffer(int w, int h);
134void offscreen_buffer_capture();
135void bmp_header_add_int(unsigned char* header, int& pos, int data);
136void bmp_write(const char* cmd);
137void bmp_write_to_file();
138void display();
139void display_offscreen_buffer();
140void read_screen();
141
142//ParticleSystem* psys;
143//float psys_x=0.4, psys_y=0, psys_z=0;
144
145Lic* lic;
146
147//frame buffer for final rendering
148unsigned char* screen_buffer = NULL;
149NVISid final_fbo, final_color_tex, final_depth_rb;
150
151//bool advect=false;
152float vert[NMESH*NMESH*3];              //particle positions in main memory
153float slice_vector[NMESH*NMESH*4];      //per slice vectors in main memory
154
155int n_volumes = 0;
156// pointers to volumes, currently handle up to 10 volumes
157vector<Volume*> volume;
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// FOR NANOVIS.H
182//variables for mouse events
183float live_rot_x = 90.;         //object rotation angles
184float live_rot_y = 180.;
185float live_rot_z = -135;
186
187float live_obj_x = -0.0;        //object translation location from the origin
188float live_obj_y = -0.0;
189float live_obj_z = -2.5;
190
191float live_diffuse = 1.;
192float live_specular = 3.;
193
194int left_last_x, left_last_y, right_last_x, right_last_y;       //last locations mouse events
195bool left_down = false;                                         
196bool right_down = false;
197
198float lic_slice_x=0, lic_slice_y=0, lic_slice_z=0.3;//image based flow visualization slice location
199
200int win_width = NPIX;                   //size of the render window
201int win_height = NPIX;                  //size of the render window
202
203
204//image based flow visualization variables
205int    iframe = 0;
206int    Npat   = 64;
207int    alpha  = (int)round(0.12*255);
208float  sa;
209float  tmax   = NPIX/(SCALE*NPN);
210float  dmax   = SCALE/NPIX;
211
212
213//currently active shader, default renders one volume only
214int cur_shader = 0;
215
216/*
217CGprogram m_copy_texcoord_fprog;
218CGprogram m_one_volume_fprog;
219CGparameter m_vol_one_volume_param;
220CGparameter m_tf_one_volume_param;
221CGparameter m_mvi_one_volume_param;
222CGparameter m_mv_one_volume_param;
223CGparameter m_render_param_one_volume_param;
224*/
225
226/*
227CGprogram m_vert_std_vprog;
228CGparameter m_mvp_vert_std_param;
229CGparameter m_mvi_vert_std_param;
230*/
231
232using namespace std;
233
234#define RM_VOLUME 1
235#define RM_POINT 2
236
237int renderMode = RM_VOLUME;
238
239
240
241/*
242 * ----------------------------------------------------------------------
243 * USAGE: debug("string", ...)
244 *
245 * Use this anywhere within the package to send a debug message
246 * back to the client.  The string can have % fields as used by
247 * the printf() package.  Any remaining arguments are treated as
248 * field substitutions on that.
249 * ----------------------------------------------------------------------
250 */
251void
252debug(char *str)
253{
254    write(0, str, strlen(str));
255}
256
257void
258debug(char *str, double v1)
259{
260    char buffer[512];
261    sprintf(buffer, str, v1);
262    write(0, buffer, strlen(buffer));
263}
264
265void
266debug(char *str, double v1, double v2)
267{
268    char buffer[512];
269    sprintf(buffer, str, v1, v2);
270    write(0, buffer, strlen(buffer));
271}
272
273void
274debug(char *str, double v1, double v2, double v3)
275{
276    char buffer[512];
277    sprintf(buffer, str, v1, v2, v3);
278    write(0, buffer, strlen(buffer));
279}
280
281//report errors related to CG shaders
282void cgErrorCallback(void)
283{
284    CGerror lastError = cgGetError();
285    if(lastError) {
286        const char *listing = cgGetLastListing(g_context);
287        printf("\n---------------------------------------------------\n");
288        printf("%s\n\n", cgGetErrorString(lastError));
289        printf("%s\n", listing);
290        printf("-----------------------------------------------------\n");
291        printf("Cg error, exiting...\n");
292        cgDestroyContext(g_context);
293        exit(-1);
294    }
295}
296
297/* Load a 3D vector field from a dx-format file
298 */
299void
300load_vector_stream(int index, std::iostream& fin) {
301    int dummy, nx, ny, nz, nxy, npts;
302    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
303    char line[128], type[128], *start;
304
305    do {
306        fin.getline(line,sizeof(line)-1);
307        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
308            ;  // skip leading blanks
309
310        if (*start != '#') {  // skip comment lines
311            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
312                // found grid size
313            }
314            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
315                // found origin
316            }
317            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
318                // found one of the delta lines
319                if (ddx != 0.0) { dx = ddx; }
320                else if (ddy != 0.0) { dy = ddy; }
321                else if (ddz != 0.0) { dz = ddz; }
322            }
323            else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
324                if (npts != nx*ny*nz) {
325                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
326                    return;
327                }
328                break;
329            }
330            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
331                if (npts != nx*ny*nz) {
332                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
333                    return;
334                }
335                break;
336            }
337        }
338    } while (!fin.eof());
339
340    // read data points
341    if (!fin.eof()) {
342        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
343        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
344        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
345        Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
346        Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
347        Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
348
349        double vx, vy, vz;
350        int nread = 0;
351        for (int ix=0; ix < nx; ix++) {
352            for (int iy=0; iy < ny; iy++) {
353                for (int iz=0; iz < nz; iz++) {
354                    if (fin.eof() || nread > npts) {
355                        break;
356                    }
357                    fin.getline(line,sizeof(line)-1);
358                    if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
359                        int nindex = iz*nx*ny + iy*nx + ix;
360                        xfield.define(nindex, vx);
361                        yfield.define(nindex, vy);
362                        zfield.define(nindex, vz);
363                        nread++;
364                    }
365                }
366            }
367        }
368
369        // make sure that we read all of the expected points
370        if (nread != nx*ny*nz) {
371            std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << nread << " points" << std::endl;
372            return;
373        }
374
375        // figure out a good mesh spacing
376        int nsample = 30;
377        dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
378        dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
379        dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
380        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
381
382        nx = (int)ceil(dx/dmin);
383        ny = (int)ceil(dy/dmin);
384        nz = (int)ceil(dz/dmin);
385
386#ifndef NV40
387        // must be an even power of 2 for older cards
388            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
389            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
390            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
391#endif
392
393        float *data = new float[3*nx*ny*nz];
394
395        std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl;
396
397        // generate the uniformly sampled data that we need for a volume
398        double vmin = 0.0;
399        double vmax = 0.0;
400        double nzero_min = 0.0;
401        int ngen = 0;
402        for (int iz=0; iz < nz; iz++) {
403            double zval = z0 + iz*dmin;
404            for (int iy=0; iy < ny; iy++) {
405                double yval = y0 + iy*dmin;
406                for (int ix=0; ix < nx; ix++) {
407                    double xval = x0 + ix*dmin;
408
409                    vx = xfield.value(xval,yval,zval);
410                    vy = yfield.value(xval,yval,zval);
411                    vz = zfield.value(xval,yval,zval);
412                    //vx = 1;
413                    //vy = 1;
414                    vz = 0;
415                    double vm = sqrt(vx*vx + vy*vy + vz*vz);
416
417                    if (vm < vmin) { vmin = vm; }
418                    if (vm > vmax) { vmax = vm; }
419                    if (vm != 0.0f && vm < nzero_min)
420                    {
421                        nzero_min = vm;
422                    }
423
424                    data[ngen++] = vx;
425                    data[ngen++] = vy;
426                    data[ngen++] = vz;
427                }
428            }
429        }
430
431        ngen = 0;
432        for (ngen=0; ngen < npts; ngen++) {
433            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
434        }
435
436        load_volume(index, nx, ny, nz, 3, data, vmin, vmax, nzero_min);
437        delete [] data;
438    } else {
439        std::cerr << "WARNING: data not found in stream" << std::endl;
440    }
441}
442
443
444/* Load a 3D volume from a dx-format file
445 */
446Rappture::Outcome
447load_volume_stream2(int index, std::iostream& fin) {
448    Rappture::Outcome result;
449
450    Rappture::MeshTri2D xymesh;
451    int dummy, nx, ny, nz, nxy, npts;
452    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
453    char line[128], type[128], *start;
454
455    int isrect = 1;
456
457    float* voldata = 0;
458    do {
459        fin.getline(line,sizeof(line)-1);
460        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
461            ;  // skip leading blanks
462
463        if (*start != '#') {  // skip comment lines
464            printf("%s\n", line);
465            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
466                // found grid size
467                isrect = 1;
468            }
469            else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
470                isrect = 0;
471                double xx, yy, zz;
472                for (int i=0; i < nxy; i++) {
473                    fin.getline(line, sizeof(line));
474                    fin.getline(line,sizeof(line)-1);
475                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
476                        xymesh.addNode( Rappture::Node2D(xx,yy) );
477                    }
478                }
479                char mesg[256];
480                sprintf(mesg,"test");
481                result.error(mesg);
482                return result;
483
484                char fpts[128];
485                sprintf(fpts, "/tmp/tmppts%d", getpid());
486                char fcells[128];
487                sprintf(fcells, "/tmp/tmpcells%d", getpid());
488
489                std::ofstream ftmp(fpts);
490                // save corners of bounding box first, to work around meshing
491                // problems in voronoi utility
492                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
493                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
494                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
495                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
496                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
497                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
498                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
499                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
500                for (int i=0; i < nxy; i++) {
501                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
502               
503                }
504                ftmp.close();
505
506                char cmdstr[512];
507                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
508                if (system(cmdstr) == 0) {
509                    int cx, cy, cz;
510                    std::ifstream ftri(fcells);
511                    while (!ftri.eof()) {
512                        ftri.getline(line,sizeof(line)-1);
513                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
514                            if (cx >= 4 && cy >= 4 && cz >= 4) {
515                                // skip first 4 boundary points
516                                xymesh.addCell(cx-4, cy-4, cz-4);
517                            }
518                        }
519                    }
520                    ftri.close();
521                } else {
522                    return result.error("triangularization failed");
523                }
524
525                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
526                system(cmdstr);
527            }
528            else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
529                // found z-grid
530            }
531            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
532                // found origin
533            }
534            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
535                // found one of the delta lines
536                if (ddx != 0.0) { dx = ddx; }
537                else if (ddy != 0.0) { dy = ddy; }
538                else if (ddz != 0.0) { dz = ddz; }
539            }
540            else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
541                if (isrect && (npts != nx*ny*nz)) {
542                    char mesg[256];
543                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
544                    return result.error(mesg);
545                }
546                else if (!isrect && (npts != nxy*nz)) {
547                    char mesg[256];
548                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
549                    return result.error(mesg);
550                }
551                break;
552            }
553            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
554                if (npts != nx*ny*nz) {
555                    char mesg[256];
556                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
557                    return result.error(mesg);
558                }
559                break;
560            }
561        }
562    } while (!fin.eof());
563
564    // read data points
565    if (!fin.eof()) {
566        if (isrect) {
567            double dval[6];
568            int nread = 0;
569            int ix = 0;
570            int iy = 0;
571            int iz = 0;
572            float* data = new float[nx *  ny *  nz * 4];
573            memset(data, 0, nx*ny*nz*4);
574            double vmin = 1e21;
575            double nzero_min = 1e21;
576            double vmax = -1e21;
577
578
579            while (!fin.eof() && nread < npts) {
580                fin.getline(line,sizeof(line)-1);
581                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
582
583                for (int p=0; p < n; p++) {
584                    int nindex = (iz*nx*ny + iy*nx + ix) * 4;
585                    data[nindex] = dval[p];
586
587                    if (dval[p] < vmin) vmin = dval[p];
588                    if (dval[p] > vmax) vmax = dval[p];
589                    if (dval[p] != 0.0f && dval[p] < nzero_min)
590                    {
591                         nzero_min = dval[p];
592                    }
593
594                    nread++;
595                    if (++iz >= nz) {
596                        iz = 0;
597                        if (++iy >= ny) {
598                            iy = 0;
599                            ++ix;
600                        }
601                    }
602                }
603            }
604
605            // make sure that we read all of the expected points
606            if (nread != nx*ny*nz) {
607                char mesg[256];
608                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
609                result.error(mesg);
610                return result;
611            }
612
613            double dv = vmax - vmin;
614            int count = nx*ny*nz;
615            int ngen = 0;
616            double v;
617            printf("test2\n");
618                        fflush(stdout);
619            if (dv == 0.0) { dv = 1.0; }
620            for (int i = 0; i < count; ++i)
621            {
622                v = data[ngen];
623                // scale all values [0-1], -1 => out of bounds
624                v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
625                data[ngen] = v;
626                ngen += 4;
627            }
628
629            // Compute the gradient of this data.  BE CAREFUL: center
630            // calculation on each node to avoid skew in either direction.
631            ngen = 0;
632            for (int iz=0; iz < nz; iz++) {
633                for (int iy=0; iy < ny; iy++) {
634                    for (int ix=0; ix < nx; ix++) {
635                        // gradient in x-direction
636                        double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
637                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
638                        if (valm1 < 0 || valp1 < 0) {
639                            data[ngen+1] = 0.0;
640                        } else {
641                            data[ngen+1] = valp1-valm1; // assume dx=1
642                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
643                        }
644
645                        // gradient in y-direction
646                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
647                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
648                        if (valm1 < 0 || valp1 < 0) {
649                            data[ngen+2] = 0.0;
650                        } else {
651                            data[ngen+2] = valp1-valm1; // assume dx=1
652                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
653                        }
654
655                        // gradient in z-direction
656                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
657                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
658                        if (valm1 < 0 || valp1 < 0) {
659                            data[ngen+3] = 0.0;
660                        } else {
661                            data[ngen+3] = valp1-valm1; // assume dx=1
662                            //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
663                        }
664
665                        ngen += 4;
666                    }
667                }
668            }
669
670            dx = nx;
671            dy = ny;
672            dz = nz;
673            load_volume(index, nx, ny, nz, 4, data,
674                vmin, vmax, nzero_min);
675
676            delete [] data;
677
678        } else {
679            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
680            Rappture::FieldPrism3D field(xymesh, zgrid);
681
682            double dval;
683            int nread = 0;
684            int ixy = 0;
685            int iz = 0;
686            while (!fin.eof() && nread < npts) {
687                if (!(fin >> dval).fail()) {
688                    int nid = nxy*iz + ixy;
689                    field.define(nid, dval);
690
691                    nread++;
692                    if (++iz >= nz) {
693                        iz = 0;
694                        ixy++;
695                    }
696                }
697            }
698
699            // make sure that we read all of the expected points
700            if (nread != nxy*nz) {
701                char mesg[256];
702                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
703                return result.error(mesg);
704            }
705
706            // figure out a good mesh spacing
707            int nsample = 30;
708            x0 = field.rangeMin(Rappture::xaxis);
709            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
710            y0 = field.rangeMin(Rappture::yaxis);
711            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
712            z0 = field.rangeMin(Rappture::zaxis);
713            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
714            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
715
716            nx = (int)ceil(dx/dmin);
717            ny = (int)ceil(dy/dmin);
718            nz = (int)ceil(dz/dmin);
719#ifndef NV40
720            // must be an even power of 2 for older cards
721                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
722                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
723                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
724#endif
725            float *data = new float[4*nx*ny*nz];
726
727            double vmin = field.valueMin();
728            double dv = field.valueMax() - field.valueMin();
729            if (dv == 0.0) { dv = 1.0; }
730
731            // generate the uniformly sampled data that we need for a volume
732            int ngen = 0;
733            double nzero_min = 0.0;
734            for (iz=0; iz < nz; iz++) {
735                double zval = z0 + iz*dmin;
736                for (int iy=0; iy < ny; iy++) {
737                    double yval = y0 + iy*dmin;
738                    for (int ix=0; ix < nx; ix++) {
739                        double xval = x0 + ix*dmin;
740                        double v = field.value(xval,yval,zval);
741
742                        if (v != 0.0f && v < nzero_min)
743                        {
744                            nzero_min = v;
745                        }
746                        // scale all values [0-1], -1 => out of bounds
747                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
748                        data[ngen] = v;
749
750                        ngen += 4;
751                    }
752                }
753            }
754
755            // Compute the gradient of this data.  BE CAREFUL: center
756            // calculation on each node to avoid skew in either direction.
757            ngen = 0;
758            for (int iz=0; iz < nz; iz++) {
759                for (int iy=0; iy < ny; iy++) {
760                    for (int ix=0; ix < nx; ix++) {
761                        // gradient in x-direction
762                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
763                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
764                        if (valm1 < 0 || valp1 < 0) {
765                            data[ngen+1] = 0.0;
766                        } else {
767                            //data[ngen+1] = valp1-valm1; // assume dx=1
768                            data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
769                        }
770
771                        // gradient in y-direction
772                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
773                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
774                        if (valm1 < 0 || valp1 < 0) {
775                            data[ngen+2] = 0.0;
776                        } else {
777                            //data[ngen+2] = valp1-valm1; // assume dy=1
778                            data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
779                        }
780
781                        // gradient in z-direction
782                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
783                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
784                        if (valm1 < 0 || valp1 < 0) {
785                            data[ngen+3] = 0.0;
786                        } else {
787                            //data[ngen+3] = valp1-valm1; // assume dz=1
788                            data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
789                        }
790
791                        ngen += 4;
792                    }
793                }
794            }
795
796            load_volume(index, nx, ny, nz, 4, data,
797                field.valueMin(), field.valueMax(), nzero_min);
798
799            delete [] data;
800        }
801    } else {
802        return result.error("data not found in stream");
803    }
804
805    //
806    // Center this new volume on the origin.
807    //
808    float dx0 = -0.5;
809    float dy0 = -0.5*dy/dx;
810    float dz0 = -0.5*dz/dx;
811    volume[index]->move(Vector3(dx0, dy0, dz0));
812
813    return result;
814}
815
816Rappture::Outcome
817load_volume_stream(int index, std::iostream& fin) {
818    Rappture::Outcome result;
819
820    Rappture::MeshTri2D xymesh;
821    int dummy, nx, ny, nz, nxy, npts;
822    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
823    char line[128], type[128], *start;
824
825    int isrect = 1;
826
827    do {
828        fin.getline(line,sizeof(line)-1);
829        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
830            ;  // skip leading blanks
831
832        if (*start != '#') {  // skip comment lines
833            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
834                // found grid size
835                isrect = 1;
836            }
837            else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
838                isrect = 0;
839                double xx, yy, zz;
840                for (int i=0; i < nxy; i++) {
841                    fin.getline(line,sizeof(line)-1);
842                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
843                        xymesh.addNode( Rappture::Node2D(xx,yy) );
844                    }
845                }
846
847                char fpts[128];
848                sprintf(fpts, "/tmp/tmppts%d", getpid());
849                char fcells[128];
850                sprintf(fcells, "/tmp/tmpcells%d", getpid());
851
852                std::ofstream ftmp(fpts);
853                // save corners of bounding box first, to work around meshing
854                // problems in voronoi utility
855                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
856                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
857                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
858                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
859                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
860                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
861                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
862                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
863                for (int i=0; i < nxy; i++) {
864                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
865               
866                }
867                ftmp.close();
868
869                char cmdstr[512];
870                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
871                if (system(cmdstr) == 0) {
872                    int cx, cy, cz;
873                    std::ifstream ftri(fcells);
874                    while (!ftri.eof()) {
875                        ftri.getline(line,sizeof(line)-1);
876                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
877                            if (cx >= 4 && cy >= 4 && cz >= 4) {
878                                // skip first 4 boundary points
879                                xymesh.addCell(cx-4, cy-4, cz-4);
880                            }
881                        }
882                    }
883                    ftri.close();
884                } else {
885                    return result.error("triangularization failed");
886                }
887
888                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
889                system(cmdstr);
890            }
891            else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
892                // found z-grid
893            }
894            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
895                // found origin
896            }
897            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
898                // found one of the delta lines
899                if (ddx != 0.0) { dx = ddx; }
900                else if (ddy != 0.0) { dy = ddy; }
901                else if (ddz != 0.0) { dz = ddz; }
902            }
903            else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
904                if (isrect && (npts != nx*ny*nz)) {
905                    char mesg[256];
906                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
907                    return result.error(mesg);
908                }
909                else if (!isrect && (npts != nxy*nz)) {
910                    char mesg[256];
911                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
912                    return result.error(mesg);
913                }
914                break;
915            }
916            else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
917                if (npts != nx*ny*nz) {
918                    char mesg[256];
919                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
920                    return result.error(mesg);
921                }
922                break;
923            }
924        }
925    } while (!fin.eof());
926
927    // read data points
928    if (!fin.eof()) {
929        if (isrect) {
930            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
931            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
932            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
933            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
934
935            double dval[6];
936            int nread = 0;
937            int ix = 0;
938            int iy = 0;
939            int iz = 0;
940            while (!fin.eof() && nread < npts) {
941                fin.getline(line,sizeof(line)-1);
942                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
943
944                for (int p=0; p < n; p++) {
945                    int nindex = iz*nx*ny + iy*nx + ix;
946                    field.define(nindex, dval[p]);
947                    nread++;
948                    if (++iz >= nz) {
949                        iz = 0;
950                        if (++iy >= ny) {
951                            iy = 0;
952                            ++ix;
953                        }
954                    }
955                }
956            }
957
958            // make sure that we read all of the expected points
959            if (nread != nx*ny*nz) {
960                char mesg[256];
961                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
962                result.error(mesg);
963                return result;
964            }
965
966            // figure out a good mesh spacing
967            int nsample = 30;
968            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
969            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
970            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
971            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
972
973            nx = (int)ceil(dx/dmin);
974            ny = (int)ceil(dy/dmin);
975            nz = (int)ceil(dz/dmin);
976
977#ifndef NV40
978            // must be an even power of 2 for older cards
979                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
980                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
981                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
982#endif
983
984            float *data = new float[4*nx*ny*nz];
985
986            double vmin = field.valueMin();
987            double dv = field.valueMax() - field.valueMin();
988            if (dv == 0.0) { dv = 1.0; }
989
990            // generate the uniformly sampled data that we need for a volume
991            int ngen = 0;
992            double nzero_min = 0.0;
993            for (int iz=0; iz < nz; iz++) {
994                double zval = z0 + iz*dmin;
995                for (int iy=0; iy < ny; iy++) {
996                    double yval = y0 + iy*dmin;
997                    for (int ix=0; ix < nx; ix++) {
998                        double xval = x0 + ix*dmin;
999                        double v = field.value(xval,yval,zval);
1000
1001                        if (v != 0.0f && v < nzero_min)
1002                        {
1003                            nzero_min = v;
1004                        }
1005
1006                        // scale all values [0-1], -1 => out of bounds
1007                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1008
1009                        data[ngen] = v;
1010                        ngen += 4;
1011                    }
1012                }
1013            }
1014
1015            // Compute the gradient of this data.  BE CAREFUL: center
1016            // calculation on each node to avoid skew in either direction.
1017            ngen = 0;
1018            for (int iz=0; iz < nz; iz++) {
1019                for (int iy=0; iy < ny; iy++) {
1020                    for (int ix=0; ix < nx; ix++) {
1021                        // gradient in x-direction
1022                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
1023                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
1024                        if (valm1 < 0 || valp1 < 0) {
1025                            data[ngen+1] = 0.0;
1026                        } else {
1027                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
1028                        }
1029
1030                        // gradient in y-direction
1031                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1032                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1033                        if (valm1 < 0 || valp1 < 0) {
1034                            data[ngen+2] = 0.0;
1035                        } else {
1036                            //data[ngen+2] = valp1-valm1; // assume dy=1
1037                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
1038                        }
1039
1040                        // gradient in z-direction
1041                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1042                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1043                        if (valm1 < 0 || valp1 < 0) {
1044                            data[ngen+3] = 0.0;
1045                        } else {
1046                            //data[ngen+3] = valp1-valm1; // assume dz=1
1047                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
1048                        }
1049
1050                        ngen += 4;
1051                    }
1052                }
1053            }
1054
1055            load_volume(index, nx, ny, nz, 4, data,
1056                field.valueMin(), field.valueMax(), nzero_min);
1057
1058            delete [] data;
1059
1060        } else {
1061            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1062            Rappture::FieldPrism3D field(xymesh, zgrid);
1063
1064            double dval;
1065            int nread = 0;
1066            int ixy = 0;
1067            int iz = 0;
1068            while (!fin.eof() && nread < npts) {
1069                if (!(fin >> dval).fail()) {
1070                    int nid = nxy*iz + ixy;
1071                    field.define(nid, dval);
1072
1073                    nread++;
1074                    if (++iz >= nz) {
1075                        iz = 0;
1076                        ixy++;
1077                    }
1078                }
1079            }
1080
1081            // make sure that we read all of the expected points
1082            if (nread != nxy*nz) {
1083                char mesg[256];
1084                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
1085                return result.error(mesg);
1086            }
1087
1088            // figure out a good mesh spacing
1089            int nsample = 30;
1090            x0 = field.rangeMin(Rappture::xaxis);
1091            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1092            y0 = field.rangeMin(Rappture::yaxis);
1093            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1094            z0 = field.rangeMin(Rappture::zaxis);
1095            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1096            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1097
1098            nx = (int)ceil(dx/dmin);
1099            ny = (int)ceil(dy/dmin);
1100            nz = (int)ceil(dz/dmin);
1101#ifndef NV40
1102            // must be an even power of 2 for older cards
1103                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1104                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1105                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1106#endif
1107            float *data = new float[4*nx*ny*nz];
1108
1109            double vmin = field.valueMin();
1110            double dv = field.valueMax() - field.valueMin();
1111            if (dv == 0.0) { dv = 1.0; }
1112
1113            // generate the uniformly sampled data that we need for a volume
1114            int ngen = 0;
1115            double nzero_min = 0.0;
1116            for (iz=0; iz < nz; iz++) {
1117                double zval = z0 + iz*dmin;
1118                for (int iy=0; iy < ny; iy++) {
1119                    double yval = y0 + iy*dmin;
1120                    for (int ix=0; ix < nx; ix++) {
1121                        double xval = x0 + ix*dmin;
1122                        double v = field.value(xval,yval,zval);
1123
1124                        if (v != 0.0f && v < nzero_min)
1125                        {
1126                            nzero_min = v;
1127                        }
1128                        // scale all values [0-1], -1 => out of bounds
1129                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1130                        data[ngen] = v;
1131
1132                        ngen += 4;
1133                    }
1134                }
1135            }
1136
1137            // Compute the gradient of this data.  BE CAREFUL: center
1138            // calculation on each node to avoid skew in either direction.
1139            ngen = 0;
1140            for (int iz=0; iz < nz; iz++) {
1141                for (int iy=0; iy < ny; iy++) {
1142                    for (int ix=0; ix < nx; ix++) {
1143                        // gradient in x-direction
1144                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1145                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1146                        if (valm1 < 0 || valp1 < 0) {
1147                            data[ngen+1] = 0.0;
1148                        } else {
1149                            //data[ngen+1] = valp1-valm1; // assume dx=1
1150                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
1151                        }
1152
1153                        // gradient in y-direction
1154                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1155                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1156                        if (valm1 < 0 || valp1 < 0) {
1157                            data[ngen+2] = 0.0;
1158                        } else {
1159                            //data[ngen+2] = valp1-valm1; // assume dy=1
1160                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
1161                        }
1162
1163                        // gradient in z-direction
1164                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1165                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1166                        if (valm1 < 0 || valp1 < 0) {
1167                            data[ngen+3] = 0.0;
1168                        } else {
1169                            //data[ngen+3] = valp1-valm1; // assume dz=1
1170                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
1171                        }
1172
1173                        ngen += 4;
1174                    }
1175                }
1176            }
1177
1178            load_volume(index, nx, ny, nz, 4, data,
1179                field.valueMin(), field.valueMax(), nzero_min);
1180
1181            delete [] data;
1182        }
1183    } else {
1184        return result.error("data not found in stream");
1185    }
1186
1187    //
1188    // Center this new volume on the origin.
1189    //
1190    float dx0 = -0.5;
1191    float dy0 = -0.5*dy/dx;
1192    float dz0 = -0.5*dz/dx;
1193    volume[index]->move(Vector3(dx0, dy0, dz0));
1194
1195    return result;
1196}
1197
1198/* Load a 3D volume
1199 * index: the index into the volume array, which stores pointers to 3D volume instances
1200 * data: pointer to an array of floats. 
1201 * n_component: the number of scalars for each space point.
1202 *              All component scalars for a point are placed consequtively in data array
1203 * width, height and depth: number of points in each dimension
1204 */
1205void load_volume(int index, int width, int height, int depth,
1206    int n_component, float* data, double vmin, double vmax, double nzero_min)
1207{
1208    while (n_volumes <= index) {
1209        volume.push_back(NULL);
1210        n_volumes++;
1211    }
1212
1213    if (volume[index] != NULL) {
1214        delete volume[index];
1215        volume[index] = NULL;
1216    }
1217
1218    volume[index] = new Volume(0.f, 0.f, 0.f, width, height, depth, 1.,
1219                                 n_component, data, vmin, vmax, nzero_min);
1220    assert(volume[index]!=0);
1221}
1222
1223// get a colormap 1D texture by name
1224TransferFunction*
1225get_transfunc(char *name) {
1226    Tcl_HashEntry *entryPtr;
1227    entryPtr = Tcl_FindHashEntry(&tftable, name);
1228    if (entryPtr) {
1229        return (TransferFunction*)Tcl_GetHashValue(entryPtr);
1230    }
1231    return NULL;
1232}
1233
1234
1235//Update the transfer function using local GUI in the non-server mode
1236void update_tf_texture()
1237{
1238  glutSetWindow(render_window);
1239
1240  //fprintf(stderr, "tf update\n");
1241  TransferFunction *tf = get_transfunc("default");
1242  if (tf == NULL) return;
1243
1244  float data[256*4];
1245  for(int i=0; i<256; i++)
1246  {
1247    data[4*i+0] = color_table[i][0];
1248    data[4*i+1] = color_table[i][1];
1249    data[4*i+2] = color_table[i][2];
1250    data[4*i+3] = color_table[i][3];
1251    //fprintf(stderr, "(%f,%f,%f,%f) ", data[4*i+0], data[4*i+1], data[4*i+2], data[4*i+3]);
1252  }
1253
1254  tf->update(data);
1255
1256#ifdef EVENTLOG
1257  float param[3] = {0,0,0};
1258  Event* tmp = new Event(EVENT_ROTATE, param, NvGetTimeInterval());
1259  tmp->write(event_log);
1260  delete tmp;
1261#endif
1262}
1263
1264
1265int renderLegend(int ivol, int width, int height, const char* volArg)
1266{
1267    TransferFunction *tf = NULL;
1268
1269    if (ivol < n_volumes) {
1270        tf = g_vol_render->get_volume_shading(volume[ivol]);
1271    }
1272
1273    if (tf == NULL) {
1274        return TCL_ERROR;
1275    }
1276
1277    int old_width = win_width;
1278    int old_height = win_height;
1279
1280    plane_render->set_screen_size(width, height);
1281    resize_offscreen_buffer(width, height);
1282
1283    // generate data for the legend
1284    float data[512];
1285    for (int i=0; i < 256; i++) {
1286        data[i] = data[i+256] = (float)(i/255.0);
1287    }
1288    plane[0] = new Texture2D(256, 2, GL_FLOAT, GL_LINEAR, 1, data);
1289    int index = plane_render->add_plane(plane[0], tf);
1290    plane_render->set_active_plane(index);
1291
1292    offscreen_buffer_capture();
1293    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear screen
1294    plane_render->render();
1295
1296    // INSOO
1297    glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, screen_buffer);
1298    //glReadPixels(0, 0, width, height, GL_BGR, GL_UNSIGNED_BYTE, screen_buffer); // INSOO's
1299
1300    std::ostringstream result;
1301    result << "nv>legend " << volArg;
1302    result << " " << volume[ivol]->range_min();
1303    result << " " << volume[ivol]->range_max();
1304    bmp_write(result.str().c_str());
1305    write(0, "\n", 1);
1306
1307    plane_render->remove_plane(index);
1308    resize_offscreen_buffer(old_width, old_height);
1309
1310    return TCL_OK;
1311}
1312
1313//initialize frame buffer objects for offscreen rendering
1314void init_offscreen_buffer()
1315{
1316
1317  //initialize final fbo for final display
1318  glGenFramebuffersEXT(1, &final_fbo);
1319  glGenTextures(1, &final_color_tex);
1320  glGenRenderbuffersEXT(1, &final_depth_rb);
1321
1322  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
1323
1324  //initialize final color texture
1325  glBindTexture(GL_TEXTURE_2D, final_color_tex);
1326  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1327  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
1328#ifdef NV40
1329  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, win_width, win_height, 0,
1330               GL_RGB, GL_INT, NULL);
1331#else
1332  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, win_width, win_height, 0,
1333               GL_RGB, GL_INT, NULL);
1334#endif
1335  glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,
1336                            GL_COLOR_ATTACHMENT0_EXT,
1337                            GL_TEXTURE_2D, final_color_tex, 0);
1338
1339  // initialize final depth renderbuffer
1340  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
1341  glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT,
1342                           GL_DEPTH_COMPONENT24, win_width, win_height);
1343  glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT,
1344                               GL_DEPTH_ATTACHMENT_EXT,
1345                               GL_RENDERBUFFER_EXT, final_depth_rb);
1346
1347  // Check framebuffer completeness at the end of initialization.
1348  CHECK_FRAMEBUFFER_STATUS();
1349  assert(glGetError()==0);
1350}
1351
1352
1353//resize the offscreen buffer
1354void resize_offscreen_buffer(int w, int h){
1355  win_width = w;
1356  win_height = h;
1357
1358    if (g_fonts)
1359    {
1360        g_fonts->resize(w, h);
1361    }
1362
1363  //fprintf(stderr, "screen_buffer size: %d\n", sizeof(screen_buffer));
1364  printf("screen_buffer size: %d %d\n", w, h);
1365
1366  if (screen_buffer) {
1367      delete[] screen_buffer;
1368      screen_buffer = NULL;
1369  }
1370
1371  screen_buffer = new unsigned char[4*win_width*win_height];
1372  assert(screen_buffer != NULL);
1373
1374  //delete the current render buffer resources
1375  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
1376  glDeleteTextures(1, &final_color_tex);
1377  glDeleteFramebuffersEXT(1, &final_fbo);
1378
1379  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
1380  glDeleteRenderbuffersEXT(1, &final_depth_rb);
1381fprintf(stdin,"  after glDeleteRenderbuffers\n");
1382
1383  //change the camera setting
1384  cam->set_screen_size(0, 0, win_width, win_height);
1385  plane_render->set_screen_size(win_width, win_height);
1386
1387  //Reinitialize final fbo for final display
1388  glGenFramebuffersEXT(1, &final_fbo);
1389  glGenTextures(1, &final_color_tex);
1390  glGenRenderbuffersEXT(1, &final_depth_rb);
1391fprintf(stdin,"  after glGenRenderbuffers\n");
1392
1393  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
1394
1395  //initialize final color texture
1396  glBindTexture(GL_TEXTURE_2D, final_color_tex);
1397  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1398  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
1399#ifdef NV40
1400  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, win_width, win_height, 0,
1401               GL_RGB, GL_INT, NULL);
1402#else
1403  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, win_width, win_height, 0,
1404               GL_RGB, GL_INT, NULL);
1405#endif
1406  glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,
1407                            GL_COLOR_ATTACHMENT0_EXT,
1408                            GL_TEXTURE_2D, final_color_tex, 0);
1409fprintf(stdin,"  after glFramebufferTexture2DEXT\n");
1410       
1411  // initialize final depth renderbuffer
1412  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, final_depth_rb);
1413  glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT,
1414                           GL_DEPTH_COMPONENT24, win_width, win_height);
1415  glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT,
1416                               GL_DEPTH_ATTACHMENT_EXT,
1417                               GL_RENDERBUFFER_EXT, final_depth_rb);
1418fprintf(stdin,"  after glFramebufferRenderEXT\n");
1419
1420  // Check framebuffer completeness at the end of initialization.
1421  CHECK_FRAMEBUFFER_STATUS();
1422  assert(glGetError()==0);
1423fprintf(stdin,"  after assert\n");
1424}
1425
1426
1427//init line integral convolution
1428void init_lic(){
1429  lic = new Lic(NMESH, win_width, win_height, 0.3, g_context, volume[1]->id,
1430                  volume[1]->aspect_ratio_width,
1431                  volume[1]->aspect_ratio_height,
1432                  volume[1]->aspect_ratio_depth);
1433}
1434
1435//init the particle system using vector field volume->[1]
1436/*
1437void init_particle_system()
1438{
1439    psys = new ParticleSystem(NMESH, NMESH, g_context, volume[1]->id,
1440        1./volume[1]->aspect_ratio_width,
1441        1./volume[1]->aspect_ratio_height,
1442        1./volume[1]->aspect_ratio_depth);
1443
1444    init_particles();   //fill initial particles
1445}
1446*/
1447
1448
1449void make_test_2D_data()
1450{
1451
1452  int w = 300;
1453  int h = 200;
1454  float* data = new float[w*h];
1455
1456  //procedurally make a gradient plane
1457  for(int j=0; j<h; j++){
1458    for(int i=0; i<w; i++){
1459      data[w*j+i] = float(i)/float(w);
1460    }
1461  }
1462
1463  plane[0] = new Texture2D(w, h, GL_FLOAT, GL_LINEAR, 1, data);
1464
1465  delete[] data;
1466}
1467
1468
1469/*----------------------------------------------------*/
1470void initGL(void)
1471{
1472   //buffer to store data read from the screen
1473   if (screen_buffer) {
1474       delete[] screen_buffer;
1475       screen_buffer = NULL;
1476   }
1477   screen_buffer = new unsigned char[4*win_width*win_height];
1478   assert(screen_buffer != NULL);
1479
1480   //create the camera with default setting
1481   cam = new Camera(0, 0, win_width, win_height,
1482                   live_obj_x, live_obj_y, live_obj_z,
1483                   0., 0., 100.,
1484                   (int)live_rot_x, (int)live_rot_y, (int)live_rot_z);
1485
1486   glEnable(GL_TEXTURE_2D);
1487   glShadeModel(GL_FLAT);
1488   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
1489
1490   glClearColor(0,0,0,1);
1491   glClear(GL_COLOR_BUFFER_BIT);
1492
1493   //initialize lighting
1494   GLfloat mat_specular[] = {1.0, 1.0, 1.0, 1.0};
1495   GLfloat mat_shininess[] = {30.0};
1496   GLfloat white_light[] = {1.0, 1.0, 1.0, 1.0};
1497   GLfloat green_light[] = {0.1, 0.5, 0.1, 1.0};
1498
1499   glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
1500   glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
1501   glLightfv(GL_LIGHT0, GL_DIFFUSE, white_light);
1502   glLightfv(GL_LIGHT0, GL_SPECULAR, white_light);     
1503   glLightfv(GL_LIGHT1, GL_DIFFUSE, green_light);
1504   glLightfv(GL_LIGHT1, GL_SPECULAR, white_light);     
1505
1506   // init table of transfer functions
1507   Tcl_InitHashTable(&tftable, TCL_STRING_KEYS);
1508
1509   //check if performance query is supported
1510   if(check_query_support()){
1511     //create queries to count number of rendered pixels
1512     perf = new PerfQuery();
1513   }
1514
1515   init_offscreen_buffer();    //frame buffer object for offscreen rendering
1516
1517   //create volume renderer and add volumes to it
1518   g_vol_render = new VolumeRenderer();
1519
1520   /*
1521   //I added this to debug : Wei
1522   float tmp_data[4*124];
1523   memset(tmp_data, 0, 4*4*124);
1524   TransferFunction* tmp_tf = new TransferFunction(124, tmp_data);
1525   g_vol_render->add_volume(volume[0], tmp_tf);
1526   volume[0]->get_cutplane(0)->enabled = false;
1527   volume[0]->get_cutplane(1)->enabled = false;
1528   volume[0]->get_cutplane(2)->enabled = false;
1529
1530   //volume[1]->move(Vector3(0.5, 0.6, 0.7));
1531   //g_vol_render->add_volume(volume[1], tmp_tf);
1532   */
1533
1534   //create an 2D plane renderer
1535   plane_render = new PlaneRenderer(g_context, win_width, win_height);
1536   make_test_2D_data();
1537
1538   plane_render->add_plane(plane[0], get_transfunc("default"));
1539
1540   assert(glGetError()==0);
1541
1542   //init_particle_system();
1543   //init_lic();
1544}
1545
1546void read_screen()
1547{
1548  glReadPixels(0, 0, win_width, win_height, GL_RGB, GL_UNSIGNED_BYTE, screen_buffer);
1549  assert(glGetError()==0);
1550}
1551
1552void display();
1553
1554
1555#if DO_RLE
1556char rle[512*512*3];
1557int rle_size;
1558
1559short offsets[512*512*3];
1560int offsets_size;
1561
1562void do_rle(){
1563  int len = win_width*win_height*3;
1564  rle_size = 0;
1565  offsets_size = 0;
1566
1567  int i=0;
1568  while(i<len){
1569    if (screen_buffer[i] == 0) {
1570      int pos = i+1;
1571      while ( (pos<len) && (screen_buffer[pos] == 0)) {
1572        pos++;
1573      }
1574      offsets[offsets_size++] = -(pos - i);
1575      i = pos;
1576    }
1577
1578    else {
1579      int pos;
1580      for (pos = i; (pos<len) && (screen_buffer[pos] != 0); pos++) {
1581        rle[rle_size++] = screen_buffer[pos];
1582      }
1583      offsets[offsets_size++] = (pos - i);
1584      i = pos;
1585    }
1586
1587  }
1588}
1589#endif
1590
1591// used internally to build up the BMP file header
1592// writes an integer value into the header data structure at pos
1593void
1594bmp_header_add_int(unsigned char* header, int& pos, int data)
1595{
1596    header[pos++] = data & 0xff;
1597    header[pos++] = (data >> 8) & 0xff;
1598    header[pos++] = (data >> 16) & 0xff;
1599    header[pos++] = (data >> 24) & 0xff;
1600}
1601
1602// INSOO
1603// FOR DEBUGGING
1604void
1605bmp_write_to_file()
1606{
1607    unsigned char header[54];
1608    int pos = 0;
1609    header[pos++] = 'B';
1610    header[pos++] = 'M';
1611
1612    // BE CAREFUL:  BMP files must have an even multiple of 4 bytes
1613    // on each scan line.  If need be, we add padding to each line.
1614    int pad = 0;
1615    if ((3*win_width) % 4 > 0) {
1616        pad = 4 - ((3*win_width) % 4);
1617    }
1618
1619    // file size in bytes
1620    int fsize = (3*win_width+pad)*win_height + sizeof(header);
1621    bmp_header_add_int(header, pos, fsize);
1622
1623    // reserved value (must be 0)
1624    bmp_header_add_int(header, pos, 0);
1625
1626    // offset in bytes to start of bitmap data
1627    bmp_header_add_int(header, pos, sizeof(header));
1628
1629    // size of the BITMAPINFOHEADER
1630    bmp_header_add_int(header, pos, 40);
1631
1632    // width of the image in pixels
1633    bmp_header_add_int(header, pos, win_width);
1634
1635    // height of the image in pixels
1636    bmp_header_add_int(header, pos, win_height);
1637
1638    // 1 plane + (24 bits/pixel << 16)
1639    bmp_header_add_int(header, pos, 1572865);
1640
1641    // no compression
1642    // size of image for compression
1643    bmp_header_add_int(header, pos, 0);
1644    bmp_header_add_int(header, pos, 0);
1645
1646    // x pixels per meter
1647    // y pixels per meter
1648    bmp_header_add_int(header, pos, 0);
1649    bmp_header_add_int(header, pos, 0);
1650
1651    // number of colors used (0 = compute from bits/pixel)
1652    // number of important colors (0 = all colors important)
1653    bmp_header_add_int(header, pos, 0);
1654    bmp_header_add_int(header, pos, 0);
1655
1656    // BE CAREFUL: BMP format wants BGR ordering for screen data
1657    unsigned char* scr = screen_buffer;
1658    for (int row=0; row < win_height; row++) {
1659        for (int col=0; col < win_width; col++) {
1660            unsigned char tmp = scr[2];
1661            scr[2] = scr[0];  // B
1662            scr[0] = tmp;     // R
1663            scr += 3;
1664        }
1665        scr += pad;  // skip over padding already in screen data
1666    }
1667
1668    FILE* f;
1669    f = fopen("/tmp/image.bmp", "wb");
1670    fwrite((void*) header, sizeof(header), 1, f);
1671    fwrite((void*) screen_buffer, (3*win_width+pad)*win_height, 1, f);
1672    fclose(f);
1673}
1674
1675void
1676bmp_write(const char* cmd)
1677{
1678    unsigned char header[54];
1679    int pos = 0;
1680    header[pos++] = 'B';
1681    header[pos++] = 'M';
1682
1683    // BE CAREFUL:  BMP files must have an even multiple of 4 bytes
1684    // on each scan line.  If need be, we add padding to each line.
1685    int pad = 0;
1686    if ((3*win_width) % 4 > 0) {
1687        pad = 4 - ((3*win_width) % 4);
1688    }
1689
1690    // file size in bytes
1691    int fsize = (3*win_width+pad)*win_height + sizeof(header);
1692    bmp_header_add_int(header, pos, fsize);
1693
1694    // reserved value (must be 0)
1695    bmp_header_add_int(header, pos, 0);
1696
1697    // offset in bytes to start of bitmap data
1698    bmp_header_add_int(header, pos, sizeof(header));
1699
1700    // size of the BITMAPINFOHEADER
1701    bmp_header_add_int(header, pos, 40);
1702
1703    // width of the image in pixels
1704    bmp_header_add_int(header, pos, win_width);
1705
1706    // height of the image in pixels
1707    bmp_header_add_int(header, pos, win_height);
1708
1709    // 1 plane + (24 bits/pixel << 16)
1710    bmp_header_add_int(header, pos, 1572865);
1711
1712    // no compression
1713    // size of image for compression
1714    bmp_header_add_int(header, pos, 0);
1715    bmp_header_add_int(header, pos, 0);
1716
1717    // x pixels per meter
1718    // y pixels per meter
1719    bmp_header_add_int(header, pos, 0);
1720    bmp_header_add_int(header, pos, 0);
1721
1722    // number of colors used (0 = compute from bits/pixel)
1723    // number of important colors (0 = all colors important)
1724    bmp_header_add_int(header, pos, 0);
1725    bmp_header_add_int(header, pos, 0);
1726
1727    // BE CAREFUL: BMP format wants BGR ordering for screen data
1728    unsigned char* scr = screen_buffer;
1729    for (int row=0; row < win_height; row++) {
1730        for (int col=0; col < win_width; col++) {
1731            unsigned char tmp = scr[2];
1732            scr[2] = scr[0];  // B
1733            scr[0] = tmp;     // R
1734            scr += 3;
1735        }
1736        scr += pad;  // skip over padding already in screen data
1737    }
1738
1739    std::ostringstream result;
1740    result << cmd << " " << fsize << "\n";
1741    write(0, result.str().c_str(), result.str().size());
1742
1743    write(0, header, sizeof(header));
1744    write(0, screen_buffer, (3*win_width+pad)*win_height);
1745
1746}
1747
1748/*
1749//draw vectors
1750void draw_arrows(){
1751  glColor4f(0.,0.,1.,1.);
1752  for(int i=0; i<NMESH; i++){
1753    for(int j=0; j<NMESH; j++){
1754      Vector2 v = grid.get(i, j);
1755
1756      int x1 = i*DM;
1757      int y1 = j*DM;
1758
1759      int x2 = x1 + v.x;
1760      int y2 = y1 + v.y;
1761
1762      glBegin(GL_LINES);
1763        glVertex2d(x1, y1);
1764        glVertex2d(x2, y2);
1765      glEnd();
1766    }
1767  }
1768}
1769*/
1770
1771
1772/*----------------------------------------------------*/
1773void idle()
1774{
1775    glutSetWindow(render_window);
1776 
1777  /*
1778  struct timespec ts;
1779  ts.tv_sec = 0;
1780  ts.tv_nsec = 300000000;
1781  nanosleep(&ts, 0);
1782  */
1783
1784#ifdef XINETD
1785    // in Command.cpp
1786    xinetd_listen();
1787#else
1788    glutPostRedisplay();
1789#endif
1790}
1791
1792
1793void display_offscreen_buffer()
1794{
1795   glEnable(GL_TEXTURE_2D);
1796   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
1797   glBindTexture(GL_TEXTURE_2D, final_color_tex);
1798
1799   glViewport(0, 0, win_width, win_height);
1800   glMatrixMode(GL_PROJECTION);
1801   glLoadIdentity();
1802   gluOrtho2D(0, win_width, 0, win_height);
1803   glMatrixMode(GL_MODELVIEW);
1804   glLoadIdentity();
1805
1806   glColor3f(1.,1.,1.);         //MUST HAVE THIS LINE!!!
1807   glBegin(GL_QUADS);
1808   glTexCoord2f(0, 0); glVertex2f(0, 0);
1809   glTexCoord2f(1, 0); glVertex2f(win_width, 0);
1810   glTexCoord2f(1, 1); glVertex2f(win_width, win_height);
1811   glTexCoord2f(0, 1); glVertex2f(0, win_height);
1812   glEnd();
1813
1814}
1815
1816void offscreen_buffer_capture()
1817{
1818   glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, final_fbo);
1819}
1820
1821void draw_bounding_box(float x0, float y0, float z0,
1822                float x1, float y1, float z1,
1823                float r, float g, float b, float line_width)
1824{
1825        glDisable(GL_TEXTURE_2D);
1826
1827        glColor4d(r, g, b, 1.0);
1828        glLineWidth(line_width);
1829       
1830        glBegin(GL_LINE_LOOP);
1831
1832                glVertex3d(x0, y0, z0);
1833                glVertex3d(x1, y0, z0);
1834                glVertex3d(x1, y1, z0);
1835                glVertex3d(x0, y1, z0);
1836               
1837        glEnd();
1838
1839        glBegin(GL_LINE_LOOP);
1840
1841                glVertex3d(x0, y0, z1);
1842                glVertex3d(x1, y0, z1);
1843                glVertex3d(x1, y1, z1);
1844                glVertex3d(x0, y1, z1);
1845               
1846        glEnd();
1847
1848
1849        glBegin(GL_LINE_LOOP);
1850
1851                glVertex3d(x0, y0, z0);
1852                glVertex3d(x0, y0, z1);
1853                glVertex3d(x0, y1, z1);
1854                glVertex3d(x0, y1, z0);
1855               
1856        glEnd();
1857
1858        glBegin(GL_LINE_LOOP);
1859
1860                glVertex3d(x1, y0, z0);
1861                glVertex3d(x1, y0, z1);
1862                glVertex3d(x1, y1, z1);
1863                glVertex3d(x1, y1, z0);
1864               
1865        glEnd();
1866
1867        glEnable(GL_TEXTURE_2D);
1868}
1869
1870
1871
1872int particle_distance_sort(const void* a, const void* b){
1873  if((*((Particle*)a)).aux > (*((Particle*)b)).aux)
1874    return -1;
1875  else
1876    return 1;
1877}
1878
1879/*
1880void soft_read_verts()
1881{
1882  glReadPixels(0, 0, psys->psys_width, psys->psys_height, GL_RGB, GL_FLOAT, vert);
1883  //fprintf(stderr, "soft_read_vert");
1884
1885  //cpu sort the distance 
1886  Particle* p = (Particle*) malloc(sizeof(Particle)* psys->psys_width * psys->psys_height);
1887  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
1888    float x = vert[3*i];
1889    float y = vert[3*i+1];
1890    float z = vert[3*i+2];
1891
1892    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);
1893    p[i].x = x;
1894    p[i].y = y;
1895    p[i].z = z;
1896    p[i].aux = dis;
1897  }
1898
1899  qsort(p, psys->psys_width * psys->psys_height, sizeof(Particle), particle_distance_sort);
1900
1901  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
1902    vert[3*i] = p[i].x;
1903    vert[3*i+1] = p[i].y;
1904    vert[3*i+2] = p[i].z;
1905  }
1906
1907  free(p);
1908}
1909*/
1910
1911/*
1912//display the content of a texture as a screen aligned quad
1913void display_texture(NVISid tex, int width, int height)
1914{
1915   glPushMatrix();
1916
1917   glEnable(GL_TEXTURE_2D);
1918   glBindTexture(GL_TEXTURE_RECTANGLE_NV, tex);
1919
1920   glViewport(0, 0, width, height);
1921   glMatrixMode(GL_PROJECTION);
1922   glLoadIdentity();
1923   gluOrtho2D(0, width, 0, height);
1924   glMatrixMode(GL_MODELVIEW);
1925   glLoadIdentity();
1926
1927   cgGLBindProgram(m_passthru_fprog);
1928   cgGLEnableProfile(CG_PROFILE_FP30);
1929
1930   cgGLSetParameter4f(m_passthru_scale_param, 1.0, 1.0, 1.0, 1.0);
1931   cgGLSetParameter4f(m_passthru_bias_param, 0.0, 0.0, 0.0, 0.0);
1932   
1933   draw_quad(width, height, width, height);
1934   cgGLDisableProfile(CG_PROFILE_FP30);
1935
1936   glPopMatrix();
1937
1938   assert(glGetError()==0);
1939}
1940*/
1941
1942
1943//draw vertices in the main memory
1944/*
1945void soft_display_verts()
1946{
1947  glDisable(GL_TEXTURE_2D);
1948  glEnable(GL_BLEND);
1949
1950  glPointSize(0.5);
1951  glColor4f(0,0.8,0.8,0.6);
1952  glBegin(GL_POINTS);
1953  for(int i=0; i<psys->psys_width * psys->psys_height; i++){
1954    glVertex3f(vert[3*i], vert[3*i+1], vert[3*i+2]);
1955  }
1956  glEnd();
1957  //fprintf(stderr, "soft_display_vert");
1958}
1959*/
1960
1961#if 0
1962
1963//oddeven sort on GPU
1964void sortstep()
1965{
1966    // perform one step of the current sorting algorithm
1967
1968        /*
1969    // swap buffers
1970    int sourceBuffer = targetBuffer;
1971    targetBuffer = (targetBuffer+1)%2;   
1972    int pstage = (1<<stage);
1973    int ppass  = (1<<pass);
1974    int activeBitonicShader = 0;
1975
1976#ifdef _WIN32
1977    buffer->BindBuffer(wglTargets[sourceBuffer]);
1978#else
1979    buffer->BindBuffer(glTargets[sourceBuffer]);
1980#endif
1981    if (buffer->IsDoubleBuffered()) glDrawBuffer(glTargets[targetBuffer]);
1982    */
1983
1984    checkGLError("after db");
1985
1986    int pstage = (1<<stage);
1987    int ppass  = (1<<pass);
1988    //int activeBitonicShader = 0;
1989
1990    // switch on correct sorting shader
1991    oddevenMergeSort.bind();
1992    glUniform3fARB(oddevenMergeSort.getUniformLocation("Param1"), float(pstage+pstage),
1993                   float(ppass%pstage), float((pstage+pstage)-(ppass%pstage)-1));
1994    glUniform3fARB(oddevenMergeSort.getUniformLocation("Param2"), float(psys_width), float(psys_height), float(ppass));
1995    glUniform1iARB(oddevenMergeSort.getUniformLocation("Data"), 0);
1996    staticdebugmsg("sort","stage "<<pstage<<" pass "<<ppass);
1997
1998    // This clear is not necessary for sort to function. But if we are in interactive mode
1999    // unused parts of the texture that are visible will look bad.
2000    //if (!perfTest) glClear(GL_COLOR_BUFFER_BIT);
2001
2002    //buffer->Bind();
2003    //buffer->EnableTextureTarget();
2004
2005    // initiate the sorting step on the GPU
2006    // a full-screen quad
2007    glBegin(GL_QUADS);
2008    /*
2009    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,0.0f,0.0f,1.0f); glVertex2f(-1.0f,-1.0f);
2010    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),0.0f,1.0f,1.0f); glVertex2f(1.0f,-1.0f);
2011    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),float(psys_height),1.0f,0.0f); glVertex2f(1.0f,1.0f);       
2012    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,float(psys_height),0.0f,0.0f); glVertex2f(-1.0f,1.0f);   
2013    */
2014    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,0.0f,0.0f,1.0f); glVertex2f(0.,0.);       
2015    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),0.0f,1.0f,1.0f); glVertex2f(float(psys_width), 0.);
2016    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,float(psys_width),float(psys_height),1.0f,0.0f); glVertex2f(float(psys_width), float(psys_height));   
2017    glMultiTexCoord4fARB(GL_TEXTURE0_ARB,0.0f,float(psys_height),0.0f,0.0f); glVertex2f(0., float(psys_height));       
2018    glEnd();
2019
2020
2021    // switch off sorting shader
2022    oddevenMergeSort.release();
2023
2024    //buffer->DisableTextureTarget();
2025
2026    assert(glGetError()==0);
2027}
2028
2029#endif
2030
2031
2032void draw_3d_axis()
2033{
2034    glDisable(GL_TEXTURE_2D);
2035    glEnable(GL_DEPTH_TEST);
2036
2037        //draw axes
2038        GLUquadric *obj;
2039        obj = gluNewQuadric();
2040       
2041        glDepthFunc(GL_LESS);
2042    glEnable(GL_COLOR_MATERIAL);
2043        glEnable(GL_DEPTH_TEST);
2044        glDisable(GL_BLEND);
2045
2046        int segments = 50;
2047
2048        glColor3f(0.8, 0.8, 0.8);
2049        glPushMatrix();
2050        glTranslatef(0.4, 0., 0.);
2051        glRotatef(90, 1, 0, 0);
2052        glRotatef(180, 0, 1, 0);
2053        glScalef(0.0005, 0.0005, 0.0005);
2054        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'x');
2055        glPopMatrix(); 
2056
2057        glPushMatrix();
2058        glTranslatef(0., 0.4, 0.);
2059        glRotatef(90, 1, 0, 0);
2060        glRotatef(180, 0, 1, 0);
2061        glScalef(0.0005, 0.0005, 0.0005);
2062        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'y');
2063        glPopMatrix(); 
2064
2065        glPushMatrix();
2066        glTranslatef(0., 0., 0.4);
2067        glRotatef(90, 1, 0, 0);
2068        glRotatef(180, 0, 1, 0);
2069        glScalef(0.0005, 0.0005, 0.0005);
2070        glutStrokeCharacter(GLUT_STROKE_ROMAN, 'z');
2071        glPopMatrix(); 
2072
2073        glEnable(GL_LIGHTING);
2074        glEnable(GL_LIGHT0);
2075
2076        //glColor3f(0.2, 0.2, 0.8);
2077        glPushMatrix();
2078        glutSolidSphere(0.02, segments, segments );
2079        glPopMatrix();
2080
2081        glPushMatrix();
2082        glRotatef(-90, 1, 0, 0);       
2083        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
2084        glPopMatrix(); 
2085
2086        glPushMatrix();
2087        glTranslatef(0., 0.3, 0.);
2088        glRotatef(-90, 1, 0, 0);       
2089        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
2090        glPopMatrix(); 
2091
2092        glPushMatrix();
2093        glRotatef(90, 0, 1, 0);
2094        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
2095        glPopMatrix(); 
2096
2097        glPushMatrix();
2098        glTranslatef(0.3, 0., 0.);
2099        glRotatef(90, 0, 1, 0);
2100        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
2101        glPopMatrix(); 
2102
2103        glPushMatrix();
2104        gluCylinder(obj, 0.01, 0.01, 0.3, segments, segments);
2105        glPopMatrix(); 
2106
2107        glPushMatrix();
2108        glTranslatef(0., 0., 0.3);
2109        gluCylinder(obj, 0.02, 0.0, 0.06, segments, segments);
2110        glPopMatrix(); 
2111
2112        glDisable(GL_LIGHTING);
2113        glDisable(GL_DEPTH_TEST);
2114        gluDeleteQuadric(obj);
2115
2116    glEnable(GL_TEXTURE_2D);
2117    glDisable(GL_DEPTH_TEST);
2118}
2119
2120/*
2121void draw_axis()
2122{
2123  glDisable(GL_TEXTURE_2D);
2124  glEnable(GL_DEPTH_TEST);
2125
2126  //red x
2127  glColor3f(1,0,0);
2128  glBegin(GL_LINES);
2129    glVertex3f(0,0,0);
2130    glVertex3f(1.5,0,0);
2131  glEnd();
2132
2133  //blue y
2134  glColor3f(0,0,1);
2135  glBegin(GL_LINES);
2136    glVertex3f(0,0,0);
2137    glVertex3f(0,1.5,0);
2138  glEnd();
2139
2140  //green z
2141  glColor3f(0,1,0);
2142  glBegin(GL_LINES);
2143    glVertex3f(0,0,0);
2144    glVertex3f(0,0,1.5);
2145  glEnd();
2146
2147  glEnable(GL_TEXTURE_2D);
2148  glDisable(GL_DEPTH_TEST);
2149}
2150*/
2151
2152
2153
2154/*----------------------------------------------------*/
2155void display()
2156{
2157    assert(glGetError()==0);
2158
2159    //lic->convolve(); //flow line integral convolution
2160    //psys->advect(); //advect particles
2161
2162    //start final rendering
2163    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear screen
2164
2165    if (volume_mode)
2166    {
2167        //3D rendering mode
2168        glEnable(GL_TEXTURE_2D);
2169        glEnable(GL_DEPTH_TEST);
2170
2171        //camera setting activated
2172        cam->activate();
2173
2174        //set up the orientation of items in the scene.
2175        glPushMatrix();
2176        switch (updir) {
2177        case 1:  // x
2178            glRotatef(90, 0, 0, 1);
2179            glRotatef(90, 1, 0, 0);
2180            break;
2181
2182        case 2:  // y
2183            // this is the default
2184            break;
2185
2186        case 3:  // z
2187            glRotatef(-90, 1, 0, 0);
2188            glRotatef(-90, 0, 0, 1);
2189            break;
2190
2191        case -1:  // -x
2192            glRotatef(-90, 0, 0, 1);
2193            break;
2194
2195        case -2:  // -y
2196            glRotatef(180, 0, 0, 1);
2197            glRotatef(-90, 0, 1, 0);
2198            break;
2199
2200        case -3:  // -z
2201            glRotatef(90, 1, 0, 0);
2202            break;
2203        }
2204
2205        glPushMatrix();
2206        //now render things in the scene
2207        if (axis_on)
2208        {
2209            draw_3d_axis();
2210        }
2211
2212        if (g_grid->isVisible())
2213        {
2214            g_grid->render();
2215        }
2216
2217        //lic->render();        //display the line integral convolution result
2218        //soft_display_verts();
2219        //perf->enable();
2220        //psys->render();
2221        //perf->disable();
2222        //fprintf(stderr, "particle pixels: %d\n", perf->get_pixel_count());
2223        //perf->reset();
2224
2225        perf->enable();
2226        g_vol_render->render_all();
2227
2228        perf->disable();
2229
2230        for (int ii = 0; ii < g_heightMap.size(); ++ii)
2231        {
2232            if (g_heightMap[ii]->isVisible())
2233                g_heightMap[ii]->render();
2234        }
2235
2236        glPopMatrix();
2237
2238        float mat[16];
2239        glGetFloatv(GL_MODELVIEW_MATRIX, mat);
2240
2241        for (int i = 0; i < g_pointSet.size(); ++i)
2242        {
2243            if (g_pointSet[i]->isVisible())
2244            {
2245                g_pointset_renderer->render(g_pointSet[i]->getCluster(), mat, g_pointSet[i]->getSortLevel(),
2246                                            g_pointSet[i]->getScale(),
2247                                            g_pointSet[i]->getOrigin());
2248            }
2249        }
2250
2251        glPopMatrix();
2252   }
2253   else {
2254        //2D rendering mode
2255        perf->enable();
2256        plane_render->render();
2257        perf->disable();
2258   }
2259
2260#ifdef XINETD
2261    float cost  = perf->get_pixel_count();
2262    write(3, &cost, sizeof(cost));
2263#endif
2264    perf->reset();
2265
2266}
2267
2268
2269void mouse(int button, int state, int x, int y){
2270  if(button==GLUT_LEFT_BUTTON){
2271
2272    if(state==GLUT_DOWN){
2273      left_last_x = x;
2274      left_last_y = y;
2275      left_down = true;
2276      right_down = false;
2277    }
2278    else{
2279      left_down = false;
2280      right_down = false;
2281    }
2282  }
2283  else{
2284    //fprintf(stderr, "right mouse\n");
2285
2286    if(state==GLUT_DOWN){
2287      //fprintf(stderr, "right mouse down\n");
2288      right_last_x = x;
2289      right_last_y = y;
2290      left_down = false;
2291      right_down = true;
2292    }
2293    else{
2294      //fprintf(stderr, "right mouse up\n");
2295      left_down = false;
2296      right_down = false;
2297    }
2298  }
2299}
2300
2301
2302void update_rot(int delta_x, int delta_y){
2303        live_rot_x += delta_x;
2304        live_rot_y += delta_y;
2305
2306        if(live_rot_x > 360.0)
2307                live_rot_x -= 360.0;   
2308        else if(live_rot_x < -360.0)
2309                live_rot_x += 360.0;
2310
2311        if(live_rot_y > 360.0)
2312                live_rot_y -= 360.0;   
2313        else if(live_rot_y < -360.0)
2314                live_rot_y += 360.0;
2315
2316        cam->rotate(live_rot_x, live_rot_y, live_rot_z);
2317}
2318
2319
2320void update_trans(int delta_x, int delta_y, int delta_z){
2321        live_obj_x += delta_x*0.03;
2322        live_obj_y += delta_y*0.03;
2323        live_obj_z += delta_z*0.03;
2324}
2325
2326void end_service();
2327
2328void keyboard(unsigned char key, int x, int y){
2329 
2330   bool log = false;
2331
2332   switch (key){
2333        case 'q':
2334#ifdef XINETD
2335        //end_service();
2336        NvExitService();
2337#endif
2338                exit(0);
2339                break;
2340        case '+':
2341                lic_slice_z+=0.05;
2342                lic->set_offset(lic_slice_z);
2343                break;
2344        case '-':
2345                lic_slice_z-=0.05;
2346                lic->set_offset(lic_slice_z);
2347                break;
2348        case ',':
2349                lic_slice_x+=0.05;
2350                //init_particles();
2351                break;
2352        case '.':
2353                lic_slice_x-=0.05;
2354                //init_particles();
2355                break;
2356        case '1':
2357                //advect = true;
2358                break;
2359        case '2':
2360                //psys_x+=0.05;
2361                break;
2362        case '3':
2363                //psys_x-=0.05;
2364                break;
2365        case 'w': //zoom out
2366                live_obj_z-=0.05;
2367                log = true;
2368                cam->move(live_obj_x, live_obj_y, live_obj_z);
2369                break;
2370        case 's': //zoom in
2371                live_obj_z+=0.05;
2372                log = true;
2373                cam->move(live_obj_x, live_obj_y, live_obj_z);
2374                break;
2375        case 'a': //left
2376                live_obj_x-=0.05;
2377                log = true;
2378                cam->move(live_obj_x, live_obj_y, live_obj_z);
2379                break;
2380        case 'd': //right
2381                live_obj_x+=0.05;
2382                log = true;
2383                cam->move(live_obj_x, live_obj_y, live_obj_z);
2384                break;
2385        case 'i':
2386                //init_particles();
2387                break;
2388        case 'v':
2389                g_vol_render->switch_volume_mode();
2390                break;
2391        case 'b':
2392                g_vol_render->switch_slice_mode();
2393                break;
2394        case 'n':
2395                resize_offscreen_buffer(win_width*2, win_height*2);
2396                break;
2397        case 'm':
2398                resize_offscreen_buffer(win_width/2, win_height/2);
2399                break;
2400
2401        default:
2402                break;
2403    }   
2404
2405#ifdef EVENTLOG
2406   if(log){
2407     float param[3] = {live_obj_x, live_obj_y, live_obj_z};
2408     Event* tmp = new Event(EVENT_MOVE, param, NvGetTimeInterval());
2409     tmp->write(event_log);
2410     delete tmp;
2411   }
2412#endif
2413}
2414
2415void motion(int x, int y)
2416{
2417
2418    int old_x, old_y;   
2419
2420    if(left_down){
2421      old_x = left_last_x;
2422      old_y = left_last_y;   
2423    }
2424    else if(right_down){
2425      old_x = right_last_x;
2426      old_y = right_last_y;   
2427    }
2428
2429    int delta_x = x - old_x;
2430    int delta_y = y - old_y;
2431
2432    //more coarse event handling
2433    //if(abs(delta_x)<10 && abs(delta_y)<10)
2434      //return;
2435
2436    if(left_down){
2437      left_last_x = x;
2438      left_last_y = y;
2439
2440      update_rot(-delta_y, -delta_x);
2441    }
2442    else if (right_down){
2443      //fprintf(stderr, "right mouse motion (%d,%d)\n", x, y);
2444
2445      right_last_x = x;
2446      right_last_y = y;
2447
2448      update_trans(0, 0, delta_x);
2449    }
2450
2451#ifdef EVENTLOG
2452    float param[3] = {live_rot_x, live_rot_y, live_rot_z};
2453    Event* tmp = new Event(EVENT_ROTATE, param, NvGetTimeInterval());
2454    tmp->write(event_log);
2455    delete tmp;
2456#endif
2457
2458    glutPostRedisplay();
2459}
2460
2461/*
2462#ifdef XINETD
2463void init_service(){
2464  //open log and map stderr to log file
2465  xinetd_log = fopen("/tmp/log.txt", "w");
2466  close(2);
2467  dup2(fileno(xinetd_log), 2);
2468  dup2(2,1);
2469
2470  //flush junk
2471  fflush(stdout);
2472  fflush(stderr);
2473}
2474
2475void end_service(){
2476  //close log file
2477  fclose(xinetd_log);
2478}
2479#endif
2480
2481void init_event_log(){
2482  event_log = fopen("event.txt", "w");
2483  assert(event_log!=0);
2484
2485  struct timeval time;
2486  gettimeofday(&time, NULL);
2487  cur_time = time.tv_sec*1000. + time.tv_usec/1000.;
2488}
2489
2490void end_event_log(){
2491  fclose(event_log);
2492}
2493
2494double get_time_interval(){
2495  struct timeval time;
2496  gettimeofday(&time, NULL);
2497  double new_time = time.tv_sec*1000. + time.tv_usec/1000.;
2498
2499  double interval = new_time - cur_time;
2500  cur_time = new_time;
2501  return interval;
2502}
2503*/
2504
2505void removeAllData()
2506{
2507    //
2508}
2509
2510
2511/*----------------------------------------------------*/
2512int main(int argc, char** argv)
2513{
2514
2515    char path [1000];
2516    while(1) {
2517        int c;
2518        int this_option_optind = optind ? optind : 1;
2519        int option_index = 0;
2520        struct option long_options[] = {
2521          // name, has_arg, flag, val
2522          { 0,0,0,0 },
2523        };
2524
2525        c = getopt_long(argc, argv, "+p:", long_options, &option_index);
2526        if (c == -1)
2527            break;
2528
2529        switch(c) {
2530            case 'p':
2531                strncpy(path, optarg, sizeof(path));
2532                break;
2533            default:
2534                fprintf(stderr,"Don't know what option '%c'.\n", c);
2535                return 1;
2536        }
2537    }
2538
2539    R2FilePath::getInstance()->setWorkingDirectory(argc, (const char**) argv);
2540
2541#ifdef XINETD
2542   signal(SIGPIPE,SIG_IGN);
2543   NvInitService();
2544#endif
2545
2546   glutInit(&argc, argv);
2547   glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
2548
2549   glutInitWindowSize(win_width, win_height);
2550
2551   glutInitWindowPosition(10, 10);
2552   render_window = glutCreateWindow(argv[0]);
2553   glutDisplayFunc(display);
2554
2555#ifndef XINETD
2556   glutMouseFunc(mouse);
2557   glutMotionFunc(motion);
2558   glutKeyboardFunc(keyboard);
2559#endif
2560
2561   glutIdleFunc(idle);
2562   glutReshapeFunc(resize_offscreen_buffer);
2563
2564   NvInit(path);
2565   initGL();
2566   initTcl();
2567
2568#ifdef EVENTLOG
2569   NvInitEventLog();
2570#endif
2571    //event loop
2572    glutMainLoop();
2573
2574    removeAllData();
2575
2576    NvExit();
2577
2578    return 0;
2579}
Note: See TracBrowser for help on using the repository browser.