source: branches/1.2/packages/vizservers/nanovis/VolumeRenderer.cpp @ 3589

Last change on this file since 3589 was 3567, checked in by ldelgass, 7 years ago

Refactor and cleanups in nanovis, mainly to switch to using STL hash tables
(TR1 required) instead of Tcl hash tables, split out Flow particles and boxes
to separate implementation files. The goal is to achieve better separation of
Tcl command parsing and the core graphics rendering objects and code.

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