source: nanovis/trunk/VolumeRenderer.cpp @ 5064

Last change on this file since 5064 was 3883, checked in by ldelgass, 11 years ago

Add cutplane visibility flag to nanovis. Will be used to add cutplanes button
in client

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