source: nanovis/branches/1.2/VolumeRenderer.cpp @ 6632

Last change on this file since 6632 was 4904, checked in by ldelgass, 9 years ago

Merge serveral changes from trunk. Does not include threading, world space
changes, etc.

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