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

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

Merge r3611:3618 from trunk

  • 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 "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", "main");
57    _cutplaneShader->loadFragmentProgram("cutplane_fp.cg", "main");
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 a 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        //volume start location
146        Vector3f volPos = volume->location();
147        Vector3f volScaling = volume->getPhysicalScaling();
148
149        TRACE("VOL POS: %g %g %g",
150              volPos.x, volPos.y, volPos.z);
151        TRACE("VOL SCALE: %g %g %g",
152              volScaling.x, volScaling.y, volScaling.z);
153
154        double x0 = 0;
155        double y0 = 0;
156        double z0 = 0;
157
158        Matrix4x4d model_view_no_trans, model_view_trans;
159        Matrix4x4d model_view_no_trans_inverse, model_view_trans_inverse;
160
161        //initialize volume plane with world coordinates
162        nv::Plane volume_planes[6];
163        volume_planes[0].setCoeffs( 1,  0,  0, -x0);
164        volume_planes[1].setCoeffs(-1,  0,  0,  x0+1);
165        volume_planes[2].setCoeffs( 0,  1,  0, -y0);
166        volume_planes[3].setCoeffs( 0, -1,  0,  y0+1);
167        volume_planes[4].setCoeffs( 0,  0,  1, -z0);
168        volume_planes[5].setCoeffs( 0,  0, -1,  z0+1);
169
170        //get modelview matrix with no translation
171        glPushMatrix();
172        glScalef(volScaling.x, volScaling.y, volScaling.z);
173
174        glEnable(GL_DEPTH_TEST);
175
176        GLdouble mv_no_trans[16];
177        glGetDoublev(GL_MODELVIEW_MATRIX, mv_no_trans);
178
179        model_view_no_trans = Matrix4x4d(mv_no_trans);
180        model_view_no_trans_inverse = model_view_no_trans.inverse();
181
182        glPopMatrix();
183
184        //get modelview matrix with translation
185        glPushMatrix();
186        glTranslatef(volPos.x, volPos.y, volPos.z);
187        glScalef(volScaling.x, volScaling.y, volScaling.z);
188
189        GLdouble mv_trans[16];
190        glGetDoublev(GL_MODELVIEW_MATRIX, mv_trans);
191
192        model_view_trans = Matrix4x4d(mv_trans);
193        model_view_trans_inverse = model_view_trans.inverse();
194
195        model_view_trans.print();
196
197        //draw volume bounding box with translation (the correct location in
198        //space)
199        if (volume->outline()) {
200            float olcolor[3];
201            volume->getOutlineColor(olcolor);
202            drawBoundingBox(x0, y0, z0, x0+1, y0+1, z0+1,
203                (double)olcolor[0], (double)olcolor[1], (double)olcolor[2],
204                1.5);
205        }
206        glPopMatrix();
207
208        // transform volume_planes to eye coordinates.
209        // Need to transform without translation since we don't want
210        // to translate plane normals, just rotate them
211        for (size_t j = 0; j < 6; j++) {
212            volume_planes[j].transform(model_view_no_trans);
213        }
214        double eyeMinX, eyeMaxX, eyeMinY, eyeMaxY, zNear, zFar;
215        getEyeSpaceBounds(model_view_no_trans, 
216                          eyeMinX, eyeMaxX,
217                          eyeMinY, eyeMaxY,
218                          zNear, zFar);
219
220        //compute actual rendering slices
221        float z_step = fabs(zNear-zFar)/n_slices;
222        z_steps[i] = z_step;
223        size_t n_actual_slices;
224
225        if (volume->dataEnabled()) {
226            if (z_step == 0.0f)
227                n_actual_slices = 1;
228            else
229                n_actual_slices = (int)(fabs(zNear-zFar)/z_step + 1);
230            polys[i] = new ConvexPolygon*[n_actual_slices];
231        } else {
232            n_actual_slices = 0;
233            polys[i] = NULL;
234        }
235        actual_slices[i] = n_actual_slices;
236
237        TRACE("near: %g far: %g eye space bounds: (%g,%g)-(%g,%g) z_step: %g slices: %d",
238              zNear, zFar, eyeMinX, eyeMaxX, eyeMinY, eyeMaxY, z_step, n_actual_slices);
239
240        Vector4f vert1, vert2, vert3, vert4;
241
242        // Render cutplanes first with depth test enabled.  They will mark the
243        // image with their depth values.  Then we render other volume slices.
244        // These volume slices will be occluded correctly by the cutplanes and
245        // vice versa.
246
247        for (int j = 0; j < volume->getCutplaneCount(); j++) {
248            if (!volume->cutplanesVisible() || !volume->isCutplaneEnabled(j)) {
249                continue;
250            }
251            float offset = volume->getCutplane(j)->offset;
252            int axis = volume->getCutplane(j)->orient;
253
254            switch (axis) {
255            case 1:
256                vert1 = Vector4f(offset, 0, 0, 1);
257                vert2 = Vector4f(offset, 1, 0, 1);
258                vert3 = Vector4f(offset, 1, 1, 1);
259                vert4 = Vector4f(offset, 0, 1, 1);
260                break;
261            case 2:
262                vert1 = Vector4f(0, offset, 0, 1);
263                vert2 = Vector4f(1, offset, 0, 1);
264                vert3 = Vector4f(1, offset, 1, 1);
265                vert4 = Vector4f(0, offset, 1, 1);
266                break;
267            case 3:
268            default:
269                vert1 = Vector4f(0, 0, offset, 1);
270                vert2 = Vector4f(1, 0, offset, 1);
271                vert3 = Vector4f(1, 1, offset, 1);
272                vert4 = Vector4f(0, 1, offset, 1);
273                break;
274            }
275
276            Vector4f texcoord1 = vert1;
277            Vector4f texcoord2 = vert2;
278            Vector4f texcoord3 = vert3;
279            Vector4f texcoord4 = vert4;
280
281            _cutplaneShader->bind();
282            _cutplaneShader->setFPTextureParameter("volume", volume->textureID());
283            _cutplaneShader->setFPTextureParameter("tf", volume->transferFunction()->id());
284
285            glPushMatrix();
286            glTranslatef(volPos.x, volPos.y, volPos.z);
287            glScalef(volScaling.x, volScaling.y, volScaling.z);
288            _cutplaneShader->setGLStateMatrixVPParameter("modelViewProjMatrix",
289                                                         Shader::MODELVIEW_PROJECTION_MATRIX);
290            glPopMatrix();
291
292            glEnable(GL_DEPTH_TEST);
293            glDisable(GL_BLEND);
294
295            glBegin(GL_QUADS);
296            glTexCoord3f(texcoord1.x, texcoord1.y, texcoord1.z);
297            glVertex3f(vert1.x, vert1.y, vert1.z);
298            glTexCoord3f(texcoord2.x, texcoord2.y, texcoord2.z);
299            glVertex3f(vert2.x, vert2.y, vert2.z);
300            glTexCoord3f(texcoord3.x, texcoord3.y, texcoord3.z);
301            glVertex3f(vert3.x, vert3.y, vert3.z);
302            glTexCoord3f(texcoord4.x, texcoord4.y, texcoord4.z);
303            glVertex3f(vert4.x, vert4.y, vert4.z);
304            glEnd();
305
306            glDisable(GL_DEPTH_TEST);
307            _cutplaneShader->disableFPTextureParameter("tf");
308            _cutplaneShader->disableFPTextureParameter("volume");
309            _cutplaneShader->unbind();
310        } //done cutplanes
311
312        // Now prepare proxy geometry slices
313
314        // Initialize view-aligned quads with eye space bounds of
315        // volume
316        vert1 = Vector4f(eyeMinX, eyeMinY, -0.5, 1);
317        vert2 = Vector4f(eyeMaxX, eyeMinY, -0.5, 1);
318        vert3 = Vector4f(eyeMaxX, eyeMaxY, -0.5, 1);
319        vert4 = Vector4f(eyeMinX, eyeMaxY, -0.5, 1);
320
321        size_t counter = 0;
322
323        // Transform slices and store them
324        float slice_z;
325        for (size_t j = 0; j < n_actual_slices; j++) {
326            slice_z = zFar + j * z_step; //back to front
327
328            ConvexPolygon *poly = new ConvexPolygon();
329            polys[i][counter] = poly;
330            counter++;
331
332            poly->vertices.clear();
333            poly->setId(i);
334
335            // Set eye space Z-coordinate of slice
336            vert1.z = slice_z;
337            vert2.z = slice_z;
338            vert3.z = slice_z;
339            vert4.z = slice_z;
340
341            poly->appendVertex(vert1);
342            poly->appendVertex(vert2);
343            poly->appendVertex(vert3);
344            poly->appendVertex(vert4);
345
346            for (size_t k = 0; k < 6; k++) {
347                if (!poly->clip(volume_planes[k], true))
348                    break;
349            }
350
351            if (poly->vertices.size() >= 3) {
352                poly->transform(model_view_no_trans_inverse);
353                poly->transform(model_view_trans);
354                total_rendered_slices++;
355            }
356        }
357    } //iterate all volumes
358    TRACE("end loop");
359
360    // We sort all the polygons according to their eye-space depth, from
361    // farthest to the closest.  This step is critical for correct blending
362
363    SortElement *slices = (SortElement *)
364        malloc(sizeof(SortElement) * total_rendered_slices);
365
366    size_t counter = 0;
367    for (size_t i = 0; i < volumes.size(); i++) {
368        for (size_t j = 0; j < actual_slices[i]; j++) {
369            if (polys[i][j]->vertices.size() >= 3) {
370                slices[counter] = SortElement(polys[i][j]->vertices[0].z, i, j);
371                counter++;
372            }
373        }
374    }
375
376    //sort them
377    qsort(slices, total_rendered_slices, sizeof(SortElement), sliceSort);
378
379    //Now we are ready to render all the slices from back to front
380    glEnable(GL_DEPTH_TEST);
381    // Non pre-multiplied alpha
382    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
383    glEnable(GL_BLEND);
384
385    for (size_t i = 0; i < total_rendered_slices; i++) {
386        Volume *volume = NULL;
387
388        int volume_index = slices[i].volumeId;
389        int slice_index = slices[i].sliceId;
390        ConvexPolygon *currentSlice = polys[volume_index][slice_index];
391        float z_step = z_steps[volume_index];
392
393        volume = volumes[volume_index];
394
395        Vector3f volScaling = volume->getPhysicalScaling();
396
397        glPushMatrix();
398        glScalef(volScaling.x, volScaling.y, volScaling.z);
399
400        // FIXME: compute view-dependent volume sample distance
401        double avgSampleDistance = 1.0 / pow(volume->width() * volScaling.x * 
402                                             volume->height() * volScaling.y * 
403                                             volume->depth() * volScaling.z, 1.0/3.0);
404        float sampleRatio = z_step / avgSampleDistance;
405
406#ifdef notdef
407        TRACE("shading slice: volume %s addr=%x slice=%d, volume=%d z_step=%g avgSD=%g", 
408              volume->name(), volume, slice_index, volume_index, z_step, avgSampleDistance);
409#endif
410        activateVolumeShader(volume, false, sampleRatio);
411        glPopMatrix();
412
413        glBegin(GL_POLYGON);
414        currentSlice->emit(true);
415        glEnd();
416
417        deactivateVolumeShader();
418    }
419
420    glPopAttrib();
421
422    //Deallocate all the memory used
423    for (size_t i = 0; i < volumes.size(); i++) {
424        for (size_t j = 0; j <actual_slices[i]; j++) {
425            delete polys[i][j];
426        }
427        if (polys[i]) {
428            delete[] polys[i];
429        }
430    }
431    delete[] polys;
432    delete[] actual_slices;
433    delete[] z_steps;
434    free(slices);
435}
436
437void 
438VolumeRenderer::drawBoundingBox(float x0, float y0, float z0, 
439                                float x1, float y1, float z1,
440                                float r, float g, float b, 
441                                float line_width)
442{
443    glPushAttrib(GL_ENABLE_BIT);
444
445    glEnable(GL_DEPTH_TEST);
446    glDisable(GL_TEXTURE_2D);
447    glEnable(GL_BLEND);
448
449    glMatrixMode(GL_MODELVIEW);
450    glPushMatrix();
451
452    glColor4d(r, g, b, 1.0);
453    glLineWidth(line_width);
454
455    glBegin(GL_LINE_LOOP);
456    {
457        glVertex3d(x0, y0, z0);
458        glVertex3d(x1, y0, z0);
459        glVertex3d(x1, y1, z0);
460        glVertex3d(x0, y1, z0);
461    }
462    glEnd();
463
464    glBegin(GL_LINE_LOOP);
465    {
466        glVertex3d(x0, y0, z1);
467        glVertex3d(x1, y0, z1);
468        glVertex3d(x1, y1, z1);
469        glVertex3d(x0, y1, z1);
470    }
471    glEnd();
472
473    glBegin(GL_LINE_LOOP);
474    {
475        glVertex3d(x0, y0, z0);
476        glVertex3d(x0, y0, z1);
477        glVertex3d(x0, y1, z1);
478        glVertex3d(x0, y1, z0);
479    }
480    glEnd();
481
482    glBegin(GL_LINE_LOOP);
483    {
484        glVertex3d(x1, y0, z0);
485        glVertex3d(x1, y0, z1);
486        glVertex3d(x1, y1, z1);
487        glVertex3d(x1, y1, z0);
488    }
489    glEnd();
490
491    glPopMatrix();
492    glPopAttrib();
493}
494
495void 
496VolumeRenderer::activateVolumeShader(Volume *volume, bool sliceMode,
497                                     float sampleRatio)
498{
499    //vertex shader
500    _stdVertexShader->bind();
501    TransferFunction *transferFunc  = volume->transferFunction();
502    if (volume->volumeType() == Volume::CUBIC) {
503        _regularVolumeShader->bind(transferFunc->id(), volume, sliceMode, sampleRatio);
504    } else if (volume->volumeType() == Volume::ZINCBLENDE) {
505        _zincBlendeShader->bind(transferFunc->id(), volume, sliceMode, sampleRatio);
506    }
507}
508
509void VolumeRenderer::deactivateVolumeShader()
510{
511    _stdVertexShader->unbind();
512    _regularVolumeShader->unbind();
513    _zincBlendeShader->unbind();
514}
515
516void VolumeRenderer::getEyeSpaceBounds(const Matrix4x4d& mv,
517                                       double& xMin, double& xMax,
518                                       double& yMin, double& yMax,
519                                       double& zNear, double& zFar)
520{
521    double x0 = 0;
522    double y0 = 0;
523    double z0 = 0;
524    double x1 = 1;
525    double y1 = 1;
526    double z1 = 1;
527
528    double zMin, zMax;
529    xMin = DBL_MAX;
530    xMax = -DBL_MAX;
531    yMin = DBL_MAX;
532    yMax = -DBL_MAX;
533    zMin = DBL_MAX;
534    zMax = -DBL_MAX;
535
536    double vertex[8][4];
537
538    vertex[0][0]=x0; vertex[0][1]=y0; vertex[0][2]=z0; vertex[0][3]=1.0;
539    vertex[1][0]=x1; vertex[1][1]=y0; vertex[1][2]=z0; vertex[1][3]=1.0;
540    vertex[2][0]=x0; vertex[2][1]=y1; vertex[2][2]=z0; vertex[2][3]=1.0;
541    vertex[3][0]=x0; vertex[3][1]=y0; vertex[3][2]=z1; vertex[3][3]=1.0;
542    vertex[4][0]=x1; vertex[4][1]=y1; vertex[4][2]=z0; vertex[4][3]=1.0;
543    vertex[5][0]=x1; vertex[5][1]=y0; vertex[5][2]=z1; vertex[5][3]=1.0;
544    vertex[6][0]=x0; vertex[6][1]=y1; vertex[6][2]=z1; vertex[6][3]=1.0;
545    vertex[7][0]=x1; vertex[7][1]=y1; vertex[7][2]=z1; vertex[7][3]=1.0;
546
547    for (int i = 0; i < 8; i++) {
548        Vector4f eyeVert = mv.transform(Vector4f(vertex[i][0],
549                                                 vertex[i][1],
550                                                 vertex[i][2],
551                                                 vertex[i][3]));
552        if (eyeVert.x < xMin) xMin = eyeVert.x;
553        if (eyeVert.x > xMax) xMax = eyeVert.x;
554        if (eyeVert.y < yMin) yMin = eyeVert.y;
555        if (eyeVert.y > yMax) yMax = eyeVert.y;
556        if (eyeVert.z < zMin) zMin = eyeVert.z;
557        if (eyeVert.z > zMax) zMax = eyeVert.z;
558    }
559
560    zNear = zMax;
561    zFar = zMin;
562}
Note: See TracBrowser for help on using the repository browser.