source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/ShadowMapping.cpp @ 3258

Revision 3258, 17.6 KB checked in by mattausch, 16 years ago (diff)

worked on new method

Line 
1#include "ShadowMapping.h"
2#include "FrameBufferObject.h"
3#include "RenderState.h"
4#include "RenderTraverser.h"
5#include "Light.h"
6#include "Polygon3.h"
7#include "Polyhedron.h"
8#include "ResourceManager.h"
9#include "Camera.h"
10
11#include <IL/il.h>
12#include <assert.h>
13
14
15using namespace std;
16
17
18namespace CHCDemoEngine
19{
20
21static Polyhedron *polyhedron = NULL;
22static Polyhedron *lightPoly = NULL;
23
24
25static void PrintGLerror(char *msg)
26{
27        GLenum errCode;
28        const GLubyte *errStr;
29       
30        if ((errCode = glGetError()) != GL_NO_ERROR)
31        {
32                errStr = gluErrorString(errCode);
33                fprintf(stderr,"OpenGL ERROR: %s: %s\n", errStr, msg);
34        }
35}
36
37
38static Polyhedron *CreatePolyhedron(const Matrix4x4 &lightMatrix,
39                                                                        const AxisAlignedBox3 &sceneBox)
40{
41        Frustum frustum(lightMatrix);
42
43        vector<Plane3> clipPlanes;
44
45        for (int i = 0; i < 6; ++ i)
46        {
47                ////////////
48                //-- normalize the coefficients
49
50                // the clipping planes look outward the frustum,
51                // so distances > 0 mean that a point is outside
52                const float invLength = -1.0f / Magnitude(frustum.mClipPlanes[i].mNormal);
53
54                frustum.mClipPlanes[i].mD *= invLength;
55                frustum.mClipPlanes[i].mNormal *= invLength;
56        }
57
58        // first create near plane because of precision issues
59        clipPlanes.push_back(frustum.mClipPlanes[4]);
60
61        clipPlanes.push_back(frustum.mClipPlanes[0]);
62        clipPlanes.push_back(frustum.mClipPlanes[1]);
63        clipPlanes.push_back(frustum.mClipPlanes[2]);
64        clipPlanes.push_back(frustum.mClipPlanes[3]);
65        clipPlanes.push_back(frustum.mClipPlanes[5]);
66
67        return Polyhedron::CreatePolyhedron(clipPlanes, sceneBox);
68}
69
70
71static void GrabDepthBuffer(float *data, GLuint depthTexture)
72{
73        glEnable(GL_TEXTURE_2D);
74        glBindTexture(GL_TEXTURE_2D, depthTexture);
75
76        glGetTexImage(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT, GL_FLOAT, data);
77
78        glBindTexture(GL_TEXTURE_2D, 0);
79        glDisable(GL_TEXTURE_2D);
80}
81
82
83static void ExportDepthBuffer(float *data, int size)
84{
85        ilInit();
86        assert(ilGetError() == IL_NO_ERROR);
87
88        ILstring filename = ILstring("shadow.tga");
89        ilRegisterType(IL_FLOAT);
90
91        const int depth = 1;
92        const int bpp = 1;
93
94        if (!ilTexImage(size, size, depth, bpp, IL_LUMINANCE, IL_FLOAT, data))
95        {
96                cerr << "IL error " << ilGetError() << endl;
97       
98                ilShutDown();
99                assert(ilGetError() == IL_NO_ERROR);
100
101                return;
102        }
103
104        ilEnable(IL_FILE_OVERWRITE);
105        if (!ilSaveImage(filename))
106        {
107                cerr << "TGA write error " << ilGetError() << endl;
108        }
109
110        ilShutDown();
111        assert(ilGetError() == IL_NO_ERROR);
112
113        cout << "exported depth buffer" << endl;
114}
115
116
117
118static AxisAlignedBox3 GetExtremalPoints(const Matrix4x4 &m,
119                                                                                 const VertexArray &pts)
120{
121        AxisAlignedBox3 extremalPoints;
122        extremalPoints.Initialize();
123
124        VertexArray::const_iterator it, it_end = pts.end();
125               
126        for (it = pts.begin(); it != it_end; ++ it)
127        {
128                Vector3 pt = *it;
129                pt = m * pt;
130
131                extremalPoints.Include(pt);
132        }
133
134        return extremalPoints;
135}
136
137
138ShadowMap::ShadowMap(DirectionalLight *light,
139                                         int size,
140                                         const AxisAlignedBox3 &sceneBox,
141                                         PerspectiveCamera *cam):
142mSceneBox(sceneBox), mSize(size), mCamera(cam), mLight(light)
143{
144        mFbo = new FrameBufferObject(size, size, FrameBufferObject::DEPTH_32, true);
145
146        // need a color buffer to keep driver happy
147        mFbo->AddColorBuffer(ColorBufferObject::RGB_UBYTE,
148                                 ColorBufferObject::WRAP_CLAMP_TO_EDGE,
149                                                 ColorBufferObject::FILTER_NEAREST);
150
151        mShadowCam = new PerspectiveCamera(1);
152}
153
154
155ShadowMap::~ShadowMap()
156{
157        DEL_PTR(mFbo);
158        DEL_PTR(mShadowCam);
159
160        DEL_PTR(lightPoly);
161        DEL_PTR(polyhedron);
162}
163
164
165static void DrawPolyhedron(Polyhedron *poly, const Vector3 &color)
166{
167        if (!poly) return;
168
169        for (size_t i = 0; i < poly->NumPolygons(); ++ i)
170        {
171                glColor3f(color.x, color.y, color.z);
172
173                glBegin(GL_LINE_LOOP);
174
175                Polygon3 *p = poly->GetPolygons()[i];
176
177                for (size_t j = 0; j < p->mVertices.size(); ++ j)
178                {
179                        Vector3 v = p->mVertices[j];
180                        glVertex3d(v.x, v.y, v.z);
181                }
182
183                glEnd();
184        }
185}
186
187
188void ShadowMap::VisualizeFrustra()
189{
190        DrawPolyhedron(lightPoly, Vector3(1, 0, 1));
191        DrawPolyhedron(polyhedron, Vector3(0, 1, 0));
192}
193
194
195// z0 is the point that lies on the parallel plane to the near plane through e (A)
196//and on the near plane of the C frustum (the plane z = bZmax) and on the line x = e.x
197Vector3 ShadowMap::GetLightSpaceZ0(const Matrix4x4 &lightSpace,
198                                                                   const Vector3 &e,
199                                                                   const float maxZ,
200                                                                   const Vector3 &eyeDir) const
201{
202        // to calculate the parallel plane to the near plane through e we
203        // calculate the plane A with the three points
204        Plane3 planeA(e, eyeDir);
205
206        planeA.Transform(lightSpace);
207       
208        // get the parameters of A from the plane equation n dot d = 0
209        const float d = planeA.mD;
210        const Vector3 n = planeA.mNormal;
211       
212        // transform to light space
213        const Vector3 e_ls = lightSpace * e;
214
215        Vector3 z0;
216
217        z0.x = e_ls.x;
218        z0.y = (d - n.z * maxZ - n.x * e_ls.x) / n.y;
219        z0.z = maxZ;
220
221        return z0;
222        //return V3(e_ls.x(),(d-n.z()*b_lsZmax-n.x()*e_ls.x())/n.y(),b_lsZmax);
223}
224
225
226float ShadowMap::ComputeNOpt(const Matrix4x4 &lightSpace,
227                                                         const AxisAlignedBox3 &extremalPoints,
228                                                         const VertexArray &body) const
229{
230        const Vector3 nearPt = GetNearCameraPointE(body);
231        const Vector3 eyeDir = mCamera->GetDirection();
232
233        Matrix4x4 eyeView;
234        mCamera->GetModelViewMatrix(eyeView);
235
236        const Matrix4x4 invLightSpace = Invert(lightSpace);
237
238        const Vector3 z0_ls = GetLightSpaceZ0(lightSpace, nearPt, extremalPoints.Max().z, eyeDir);
239        const Vector3 z1_ls = Vector3(z0_ls.x, z0_ls.y, extremalPoints.Min().z);
240       
241        // transform back to world space
242        const Vector3 z0_ws = invLightSpace * z0_ls;
243        const Vector3 z1_ws = invLightSpace * z1_ls;
244
245        // transform to eye space
246        const Vector3 z0_es = eyeView * z0_ws;
247        const Vector3 z1_es = eyeView * z1_ws;
248
249        const float z0 = z0_es.z;
250        const float z1 = z1_es.z;
251
252        cout << "z0 ls: " << z0_ls << " z1 ls: " << z1_ls << endl;
253        cout << "z0: " << z0_es << " z1: " << z1_es << endl;
254
255        const float d = fabs(extremalPoints.Max()[2] - extremalPoints.Min()[2]);
256
257        const float n = d / (sqrt(z1 / z0) - 1.0f);
258
259        return n;
260}
261
262
263float ShadowMap::ComputeN(const AxisAlignedBox3 &extremalPoints) const
264{
265        const float nearPlane = mCamera->GetNear();
266       
267        const float d = fabs(extremalPoints.Max()[2] - extremalPoints.Min()[2]);
268       
269        const float dotProd = DotProd(mCamera->GetDirection(), mShadowCam->GetDirection());
270        const float sinGamma = sin(fabs(acos(dotProd)));
271
272        // test for values close to zero
273        if (sinGamma < 1e-6f) return 1e6f;
274       
275        const float scale = 2.0f;
276        return scale * (nearPlane + sqrt(nearPlane * (nearPlane + d * sinGamma))) /  sinGamma;
277}
278
279
280Matrix4x4 ShadowMap::CalcLispSMTransform(const Matrix4x4 &lightSpace,
281                                                                                 const AxisAlignedBox3 &extremalPoints,
282                                                                                 const VertexArray &body
283                                                                                 )
284{
285        AxisAlignedBox3 bounds_ls = GetExtremalPoints(lightSpace, body);
286
287        ///////////////
288        //-- We apply the lispsm algorithm in order to calculate an optimal light projection matrix
289        //-- first find the free parameter values n, and P (the projection center), and the projection depth
290
291        const float n = ComputeN(bounds_ls);
292        //const float n = ComputeNOpt(lightSpace, extremalPoints, body); cout << "n: " << n << endl;
293
294        if (n >= 1e6f) // light direction nearly parallel to view => switch to uniform
295                return IdentityMatrix();
296
297        const Vector3 nearPt = GetNearCameraPointE(body);
298       
299        //get the coordinates of the near camera point in light space
300        const Vector3 lsNear = lightSpace * nearPt;
301
302        // the start point has the x and y coordinate of e, the z coord of the near plane of the light volume
303        const Vector3 startPt = Vector3(lsNear.x, lsNear.y, bounds_ls.Max().z);
304       
305        // the new projection center
306        const Vector3 projCenter = startPt + Vector3::UNIT_Z() * n;
307
308        //construct a translation that moves to the projection center
309        const Matrix4x4 projectionCenter = TranslationMatrix(-projCenter);
310
311        // light space y size
312        const float d = fabs(bounds_ls.Max()[2] - bounds_ls.Min()[2]);
313
314        const float dy = fabs(bounds_ls.Max()[1] - bounds_ls.Min()[1]);
315        const float dx = fabs(bounds_ls.Max()[0] - bounds_ls.Min()[0]);
316
317       
318
319        //////////
320        //-- now apply these values to construct the perspective lispsm matrix
321
322        Matrix4x4 matLispSM;
323       
324        matLispSM = GetFrustum(-1.0, 1.0, -1.0, 1.0, n, n + d);
325
326        // translate to the projection center
327        matLispSM = projectionCenter * matLispSM;
328
329        // transform into OpenGL right handed system
330        Matrix4x4 refl = ScaleMatrix(1.0f, 1.0f, -1.0f);
331        matLispSM *= refl;
332       
333        return matLispSM;
334}
335
336#if 0
337
338Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
339{
340        float maxDist = -1e25f;
341        Vector3 nearest = Vector3::ZERO();
342
343        Matrix4x4 eyeView;
344        mCamera->GetModelViewMatrix(eyeView);
345
346        VertexArray newPts;
347        polyhedron->CollectVertices(newPts);
348       
349        //the LVS volume is always in front of the camera
350        VertexArray::const_iterator it, it_end = pts.end();     
351
352        for (it = pts.begin(); it != it_end; ++ it)
353        {
354                Vector3 pt = *it;
355                Vector3 ptE = eyeView * pt;
356               
357                if (ptE.z > 0) cerr <<"should not happen " << ptE.z << endl;
358                else
359                if (ptE.z > maxDist)
360                {
361                        cout << " d " << ptE.z;
362       
363                        maxDist = ptE.z;
364                        nearest = pt;
365                }
366        }
367
368        //      return Invert(eyeView) * nearest;
369        return nearest;
370}
371
372#else
373
374Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
375{
376        VertexArray newPts;
377        polyhedron->CollectVertices(newPts);
378
379        Vector3 nearest = Vector3::ZERO();
380        float minDist = 1e25f;
381
382        const Vector3 camPos = mCamera->GetPosition();
383
384        VertexArray::const_iterator it, it_end = newPts.end();
385
386        for (it = newPts.begin(); it != it_end; ++ it)
387        {
388                Vector3 pt = *it;
389
390                const float dist = SqrDistance(pt, camPos);
391
392                if (dist < minDist)
393                {
394                        minDist = dist;
395                        nearest = pt;
396                }
397        }
398
399        return nearest;
400}
401
402#endif
403
404Vector3 ShadowMap::GetProjViewDir(const Matrix4x4 &lightSpace,
405                                                                  const VertexArray &pts) const
406{
407        //get the point in the LVS volume that is nearest to the camera
408        const Vector3 e = GetNearCameraPointE(pts);
409
410        //construct edge to transform into light-space
411        const Vector3 b = e + mCamera->GetDirection();
412        //transform to light-space
413        const Vector3 e_lp = lightSpace * e;
414        const Vector3 b_lp = lightSpace * b;
415
416        Vector3 projDir(b_lp - e_lp);
417
418        //project the view direction into the shadow map plane
419        projDir.y = .0f;
420
421        return Normalize(projDir);
422}
423
424
425bool ShadowMap::CalcLightProjection(Matrix4x4 &lightProj)
426{
427        ///////////////////
428        //-- First step: calc frustum clipped by scene box
429
430        DEL_PTR(polyhedron);
431        polyhedron = CalcClippedFrustum(mSceneBox);
432
433        if (!polyhedron) return false; // something is wrong
434
435        // include the part of the light volume that "sees" the frustum
436        // we only require frustum vertices
437
438        VertexArray frustumPoints;
439        IncludeLightVolume(*polyhedron, frustumPoints, mLight->GetDirection(), mSceneBox);
440
441
442        ///////////////
443        //-- transform points from world view to light view and calculate extremal points
444
445        Matrix4x4 lightView;
446        mShadowCam->GetModelViewMatrix(lightView);
447
448        const AxisAlignedBox3 extremalPoints = GetExtremalPoints(lightView, frustumPoints);
449
450        // we use directional lights, so the projection can be set to identity
451        lightProj = IdentityMatrix();
452
453        // switch coordinate system to that used in the lispsm algorithm for calculations
454        Matrix4x4 transform2LispSM = ZeroMatrix();
455
456        transform2LispSM.x[0][0] =  1.0f;
457        transform2LispSM.x[1][2] = -1.0f; // y => -z
458        transform2LispSM.x[2][1] =  1.0f; // z => y
459        transform2LispSM.x[3][3] =  1.0f;
460
461
462        //switch to the lightspace used in the article
463        lightProj *= transform2LispSM;
464
465        const Vector3 projViewDir = GetProjViewDir(lightView * lightProj, frustumPoints);
466
467        //do DirectionalLight Space Perspective shadow mapping
468        //rotate the lightspace so that the projected light view always points upwards
469        //calculate a frame matrix that uses the projViewDir[lightspace] as up vector
470        //look(from position, into the direction of the projected direction, with unchanged up-vector)
471        //const Matrix4x4 frame = MyLookAt2(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());
472        const Matrix4x4 frame = LookAt(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());
473
474        lightProj *= frame;
475
476        const Matrix4x4 matLispSM =
477                CalcLispSMTransform(lightView * lightProj, extremalPoints, frustumPoints);
478
479        lightProj *= matLispSM;
480
481        // change back to GL coordinate system
482        Matrix4x4 transformToGL = ZeroMatrix();
483       
484        transformToGL.x[0][0] =   1.0f;
485        transformToGL.x[1][2] =   1.0f; // z => y
486        transformToGL.x[2][1] =  -1.0f; // y => -z
487        transformToGL.x[3][3] =   1.0f;
488
489        lightProj *= transformToGL;
490
491        AxisAlignedBox3 lightPts = GetExtremalPoints(lightView * lightProj, frustumPoints);
492
493        // focus projection matrix on the extremal points => scale to unit cube
494        Matrix4x4 scaleTranslate = GetFittingProjectionMatrix(lightPts);
495
496        lightProj = lightProj * scaleTranslate;
497
498        Matrix4x4 mymat = lightView * lightProj;
499
500        AxisAlignedBox3 lightPtsNew = GetExtremalPoints(mymat, frustumPoints);
501
502        // we have to flip the signs in order to tranform to opengl right handed system
503        Matrix4x4 refl = ScaleMatrix(1, 1, -1);
504        lightProj *= refl;
505       
506        return true;
507}
508
509
510Polyhedron *ShadowMap::CalcClippedFrustum(const AxisAlignedBox3 &box) const
511{
512        Polyhedron *p = mCamera->ComputeFrustum();
513       
514        Polyhedron *clippedPolyhedron = box.CalcIntersection(*p);
515        DEL_PTR(p);
516       
517        return clippedPolyhedron;
518}
519
520
521//calculates the up vector for the light coordinate frame
522static Vector3 CalcUpVec(const Vector3 viewDir, const Vector3 lightDir)
523{
524        //we do what gluLookAt does...
525        //left is the normalized vector perpendicular to lightDir and viewDir
526        //this means left is the normalvector of the yz-plane from the paper
527        Vector3 left = CrossProd(lightDir, viewDir);
528       
529        //we now can calculate the rotated(in the yz-plane) viewDir vector
530        //and use it as up vector in further transformations
531        Vector3 up = CrossProd(left, lightDir);
532
533        return Normalize(up);
534}
535
536
537void ShadowMap::GetTextureMatrix(Matrix4x4 &m) const
538{
539        m = mTextureMatrix;
540}
541
542 
543unsigned int ShadowMap::GetDepthTexture() const
544{
545        return mFbo->GetDepthTex();
546}
547
548
549unsigned int ShadowMap::GetShadowColorTexture() const
550{
551        return mFbo->GetColorBuffer(0)->GetTexture();
552       
553}
554
555
556void ShadowMap::IncludeLightVolume(const Polyhedron &polyhedron,
557                                                                   VertexArray &frustumPoints,
558                                                                   const Vector3 lightDir,
559                                                                   const AxisAlignedBox3 &sceneBox
560                                                                   )
561{
562        // we don't need closed form anymore => just store vertices
563        VertexArray vertices;
564        polyhedron.CollectVertices(vertices);
565
566        // we 'look' at each point and intect rays with the scene bounding box
567        VertexArray::const_iterator it, it_end = vertices.end();
568
569        for (it = vertices.begin(); it != it_end; ++ it)
570        {
571                Vector3 v  = *it;
572
573                frustumPoints.push_back(v);
574                // hack: start at point which is guaranteed to be outside of box
575                v -= Magnitude(mSceneBox.Diagonal()) * lightDir;
576
577                SimpleRay ray(v, lightDir);
578
579                float tNear, tFar;
580
581                if (sceneBox.Intersects(ray, tNear, tFar))
582                {
583                        Vector3 newpt = ray.Extrap(tNear);
584                        frustumPoints.push_back(newpt);                 
585                }
586        }
587}
588
589
590void ShadowMap::ComputeShadowMap(RenderTraverser *renderer,
591                                                                 const Matrix4x4 &projView)
592{
593        mFbo->Bind();
594       
595        glDrawBuffers(1, mrt);
596
597        glPushAttrib(GL_VIEWPORT_BIT);
598        glViewport(0, 0, mSize, mSize);
599
600        // turn off colors + lighting (should be handled by the render state)
601        glShadeModel(GL_FLAT);
602
603
604        /////////////
605        //-- render scene into shadow map
606
607        _Render(renderer);
608
609        glPopAttrib();
610        glShadeModel(GL_SMOOTH);
611
612#if 0
613        float *data = new float[mSize * mSize];
614
615        GrabDepthBuffer(data, mFbo->GetDepthTex());
616        ExportDepthBuffer(data, mSize);
617
618        delete [] data;
619       
620        PrintGLerror("shadow map");
621#endif
622       
623
624        //////////////
625        //-- compute texture matrix
626
627        static Matrix4x4 biasMatrix(0.5f, 0.0f, 0.0f, 0.5f,
628                                                                0.0f, 0.5f, 0.0f, 0.5f,
629                                                                0.0f, 0.0f, 0.5f, 0.5f,
630                                                                0.0f, 0.0f, 0.0f, 1.0f);
631
632        mTextureMatrix = mLightProjView * biasMatrix;
633
634        FrameBufferObject::Release();
635}
636
637
638void ShadowMap::RenderShadowView(RenderTraverser *renderer,
639                                                                 const Matrix4x4 &projView)
640{
641        glEnable(GL_LIGHTING);
642       
643        _Render(renderer);
644       
645        /*glDisable(GL_LIGHTING);
646        glDisable(GL_DEPTH_TEST);
647
648        Polyhedron *hpoly = CreatePolyhedron(projView, mSceneBox);
649        DrawPoly(hpoly, Vector3(1, 1, 1));
650        DEL_PTR(hpoly);
651
652        glEnable(GL_LIGHTING);
653        glEnable(GL_DEPTH_TEST);*/
654
655        glDisable(GL_POLYGON_OFFSET_FILL);
656}
657
658
659void ShadowMap::_Render(RenderTraverser *renderer)
660{
661        const Vector3 dir = mLight->GetDirection();
662
663        mShadowCam->SetDirection(dir);
664
665        // set position so that we can see the whole scene
666        Vector3 pos = mSceneBox.Center();
667        pos -= dir * Magnitude(mSceneBox.Diagonal() * 0.5f);
668
669        mShadowCam->SetPosition(mCamera->GetPosition());
670
671        const Vector3 upVec = CalcUpVec(mCamera->GetDirection(), dir);
672        const Matrix4x4 lightView = LookAt(mShadowCam->GetPosition(), dir, upVec);
673
674        mShadowCam->mViewOrientation = lightView;
675
676        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
677
678        glPolygonOffset(5.0f, 100.0f);
679        glEnable(GL_POLYGON_OFFSET_FILL);
680       
681        Matrix4x4 lightProj;
682        CalcLightProjection(lightProj);
683
684        mLightProjView = lightView * lightProj;
685
686        DEL_PTR(lightPoly);
687        lightPoly = CreatePolyhedron(mLightProjView, mSceneBox);
688
689        glMatrixMode(GL_PROJECTION);
690        glPushMatrix();
691
692        glMatrixMode(GL_MODELVIEW);
693        glPushMatrix();
694       
695        // set projection matrix manually
696        mShadowCam->mProjection = lightProj;
697       
698        // load gl view projection
699        mShadowCam->SetupViewProjection();
700
701       
702        /////////////
703        //-- render scene into shadow map
704
705        renderer->RenderScene();
706
707
708        glMatrixMode(GL_PROJECTION);
709        glPopMatrix();
710
711        glMatrixMode(GL_MODELVIEW);
712        glPopMatrix();
713
714        glDisable(GL_POLYGON_OFFSET_FILL);
715}
716
717
718} // namespace
Note: See TracBrowser for help on using the repository browser.