source: trunk/packages/vizservers/nanovis/VolumeRenderer.cpp @ 3502

Last change on this file since 3502 was 3502, checked in by ldelgass, 7 years ago

Add basic VTK structured points reader to nanovis, update copyright dates.

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