source: trunk/packages/vizservers/nanovis/VolumeRenderer.cpp @ 3464

Last change on this file since 3464 was 3452, checked in by ldelgass, 7 years ago

Remove XINETD define from nanovis. We only support server mode now, no glut
interaction loop with mouse/keyboard handlers. Fixes for trace logging to make
output closer to vtkvis: inlcude function name for trace messages, remove
newlines from format strings in macros since newlines get added by syslog.

  • Property svn:eol-style set to native
File size: 17.9 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * ----------------------------------------------------------------------
4 * VolumeRenderer.cpp : VolumeRenderer class for volume visualization
5 *
6 * ======================================================================
7 *  AUTHOR:  Wei Qiao <qiaow@purdue.edu>
8 *           Purdue Rendering and Perceptualization Lab (PURPL)
9 *
10 *  Copyright (c) 2004-2012  HUBzero Foundation, LLC
11 *
12 *  See the file "license.terms" for information on usage and
13 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
14 * ======================================================================
15 */
16#include <stdlib.h>
17#include <float.h>
18
19#include <vector>
20
21#include <GL/glew.h>
22
23#include <tcl.h>
24
25#include "nanovis.h"
26#include "VolumeRenderer.h"
27#include "ConvexPolygon.h"
28#include "NvStdVertexShader.h"
29#include "Trace.h"
30
31VolumeRenderer::VolumeRenderer()
32{
33    initShaders();
34
35    _volumeInterpolator = new VolumeInterpolator();
36}
37
38VolumeRenderer::~VolumeRenderer()
39{
40    delete _cutplaneShader;
41    delete _zincBlendeShader;
42    delete _regularVolumeShader;
43    delete _stdVertexShader;
44    delete _volumeInterpolator;
45}
46
47//initialize the volume shaders
48void VolumeRenderer::initShaders()
49{
50    _cutplaneShader = new NvShader();
51    _cutplaneShader->loadVertexProgram("cutplane_vp.cg", "main");
52    _cutplaneShader->loadFragmentProgram("cutplane_fp.cg", "main");
53
54    //standard vertex program
55    _stdVertexShader = new NvStdVertexShader();
56
57    //volume rendering shader: one cubic volume
58    _regularVolumeShader = new NvRegularVolumeShader();
59
60    //volume rendering shader: one zincblende orbital volume.
61    //This shader renders one orbital of the simulation.
62    //A sim has S, P, D, SS orbitals. thus a full rendering requires 4 zincblende orbital volumes.
63    //A zincblende orbital volume is decomposed into 2 "interlocking" cubic 4-component volumes and passed to the shader.
64    //We render each orbital with a independent transfer functions then blend the result.
65    //
66    //The engine is already capable of rendering multiple volumes and combine them. Thus, we just invoke this shader on
67    //S, P, D and SS orbitals with different transfor functions. The result is a multi-orbital rendering.
68    _zincBlendeShader = new NvZincBlendeVolumeShader();
69}
70
71struct SortElement {
72    float z;
73    int volumeId;
74    int sliceId;
75
76    SortElement(float _z, int _v, int _s) :
77        z(_z),
78        volumeId(_v),
79        sliceId(_s)
80    {}
81};
82
83static int sliceSort(const void *a, const void *b)
84{
85    if ((*((SortElement*)a)).z > (*((SortElement*)b)).z)
86        return 1;
87    else
88        return -1;
89}
90
91void 
92VolumeRenderer::renderAll()
93{
94    size_t total_rendered_slices = 0;
95
96    if (_volumeInterpolator->isStarted()) {
97#ifdef notdef
98        ani_vol = _volumeInterpolator->getVolume();
99        ani_tf = ani_vol->transferFunction();
100#endif
101        TRACE("VOLUME INTERPOLATOR IS STARTED ----------------------------");
102    }
103    // Determine the volumes that are to be rendered.
104    std::vector<Volume *> volumes;
105    Tcl_HashEntry *hPtr;
106    Tcl_HashSearch iter;
107    for (hPtr = Tcl_FirstHashEntry(&NanoVis::volumeTable, &iter); hPtr != NULL;
108         hPtr = Tcl_NextHashEntry(&iter)) {
109        Volume *volPtr;
110        volPtr = (Volume *)Tcl_GetHashValue(hPtr);
111        if (!volPtr->visible()) {
112            continue; // Skip this volume
113        }
114        // BE CAREFUL: Set the number of slices to something slightly
115        // different for each volume.  If we have identical volumes at exactly
116        // the same position with exactly the same number of slices, the
117        // second volume will overwrite the first, so the first won't appear
118        // at all.
119        volumes.push_back(volPtr);
120        volPtr->numSlices(256 - volumes.size());
121    }
122
123    glPushAttrib(GL_ENABLE_BIT);
124
125    //two dimension pointer array
126    ConvexPolygon ***polys = new ConvexPolygon**[volumes.size()];
127    //number of actual slices for each volume
128    size_t *actual_slices = new size_t[volumes.size()];
129    float *z_steps = new float[volumes.size()];
130
131    TRACE("start loop %d", volumes.size());
132    for (size_t i = 0; i < volumes.size(); i++) {
133        Volume *volPtr = volumes[i];
134        polys[i] = NULL;
135        actual_slices[i] = 0;
136
137        int n_slices = volPtr->numSlices();
138        if (volPtr->isosurface()) {
139            // double the number of slices
140            n_slices <<= 1;
141        }
142
143        //volume start location
144        Vector3 volPos = volPtr->location();
145        Vector3 volScaling = volPtr->getPhysicalScaling();
146
147        TRACE("VOL POS: %g %g %g",
148              volPos.x, volPos.y, volPos.z);
149        TRACE("VOL SCALE: %g %g %g",
150              volScaling.x, volScaling.y, volScaling.z);
151
152        double x0 = 0;
153        double y0 = 0;
154        double z0 = 0;
155
156        Mat4x4 model_view_no_trans, model_view_trans;
157        Mat4x4 model_view_no_trans_inverse, model_view_trans_inverse;
158
159        //initialize volume plane with world coordinates
160        Plane volume_planes[6];
161        volume_planes[0].setCoeffs( 1,  0,  0, -x0);
162        volume_planes[1].setCoeffs(-1,  0,  0,  x0+1);
163        volume_planes[2].setCoeffs( 0,  1,  0, -y0);
164        volume_planes[3].setCoeffs( 0, -1,  0,  y0+1);
165        volume_planes[4].setCoeffs( 0,  0,  1, -z0);
166        volume_planes[5].setCoeffs( 0,  0, -1,  z0+1);
167
168        //get modelview matrix with no translation
169        glPushMatrix();
170        glScalef(volScaling.x, volScaling.y, volScaling.z);
171
172        glEnable(GL_DEPTH_TEST);
173
174        GLfloat mv_no_trans[16];
175        glGetFloatv(GL_MODELVIEW_MATRIX, mv_no_trans);
176
177        model_view_no_trans = Mat4x4(mv_no_trans);
178        model_view_no_trans_inverse = model_view_no_trans.inverse();
179
180        glPopMatrix();
181
182        //get modelview matrix with translation
183        glPushMatrix();
184        glTranslatef(volPos.x, volPos.y, volPos.z);
185        glScalef(volScaling.x, volScaling.y, volScaling.z);
186
187        GLfloat mv_trans[16];
188        glGetFloatv(GL_MODELVIEW_MATRIX, mv_trans);
189
190        model_view_trans = Mat4x4(mv_trans);
191        model_view_trans_inverse = model_view_trans.inverse();
192
193        //draw volume bounding box with translation (the correct location in
194        //space)
195        if (volPtr->outline()) {
196            float olcolor[3];
197            volPtr->getOutlineColor(olcolor);
198            drawBoundingBox(x0, y0, z0, x0+1, y0+1, z0+1,
199                (double)olcolor[0], (double)olcolor[1], (double)olcolor[2],
200                1.5);
201        }
202        glPopMatrix();
203
204        // transform volume_planes to eye coordinates.
205        // Need to transform without translation since we don't want
206        // to translate plane normals, just rotate them
207        for (size_t j = 0; j < 6; j++) {
208            volume_planes[j].transform(model_view_no_trans);
209        }
210        double eyeMinX, eyeMaxX, eyeMinY, eyeMaxY, zNear, zFar;
211        getEyeSpaceBounds(model_view_no_trans, 
212                          eyeMinX, eyeMaxX,
213                          eyeMinY, eyeMaxY,
214                          zNear, zFar);
215
216        //compute actual rendering slices
217        float z_step = fabs(zNear-zFar)/n_slices;
218        z_steps[i] = z_step;
219        size_t n_actual_slices;
220
221        if (volPtr->dataEnabled()) {
222            if (z_step == 0.0f)
223                n_actual_slices = 1;
224            else
225                n_actual_slices = (int)(fabs(zNear-zFar)/z_step + 1);
226            polys[i] = new ConvexPolygon*[n_actual_slices];
227        } else {
228            n_actual_slices = 0;
229            polys[i] = NULL;
230        }
231        actual_slices[i] = n_actual_slices;
232
233        TRACE("near: %g far: %g eye space bounds: (%g,%g)-(%g,%g) z_step: %g slices: %d",
234              zNear, zFar, eyeMinX, eyeMaxX, eyeMinY, eyeMaxY, z_step, n_actual_slices);
235
236        Vector4 vert1, vert2, vert3, vert4;
237
238        // Render cutplanes first with depth test enabled.  They will mark the
239        // image with their depth values.  Then we render other volume slices.
240        // These volume slices will be occluded correctly by the cutplanes and
241        // vice versa.
242
243        for (int j = 0; j < volPtr->getCutplaneCount(); j++) {
244            if (!volPtr->isCutplaneEnabled(j)) {
245                continue;
246            }
247            float offset = volPtr->getCutplane(j)->offset;
248            int axis = volPtr->getCutplane(j)->orient;
249
250            switch (axis) {
251            case 1:
252                vert1 = Vector4(offset, 0, 0, 1);
253                vert2 = Vector4(offset, 1, 0, 1);
254                vert3 = Vector4(offset, 1, 1, 1);
255                vert4 = Vector4(offset, 0, 1, 1);
256                break;
257            case 2:
258                vert1 = Vector4(0, offset, 0, 1);
259                vert2 = Vector4(1, offset, 0, 1);
260                vert3 = Vector4(1, offset, 1, 1);
261                vert4 = Vector4(0, offset, 1, 1);
262                break;
263            case 3:
264            default:
265                vert1 = Vector4(0, 0, offset, 1);
266                vert2 = Vector4(1, 0, offset, 1);
267                vert3 = Vector4(1, 1, offset, 1);
268                vert4 = Vector4(0, 1, offset, 1);
269                break;
270            }
271
272            Vector4 texcoord1 = vert1;
273            Vector4 texcoord2 = vert2;
274            Vector4 texcoord3 = vert3;
275            Vector4 texcoord4 = vert4;
276
277            _cutplaneShader->bind();
278            _cutplaneShader->setFPTextureParameter("volume", volPtr->textureID());
279            _cutplaneShader->setFPTextureParameter("tf", volPtr->transferFunction()->id());
280
281            glPushMatrix();
282            glTranslatef(volPos.x, volPos.y, volPos.z);
283            glScalef(volScaling.x, volScaling.y, volScaling.z);
284            _cutplaneShader->setGLStateMatrixVPParameter("modelViewProjMatrix",
285                                                         NvShader::MODELVIEW_PROJECTION_MATRIX);
286            glPopMatrix();
287
288            glEnable(GL_DEPTH_TEST);
289            glDisable(GL_BLEND);
290
291            glBegin(GL_QUADS);
292            glTexCoord3f(texcoord1.x, texcoord1.y, texcoord1.z);
293            glVertex3f(vert1.x, vert1.y, vert1.z);
294            glTexCoord3f(texcoord2.x, texcoord2.y, texcoord2.z);
295            glVertex3f(vert2.x, vert2.y, vert2.z);
296            glTexCoord3f(texcoord3.x, texcoord3.y, texcoord3.z);
297            glVertex3f(vert3.x, vert3.y, vert3.z);
298            glTexCoord3f(texcoord4.x, texcoord4.y, texcoord4.z);
299            glVertex3f(vert4.x, vert4.y, vert4.z);
300            glEnd();
301
302            glDisable(GL_DEPTH_TEST);
303            _cutplaneShader->disableFPTextureParameter("tf");
304            _cutplaneShader->disableFPTextureParameter("volume");
305            _cutplaneShader->unbind();
306        } //done cutplanes
307
308        // Now prepare proxy geometry slices
309
310        // Initialize view-aligned quads with eye space bounds of
311        // volume
312        vert1 = Vector4(eyeMinX, eyeMinY, -0.5, 1);
313        vert2 = Vector4(eyeMaxX, eyeMinY, -0.5, 1);
314        vert3 = Vector4(eyeMaxX, eyeMaxY, -0.5, 1);
315        vert4 = Vector4(eyeMinX, eyeMaxY, -0.5, 1);
316
317        size_t counter = 0;
318
319        // Transform slices and store them
320        float slice_z;
321        for (size_t j = 0; j < n_actual_slices; j++) {
322            slice_z = zFar + j * z_step; //back to front
323
324            ConvexPolygon *poly = new ConvexPolygon();
325            polys[i][counter] = poly;
326            counter++;
327
328            poly->vertices.clear();
329            poly->setId(i);
330
331            // Set eye space Z-coordinate of slice
332            vert1.z = slice_z;
333            vert2.z = slice_z;
334            vert3.z = slice_z;
335            vert4.z = slice_z;
336
337            poly->appendVertex(vert1);
338            poly->appendVertex(vert2);
339            poly->appendVertex(vert3);
340            poly->appendVertex(vert4);
341
342            for (size_t k = 0; k < 6; k++) {
343                if (!poly->clip(volume_planes[k], true))
344                    break;
345            }
346
347            if (poly->vertices.size() >= 3) {
348                poly->transform(model_view_no_trans_inverse);
349                poly->transform(model_view_trans);
350                total_rendered_slices++;
351            }
352        }
353    } //iterate all volumes
354    TRACE("end loop");
355
356    // We sort all the polygons according to their eye-space depth, from
357    // farthest to the closest.  This step is critical for correct blending
358
359    SortElement *slices = (SortElement *)
360        malloc(sizeof(SortElement) * total_rendered_slices);
361
362    size_t counter = 0;
363    for (size_t i = 0; i < volumes.size(); i++) {
364        for (size_t j = 0; j < actual_slices[i]; j++) {
365            if (polys[i][j]->vertices.size() >= 3) {
366                slices[counter] = SortElement(polys[i][j]->vertices[0].z, i, j);
367                counter++;
368            }
369        }
370    }
371
372    //sort them
373    qsort(slices, total_rendered_slices, sizeof(SortElement), sliceSort);
374
375    //Now we are ready to render all the slices from back to front
376    glEnable(GL_DEPTH_TEST);
377    // Non pre-multiplied alpha
378    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
379    glEnable(GL_BLEND);
380
381    for (size_t i = 0; i < total_rendered_slices; i++) {
382        Volume *volPtr = NULL;
383
384        int volume_index = slices[i].volumeId;
385        int slice_index = slices[i].sliceId;
386        ConvexPolygon *currentSlice = polys[volume_index][slice_index];
387        float z_step = z_steps[volume_index];
388
389        volPtr = volumes[volume_index];
390
391        Vector3 volScaling = volPtr->getPhysicalScaling();
392
393        glPushMatrix();
394        glScalef(volScaling.x, volScaling.y, volScaling.z);
395
396        // FIXME: compute view-dependent volume sample distance
397        double avgSampleDistance = 1.0 / pow(volPtr->width() * volScaling.x * 
398                                             volPtr->height() * volScaling.y * 
399                                             volPtr->depth() * volScaling.z, 1.0/3.0);
400        float sampleRatio = z_step / avgSampleDistance;
401
402#ifdef notdef
403        TRACE("shading slice: volume %s addr=%x slice=%d, volume=%d z_step=%g avgSD=%g", 
404              volPtr->name(), volPtr, slice_index, volume_index, z_step, avgSampleDistance);
405#endif
406        activateVolumeShader(volPtr, false, sampleRatio);
407        glPopMatrix();
408
409        glBegin(GL_POLYGON);
410        currentSlice->emit(true);
411        glEnd();
412
413        deactivateVolumeShader();
414    }
415
416    glPopAttrib();
417
418    //Deallocate all the memory used
419    for (size_t i = 0; i < volumes.size(); i++) {
420        for (size_t j = 0; j <actual_slices[i]; j++) {
421            delete polys[i][j];
422        }
423        if (polys[i]) {
424            delete[] polys[i];
425        }
426    }
427    delete[] polys;
428    delete[] actual_slices;
429    delete[] z_steps;
430    free(slices);
431}
432
433void 
434VolumeRenderer::drawBoundingBox(float x0, float y0, float z0, 
435                                float x1, float y1, float z1,
436                                float r, float g, float b, 
437                                float line_width)
438{
439    glPushAttrib(GL_ENABLE_BIT);
440
441    glEnable(GL_DEPTH_TEST);
442    glDisable(GL_TEXTURE_2D);
443    glEnable(GL_BLEND);
444
445    glMatrixMode(GL_MODELVIEW);
446    glPushMatrix();
447
448    glColor4d(r, g, b, 1.0);
449    glLineWidth(line_width);
450
451    glBegin(GL_LINE_LOOP);
452    {
453        glVertex3d(x0, y0, z0);
454        glVertex3d(x1, y0, z0);
455        glVertex3d(x1, y1, z0);
456        glVertex3d(x0, y1, z0);
457    }
458    glEnd();
459
460    glBegin(GL_LINE_LOOP);
461    {
462        glVertex3d(x0, y0, z1);
463        glVertex3d(x1, y0, z1);
464        glVertex3d(x1, y1, z1);
465        glVertex3d(x0, y1, z1);
466    }
467    glEnd();
468
469    glBegin(GL_LINE_LOOP);
470    {
471        glVertex3d(x0, y0, z0);
472        glVertex3d(x0, y0, z1);
473        glVertex3d(x0, y1, z1);
474        glVertex3d(x0, y1, z0);
475    }
476    glEnd();
477
478    glBegin(GL_LINE_LOOP);
479    {
480        glVertex3d(x1, y0, z0);
481        glVertex3d(x1, y0, z1);
482        glVertex3d(x1, y1, z1);
483        glVertex3d(x1, y1, z0);
484    }
485    glEnd();
486
487    glPopMatrix();
488    glPopAttrib();
489}
490
491void 
492VolumeRenderer::activateVolumeShader(Volume *volPtr, bool sliceMode,
493                                     float sampleRatio)
494{
495    //vertex shader
496    _stdVertexShader->bind();
497    TransferFunction *tfPtr  = volPtr->transferFunction();
498    if (volPtr->volumeType() == Volume::CUBIC) {
499        _regularVolumeShader->bind(tfPtr->id(), volPtr, sliceMode, sampleRatio);
500    } else if (volPtr->volumeType() == Volume::ZINCBLENDE) {
501        _zincBlendeShader->bind(tfPtr->id(), volPtr, sliceMode, sampleRatio);
502    }
503}
504
505void VolumeRenderer::deactivateVolumeShader()
506{
507    _stdVertexShader->unbind();
508    _regularVolumeShader->unbind();
509    _zincBlendeShader->unbind();
510}
511
512void VolumeRenderer::getEyeSpaceBounds(const Mat4x4& mv,
513                                       double& xMin, double& xMax,
514                                       double& yMin, double& yMax,
515                                       double& zNear, double& zFar)
516{
517    double x0 = 0;
518    double y0 = 0;
519    double z0 = 0;
520    double x1 = 1;
521    double y1 = 1;
522    double z1 = 1;
523
524    double zMin, zMax;
525    xMin = DBL_MAX;
526    xMax = -DBL_MAX;
527    yMin = DBL_MAX;
528    yMax = -DBL_MAX;
529    zMin = DBL_MAX;
530    zMax = -DBL_MAX;
531
532    double vertex[8][4];
533
534    vertex[0][0]=x0; vertex[0][1]=y0; vertex[0][2]=z0; vertex[0][3]=1.0;
535    vertex[1][0]=x1; vertex[1][1]=y0; vertex[1][2]=z0; vertex[1][3]=1.0;
536    vertex[2][0]=x0; vertex[2][1]=y1; vertex[2][2]=z0; vertex[2][3]=1.0;
537    vertex[3][0]=x0; vertex[3][1]=y0; vertex[3][2]=z1; vertex[3][3]=1.0;
538    vertex[4][0]=x1; vertex[4][1]=y1; vertex[4][2]=z0; vertex[4][3]=1.0;
539    vertex[5][0]=x1; vertex[5][1]=y0; vertex[5][2]=z1; vertex[5][3]=1.0;
540    vertex[6][0]=x0; vertex[6][1]=y1; vertex[6][2]=z1; vertex[6][3]=1.0;
541    vertex[7][0]=x1; vertex[7][1]=y1; vertex[7][2]=z1; vertex[7][3]=1.0;
542
543    for (int i = 0; i < 8; i++) {
544        Vector4 eyeVert = mv.transform(Vector4(vertex[i][0],
545                                               vertex[i][1],
546                                               vertex[i][2],
547                                               vertex[i][3]));
548        if (eyeVert.x < xMin) xMin = eyeVert.x;
549        if (eyeVert.x > xMax) xMax = eyeVert.x;
550        if (eyeVert.y < yMin) yMin = eyeVert.y;
551        if (eyeVert.y > yMax) yMax = eyeVert.y;
552        if (eyeVert.z < zMin) zMin = eyeVert.z;
553        if (eyeVert.z > zMax) zMax = eyeVert.z;
554    }
555
556    zNear = zMax;
557    zFar = zMin;
558}
Note: See TracBrowser for help on using the repository browser.