source: nanovis/trunk/VolumeRenderer.cpp @ 4794

Last change on this file since 4794 was 3883, checked in by ldelgass, 6 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.