source: nanovis/branches/1.1/VolumeRenderer.cpp @ 4818

Last change on this file since 4818 was 4818, checked in by ldelgass, 5 years ago

Merge in 'cutplane visible' command

  • Property svn:eol-style set to native
File size: 18.1 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 *volume = volumes[i];
135        polys[i] = NULL;
136        actual_slices[i] = 0;
137
138        int n_slices = volume->numSlices();
139        if (volume->isosurface()) {
140            // double the number of slices
141            n_slices <<= 1;
142        }
143
144        //volume start location
145        Vector3f volPos = volume->location();
146        Vector3f volScaling = volume->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 (volume->outline()) {
199            float olcolor[3];
200            volume->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 (volume->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 < volume->getCutplaneCount(); j++) {
247            if (!volume->cutplanesVisible() || !volume->isCutplaneEnabled(j)) {
248                continue;
249            }
250            float offset = volume->getCutplane(j)->offset;
251            int axis = volume->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", volume->textureID());
282            _cutplaneShader->setFPTextureParameter("tf", volume->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 *volume = 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        volume = volumes[volume_index];
393
394        Vector3f volScaling = volume->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(volume->width() * volScaling.x * 
401                                             volume->height() * volScaling.y * 
402                                             volume->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              volume->name(), volume, slice_index, volume_index, z_step, avgSampleDistance);
408#endif
409        activateVolumeShader(volume, 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 *volume, bool sliceMode,
496                                     float sampleRatio)
497{
498    //vertex shader
499    _stdVertexShader->bind();
500    TransferFunction *transferFunc  = volume->transferFunction();
501    if (volume->volumeType() == Volume::CUBIC) {
502        _regularVolumeShader->bind(transferFunc->id(), volume, sliceMode, sampleRatio);
503    } else if (volume->volumeType() == Volume::ZINCBLENDE) {
504        _zincBlendeShader->bind(transferFunc->id(), volume, 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.