root/OpenSceneGraph/trunk/src/osgShadow/LightSpacePerspectiveShadowMap.cpp @ 13041

Revision 13041, 30.2 kB (checked in by robert, 2 years ago)

Ran script to remove trailing spaces and tabs

  • Property svn:eol-style set to native
Line 
1/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2006 Robert Osfield
2 *
3 * This library is open source and may be redistributed and/or modified under
4 * the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
5 * (at your option) any later version.  The full license is in LICENSE file
6 * included with this distribution, and on the openscenegraph.org website.
7 *
8 * This library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11 * OpenSceneGraph Public License for more details.
12 *
13 * ViewDependentShadow codes Copyright (C) 2008 Wojciech Lewandowski
14 * Thanks to to my company http://www.ai.com.pl for allowing me free this work.
15*/
16
17
18#include <osgShadow/LightSpacePerspectiveShadowMap>
19#include <osg/io_utils>
20#include <iostream>
21#include <iomanip>
22
23#define DIRECTIONAL_ONLY     0
24#define DIRECTIONAL_ADAPTED  1
25#define DIRECTIONAL_AND_SPOT 2
26
27//#define LISPSM_ALGO DIRECTIONAL_ONLY
28#define LISPSM_ALGO DIRECTIONAL_ADAPTED
29//#define LISPSM_ALGO DIRECTIONAL_AND_SPOT
30
31#define ROBERTS_TEST_CHANGES 1
32
33#define PRINT_COMPUTED_N_OPT 0
34
35using namespace osgShadow;
36
37////////////////////////////////////////////////////////////////////////////////
38// There are two slightly differing implemetations available on
39// "Light Space Perspective Shadow Maps" page. One from 2004 and other from 2006.
40// Our implementation is written in two versions based on these solutions.
41////////////////////////////////////////////////////////////////////////////////
42// Original LisPSM authors 2004 implementation excerpt. Kept here for reference.
43// DIRECTIONAL AND DIRECTIONAL_ADAPTED versions are based on this code.
44// DIRECTIONAL_AND_SPOT version is based on later 2006 code.
45////////////////////////////////////////////////////////////////////////////////
46 #if 0
47////////////////////////////////////////////////////////////////////////////////
48// This code is copyright Vienna University of Technology, 2004.
49//
50// Please feel FREE to COPY and USE the code to include it in your own work,
51// provided you include this copyright notice.
52// This program is distributed in the hope that it will be useful,
53// but WITHOUT ANY WARRANTY; without even the implied warranty of
54// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
55//
56// Authors Code:
57// Daniel Scherzer (scherzer@cg.tuwien.ac.at)
58//
59// Authors Paper:
60// Michael Wimmer (wimmer@cg.tuwien.ac.at)
61// Daniel Scherzer (scherzer@cg.tuwien.ac.at)
62// Werner Purgathofer
63////////////////////////////////////////////////////////////////////////////////
64void calcLispSMMtx(struct VecPoint* B) {
65    Vector3 min, max;
66    Vector3 up;
67    Matrix4x4 lispMtx;
68    struct VecPoint Bcopy = VECPOINT_NULL;
69    double dotProd = dot(viewDir,lightDir);
70    double sinGamma;
71
72    sinGamma = sqrt(1.0-dotProd*dotProd);
73
74    copyMatrix(lispMtx,IDENTITY);
75
76    copyVecPoint(&Bcopy,*B);
77
78    //CHANGED
79    if(useBodyVec) {
80        Vector3 newDir;
81        calcNewDir(newDir,B);
82        calcUpVec(up,newDir,lightDir);
83    }
84    else {
85        calcUpVec(up,viewDir,lightDir);
86    }
87
88    //temporal light View
89    //look from position(eyePos)
90    //into direction(lightDir)
91    //with up vector(up)
92    look(lightView,eyePos,lightDir,up);
93
94    //transform the light volume points from world into light space
95    transformVecPoint(B,lightView);
96
97    //calculate the cubic hull (an AABB)
98    //of the light space extents of the intersection body B
99    //and save the two extreme points min and max
100    calcCubicHull(min,max,B->points,B->size);
101
102    {
103        //use the formulas of the paper to get n (and f)
104        const double factor = 1.0/sinGamma;
105        const double z_n = factor*nearDist; //often 1
106        const double d = absDouble(max[1]-min[1]); //perspective transform depth //light space y extents
107        const double z_f = z_n + d*sinGamma;
108        const double n = (z_n+sqrt(z_f*z_n))/sinGamma;
109        const double f = n+d;
110        Vector3 pos;
111
112        //new observer point n-1 behind eye position
113        //pos = eyePos-up*(n-nearDist)
114        linCombVector3(pos,eyePos,up,-(n-nearDist));
115
116        look(lightView,pos,lightDir,up);
117
118        //one possibility for a simple perspective transformation matrix
119        //with the two parameters n(near) and f(far) in y direction
120        copyMatrix(lispMtx,IDENTITY);    // a = (f+n)/(f-n); b = -2*f*n/(f-n);
121        lispMtx[ 5] = (f+n)/(f-n);        // [ 1 0 0 0]
122        lispMtx[13] = -2*f*n/(f-n);        // [ 0 a 0 b]
123        lispMtx[ 7] = 1;                // [ 0 0 1 0]
124        lispMtx[15] = 0;                // [ 0 1 0 0]
125
126        //temporal arrangement for the transformation of the points to post-perspective space
127        mult(lightProjection,lispMtx,lightView); // ligthProjection = lispMtx*lightView
128
129        //transform the light volume points from world into the distorted light space
130        transformVecPoint(&Bcopy,lightProjection);
131
132        //calculate the cubic hull (an AABB)
133        //of the light space extents of the intersection body B
134        //and save the two extreme points min and max
135        calcCubicHull(min,max,Bcopy.points,Bcopy.size);
136    }
137
138    //refit to unit cube
139    //this operation calculates a scale translate matrix that
140    //maps the two extreme points min and max into (-1,-1,-1) and (1,1,1)
141    scaleTranslateToFit(lightProjection,min,max);
142
143    //together
144    mult(lightProjection,lightProjection,lispMtx); // ligthProjection = scaleTranslate*lispMtx
145}
146#endif
147
148#if ( LISPSM_ALGO == DIRECTIONAL_ONLY )
149
150
151LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
152{
153    lispsm = NULL;
154}
155
156LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
157{
158}
159
160void LightSpacePerspectiveShadowMapAlgorithm::operator()
161    ( const osgShadow::ConvexPolyhedron* hullShadowedView,
162      const osg::Camera* cameraMain,
163      osg::Camera* cameraShadow ) const
164{
165    osg::BoundingBox bb = hullShadowedView->computeBoundingBox( cameraMain->getViewMatrix() );
166    double nearDist = -bb._max[2];
167
168    const osg::Matrix & eyeViewToWorld = cameraMain->getInverseViewMatrix();
169
170    osg::Matrix lightViewToWorld = cameraShadow->getInverseViewMatrix();
171
172    osg::Vec3d eyePos = osg::Vec3d( 0, 0, 0 ) * eyeViewToWorld;
173
174    osg::Vec3d viewDir( osg::Matrix::transform3x3( osg::Vec3d(0,0,-1), eyeViewToWorld ) );
175
176    osg::Vec3d lightDir( osg::Matrix::transform3x3( osg::Vec3d( 0,0,-1), lightViewToWorld ) );
177    osg::Vec3d up( osg::Matrix::transform3x3( osg::Vec3d(0,1,0), lightViewToWorld ) );
178
179    osg::Matrix lightView; // compute coarse light view matrix
180    lightView.makeLookAt( eyePos, eyePos + lightDir, up );
181    bb = hullShadowedView->computeBoundingBox( lightView );
182
183    const double dotProd = viewDir * lightDir;
184    const double sinGamma = sqrt(1.0- dotProd*dotProd);
185    const double factor = 1.0/sinGamma;
186    const double z_n = factor*nearDist; //often 1
187    //use the formulas of the paper to get n (and f)
188    const double d = fabs( bb._max[1]-bb._min[1]); //perspective transform depth //light space y extents
189    const double z_f = z_n + d*sinGamma;
190    const double n = (z_n+sqrt(z_f*z_n))/sinGamma;
191    const double f = n+d;
192    osg::Vec3d pos = eyePos-up*(n-nearDist);
193
194#if PRINT_COMPUTED_N_OPT
195    std::cout
196       << " N=" << std::setw(8) << n
197       << " n=" << std::setw(8) << z_n
198       << " f=" << std::setw(8) << z_f
199       << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
200       << std::flush;
201#endif
202
203    lightView.makeLookAt( pos, pos + lightDir, up );
204
205    //one possibility for a simple perspective transformation matrix
206    //with the two parameters n(near) and f(far) in y direction
207    double a = (f+n)/(f-n);
208    double b = -2*f*n/(f-n);
209
210    osg::Matrix lispProjection( 1, 0, 0, 0,
211                                0, a, 0, 1,
212                                0, 0,-1, 0,
213                                0, b, 0, 0 );
214
215//    lispProjection.makeIdentity( );
216#if 0
217    {
218        osg::Matrix mvp = _camera->getViewMatrix() *
219                      _camera->getProjectionMatrix();
220
221        extendScenePolytope( mvp, osg::Matrix::inverse( mvp ) );
222    }
223#endif
224
225    bb = hullShadowedView->computeBoundingBox( lightView * lispProjection );
226
227    osg::Matrix fitToUnitFrustum;
228    fitToUnitFrustum.makeOrtho( bb._min[0],  bb._max[0],
229                                bb._min[1],  bb._max[1],
230                               -(bb._min[2]-1), -bb._max[2] );
231
232    cameraShadow->setProjectionMatrix
233        ( lightViewToWorld * lightView * lispProjection * fitToUnitFrustum );
234
235
236#if 0 // DOUBLE CHECK!
237    bb = computeScenePolytopeBounds
238        ( cameraShadow->getViewMatrix() * cameraShadow->getProjectionMatrix() );
239
240    if( !osg::equivalent( 0.f, (bb._min - osg::Vec3d(-1,-1,-1)).length2() ) ||
241        !osg::equivalent( 0.f, (bb._max - osg::Vec3d( 1, 1, 1)).length2() ) )
242    {
243        bb = computeScenePolytopeBounds
244            ( cameraShadow->getViewMatrix() * cameraShadow->getProjectionMatrix() );
245    }
246#endif
247}
248
249#endif
250
251
252
253#if ( LISPSM_ALGO == DIRECTIONAL_ADAPTED )
254
255LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
256{
257    lispsm = NULL;
258}
259
260LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
261{
262}
263
264void LightSpacePerspectiveShadowMapAlgorithm::operator()
265    ( const osgShadow::ConvexPolyhedron* hullShadowedView,
266      const osg::Camera* cameraMain,
267      osg::Camera* cameraShadow ) const
268{
269
270    // all computations are done in post projection light space
271    // which means we are in left handed coordinate system
272    osg::Matrix mvpLight =
273        cameraShadow->getViewMatrix() * cameraShadow->getProjectionMatrix();
274
275    osg::Matrix m = cameraMain->getInverseViewMatrix() * mvpLight;
276    osg::Vec3 eye = osg::Vec3( 0, 0, 0 ) * m;
277    osg::Vec3 center = osg::Vec3( 0, 0, -1 ) * m;
278    osg::Vec3 up(0,1,0);
279    osg::Vec3 viewDir( center - eye );
280    viewDir.normalize();
281
282    m.makeLookAt( eye, center, up );
283
284    osg::BoundingBox bb = hullShadowedView->computeBoundingBox( mvpLight * m );
285    if( !bb.valid() )
286    {
287#if ROBERTS_TEST_CHANGES
288        // OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() invalid bb A"<<std::endl;
289#endif
290        return;
291    }
292
293    double nearDist = -bb._max[2];
294
295#if 1
296    // Original LiSPSM Paper suggests that algorithm should work for all light types:
297    // infinte directional, omnidirectional and spot types may be treated as directional
298    // as all computations are performed in post projection light space.
299    // Frankly, I have my doubts if their error analysis and methodology still works
300    // in non directional lights post projective space. But since I can't prove it doesn't,
301    // I assume it does ;-). So I made an effort to modify their original directional algo
302    // to work in true light post perspective space and compute all params in this space.
303    // And here is a snag. Although shadowed hull fits completely into light space,
304    // camera position may not, and after projective transform it may land outside
305    // light frustum and even on/or below infinity plane. I need camera pos to compute
306    // minimal distance to shadowed hull. If its not right rest of the computation may
307    // be completely off. So in the end this approach is not singulartity free.
308    // I guess this problem is solvable in other way but "this other
309    // way" looks like a topic for other scientific paper and I am definitely not that
310    // ambitious ;-). So for the time being I simply try to discover when this happens and
311    // apply workaround, I found works. This workaround may mean that adjusted projection
312    // may not be optimal in original LisPSM Lmax norm sense. But as I wrote above,
313    // I doubt they are optimal when Light is not directional, anyway.
314
315    // Seems that most nasty case when algorithm fails is when cam pos is
316    // below light frustum near plane but above infinity plane - when this occurs
317    // shadows simply disappear. My workaround is to find this case by
318    // checking light postperspective transform camera z ( negative value means this )
319    // and make sure min distance to shadow hull is clamped to positive value.
320
321    if( eye[2] < 0 && nearDist <= 0 )
322    {
323        float clampedNearDist = 0.0001;
324        eye += viewDir * ( clampedNearDist - nearDist );
325        nearDist = clampedNearDist;
326#if ROBERTS_TEST_CHANGES
327        // OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() adjusting eye"<<std::endl;
328#endif
329
330    }
331#endif
332
333
334#if ROBERTS_TEST_CHANGES
335    if (nearDist<0.0)
336    {
337        nearDist = 0.0;
338        // OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() nearDist<0.0, resetting to 0.0."<<std::endl;
339    }
340#endif
341
342    // Beware!!! Dirty Tricks:
343    // Light direction in light post proj space is actually (0,0,1)
344    // But since we want to pass it to std OpenGL right handed coordinate
345    // makeLookAt function we compensate the effects by also using right
346    // handed view forward vector (ie 0,0,-1) instead.
347    // So in the end we get left handed makeLookAt behaviour (D3D like)...
348    // I agree this method is bizarre. But it works so I left it as is.
349    // It sort of came out by itself through trial and error.
350    // I later understoood why it works.
351
352    osg::Vec3 lightDir(0,0,-1);
353    osg::Matrix lightView; // compute coarse light view matrix
354    lightView.makeLookAt( eye, eye + lightDir, up );
355    bb = hullShadowedView->computeBoundingBox( mvpLight * lightView );
356    if( !bb.valid() )
357    {
358#if ROBERTS_TEST_CHANGES
359        // OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() invalid bb B"<<std::endl;
360#endif
361        return;
362    }
363
364    //use the formulas from the LiSPSM paper to get n (and f)
365    const double dotProd = viewDir * lightDir;
366    const double sinGamma = sqrt(1.0- dotProd*dotProd);
367    const double factor = 1.0/sinGamma;
368    const double z_n = factor*nearDist;
369    //perspective transform depth light space y extents
370    const double d = fabs( bb._max[1]-bb._min[1]);
371    const double z_f = z_n + d*sinGamma;
372    double n = (z_n+sqrt(z_f*z_n))/sinGamma;
373
374#if ROBERTS_TEST_CHANGES
375    // clamp the localtion of p so that it isn't too close to the eye as to cause problems
376    float minRatio=0.02;
377    if (n<d*minRatio)
378    {
379        n=d*minRatio;
380        // OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() n too small, clamping n to "<<minRatio<<"*d."<<std::endl;
381    }
382#endif
383    const double f = n+d;
384    osg::Vec3d pos = eye-up*(n-nearDist);
385    //pos = eye;
386    lightView.makeLookAt( pos, pos + lightDir, up );
387
388    //one possibility for a simple perspective transformation matrix
389    //with the two parameters n(near) and f(far) in y direction
390    double a = (f+n)/(f-n);
391    double b = -2*f*n/(f-n);
392
393    osg::Matrix lispProjection( 1, 0, 0, 0,
394                                0, a, 0, 1,
395                                0, 0, 1, 0,
396                                0, b, 0, 0 );
397
398    cameraShadow->setProjectionMatrix
399        ( cameraShadow->getProjectionMatrix() * lightView * lispProjection );
400
401    //OSG_NOTICE<<"LightSpacePerspectiveShadowMapAlgorithm::operator() normal case dotProd="<<dotProd<<", sinGamma="<<sinGamma<<" d="<<d<<", pos=("<<pos<<"), (n-nearDist)="<<(n-nearDist)<<std::endl;
402    //OSG_NOTICE<<"    eye=("<<eye<<"), up=("<<up<<"), n="<<n<<", nearDist="<<nearDist<<", z_n="<<z_n<<", z_f="<<z_f<<std::endl;
403}
404
405#endif
406
407#if ( LISPSM_ALGO == DIRECTIONAL_AND_SPOT )
408
409
410// Adapted Modified version of LispSM authors implementation from 2006
411// Nopt formula differs from the paper. I adopted original authors class to
412// use with OSG
413
414
415
416//we search the point in the LVS volume that is nearest to the camera
417#include <limits.h>
418static const float OSG_INFINITY = FLT_MAX;
419
420namespace osgShadow {
421
422class LispSM {
423public:
424    typedef std::vector<osg::Vec3d> Vertices;
425
426    void setProjectionMatrix( const osg::Matrix & projectionMatrix )
427        { _projectionMatrix = projectionMatrix; }
428
429    void setViewMatrix( const osg::Matrix & viewMatrix )
430        { _viewMatrix = viewMatrix; }
431
432    void setHull( const ConvexPolyhedron & hull )
433        { _hull = hull; }
434
435    const ConvexPolyhedron & getHull( ) const
436        { return _hull; }
437
438    const osg::Matrix & getProjectionMatrix( void ) const
439        { return _projectionMatrix; }
440
441    const osg::Matrix & getViewMatrix( void ) const
442        { return _viewMatrix; }
443
444    bool getUseLiSPSM() const
445        { return _useLiSPSM; }
446
447    void setUseLiSPSM( bool use )
448        { _useLiSPSM = use; }
449
450    bool getUseFormula() const
451        { return _useFormula; }
452
453    void setUseFormula( bool use )
454        { _useFormula = use; }
455
456    bool getUseOldFormula() const
457        { return _useOldFormula; }
458
459    void setUseOldFormula( bool use )
460        { _useOldFormula = use; }
461
462    void setN(const double& n )
463        { _N = n; }
464
465    const double getN() const
466        { return _N; }
467
468    //for old LispSM formula from paper
469    const double getNearDist() const
470        { return _nearDist; }
471
472    void setNearDist( const double & nearDist )
473        { _nearDist = nearDist; }
474
475    const double getFarDist() const
476        { return _farDist; }
477
478    void setFarDist( const double & farDist )
479        { _farDist = farDist; }
480
481    const osg::Vec3d & getEyeDir() const
482        { return _eyeDir; }
483
484    const osg::Vec3d & getLightDir() const
485        { return _lightDir; }
486
487    void setEyeDir( const osg::Vec3d eyeDir )
488        { _eyeDir = eyeDir; }
489
490    void setLightDir( const osg::Vec3d lightDir )
491        { _lightDir = lightDir; }
492
493protected:
494
495    bool        _useLiSPSM;
496    bool        _useFormula;
497    bool        _useOldFormula;
498    double      _N;
499    double      _nearDist;
500    double      _farDist;
501
502    mutable osg::Vec3d  _E;
503    osg::Vec3d  _eyeDir;
504    osg::Vec3d  _lightDir;
505
506    ConvexPolyhedron _hull;
507
508    osg::Matrix _viewMatrix;
509    osg::Matrix _projectionMatrix;
510
511    double      getN(const osg::Matrix lightSpace, const osg::BoundingBox& B_ls) const;
512
513    osg::Vec3d  getNearCameraPointE() const;
514
515    osg::Vec3d  getZ0_ls
516                    (const osg::Matrix& lightSpace, const osg::Vec3d& e, const double& b_lsZmax, const osg::Vec3d& eyeDir) const;
517
518    double      calcNoptGeneral
519                    (const osg::Matrix lightSpace, const osg::BoundingBox& B_ls) const;
520
521    double      calcNoptOld
522                    ( const double gamma_ = 999) const;
523
524    osg::Matrix getLispSmMtx
525                    (const osg::Matrix& lightSpace) const;
526
527    osg::Vec3d  getProjViewDir_ls
528                    (const osg::Matrix& lightSpace) const;
529
530    void        updateLightMtx
531                    (osg::Matrix& lightView, osg::Matrix& lightProj, const std::vector<osg::Vec3d>& B) const;
532
533public:
534    LispSM( ) : _useLiSPSM( true ), _useFormula( true ), _useOldFormula( false ), _N( 1 ), _nearDist( 1 ), _farDist( 10 ) { }
535
536    virtual void updateLightMtx( osg::Matrix& lightView, osg::Matrix& lightProj ) const;
537};
538
539}
540
541osg::Vec3d LispSM::getNearCameraPointE( ) const
542{
543    const osg::Matrix& eyeView = getViewMatrix();
544
545    ConvexPolyhedron::Vertices LVS;
546    _hull.getPoints( LVS );
547
548    //the LVS volume is always in front of the camera
549    //the camera points along the neg z axis.
550    //-> so the nearest point is the maximum
551
552    unsigned max = 0;
553    for(unsigned i = 0; i < LVS.size(); i++) {
554
555        LVS[i] = LVS[i] * eyeView;
556
557        if( LVS[max].z() < LVS[i].z() ) {
558            max = i;
559        }
560    }
561    //transform back to world space
562    return LVS[max] * osg::Matrix::inverse( eyeView );
563}
564
565//z0 is the point that lies on the plane A parallel to the near plane through e
566//and on the near plane of the C frustum (the plane z = bZmax) and on the line x = e.x
567osg::Vec3d LispSM::getZ0_ls
568    (const osg::Matrix& lightSpace, const osg::Vec3d& e, const double& b_lsZmax, const osg::Vec3d& eyeDir) const
569{
570    //to calculate the parallel plane to the near plane through e we
571    //calculate the plane A with the three points
572    osg::Plane A(eyeDir,e);
573    //to transform plane A into lightSpace
574    A.transform( lightSpace );
575    //transform to light space
576    const osg::Vec3d e_ls = e * lightSpace;
577
578    //z_0 has the x coordinate of e, the z coord of B_lsZmax
579    //and the y coord of the plane A and plane (z==B_lsZmax) intersection
580#if 1
581    osg::Vec3d v = osg::Vec3d(e_ls.x(),0,b_lsZmax);
582
583    // x & z are given. We compute y from equations:
584    // A.distance( x,y,z ) == 0
585    // A.distance( x,y,z ) == A.distance( x,0,z ) + A.normal.y * y
586    // hence A.distance( x,0,z ) == -A.normal.y * y
587
588    v.y() = -A.distance( v ) / A.getNormal().y();
589#else
590    //get the parameters of A from the plane equation n dot d = 0
591    const double d = A.asVec4()[3];
592    const osg::Vec3d n = A.getNormal();
593    osg::Vec3d v(e_ls.x(),(-d-n.z()*b_lsZmax-n.x()*e_ls.x())/n.y(),b_lsZmax);
594#endif
595
596    return v;
597
598}
599
600double LispSM::calcNoptGeneral(const osg::Matrix lightSpace, const osg::BoundingBox& B_ls) const
601{
602    const osg::Matrix& eyeView = getViewMatrix();
603    const osg::Matrix invLightSpace = osg::Matrix::inverse( lightSpace );
604
605    const osg::Vec3d z0_ls = getZ0_ls(lightSpace, _E,B_ls.zMax(),getEyeDir());
606    const osg::Vec3d z1_ls = osg::Vec3d(z0_ls.x(),z0_ls.y(),B_ls.zMin());
607
608    //to world
609    const osg::Vec4d z0_ws = osg::Vec4d( z0_ls, 1 ) * invLightSpace;
610    const osg::Vec4d z1_ws = osg::Vec4d( z1_ls, 1 ) * invLightSpace;
611
612    //to eye
613    const osg::Vec4d z0_cs = z0_ws * eyeView;
614    const osg::Vec4d z1_cs = z1_ws * eyeView;
615
616    double z0 = -z0_cs.z() / z0_cs.w();
617    double z1 = -z1_cs.z() / z1_cs.w();
618
619    if( z1 / z0 <= 1.0 ) {
620
621        // solve camera pos singularity in light space problem brutally:
622        // if extreme points of B projected to Light space extend beyond
623        // camera frustum simply use B extents in camera frustum
624
625        // Its not optimal selection but ceratainly better than negative N
626        osg::BoundingBox bb = _hull.computeBoundingBox( eyeView );
627        z0 = -bb.zMax();
628        if( z0 <= 0 )
629            z0 = 0.1;
630
631        z1 = -bb.zMin();
632        if( z1 <= z0 )
633            z1 = z0 + 0.1;
634    }
635
636    const double d = osg::absolute(B_ls.zMax()-B_ls.zMin());
637
638    double N = d/( sqrt( z1 / z0 ) - 1.0 );
639#if PRINT_COMPUTED_N_OPT
640    std::cout
641       << " N=" << std::setw(8) << N
642       << " n=" << std::setw(8) << z0
643       << " f=" << std::setw(8) << z1
644       << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
645       << std::flush;
646#endif
647    return N;
648}
649
650double LispSM::calcNoptOld( const double gamma_ ) const
651{
652    const double& n = getNearDist();
653    const double& f = getFarDist();
654    const double d = abs(f-n);
655    double sinGamma(0);
656    if(999 == gamma_) {
657        double dot = getEyeDir() * getLightDir();
658        sinGamma = sqrt( 1.0 - dot * dot );
659    }
660    else {
661        sinGamma = sin(gamma_);
662    }
663
664    double N = (n+sqrt(n*(n+d*sinGamma)))/sinGamma;
665#if PRINT_COMPUTED_N_OPT
666    std::cout
667       << " N=" << std::setw(8) << N
668       << " n=" << std::setw(8) << n
669       << " f=" << std::setw(8) << f
670       << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
671       << std::flush;
672#endif
673    return N;
674}
675
676double LispSM::getN(const osg::Matrix lightSpace, const osg::BoundingBox& B_ls) const
677{
678    if( getUseFormula()) {
679        if( getUseOldFormula() )
680            return calcNoptOld();
681        else
682            return calcNoptGeneral(lightSpace,B_ls);
683    }
684    else {
685        return getN();
686    }
687}
688//this is the algorithm discussed in the article
689osg::Matrix LispSM::getLispSmMtx( const osg::Matrix& lightSpace ) const
690{
691    const osg::BoundingBox B_ls = _hull.computeBoundingBox( lightSpace );
692
693    const double n = getN(lightSpace,B_ls);
694
695    //get the coordinates of the near camera point in light space
696    const osg::Vec3d e_ls = _E * lightSpace;
697    //c start has the x and y coordinate of e, the z coord of B.min()
698    const osg::Vec3d Cstart_lp(e_ls.x(),e_ls.y(),B_ls.zMax());
699
700    if( n >= OSG_INFINITY ) {
701        //if n is inf. than we should do uniform shadow mapping
702        return osg::Matrix::identity();
703    }
704    //calc C the projection center
705    //new projection center C, n behind the near plane of P
706    //we work along a negative axis so we transform +n*<the positive axis> == -n*<neg axis>
707    const osg::Vec3d C( Cstart_lp + osg::Vec3d(0,0,1) * n );
708    //construct a translation that moves to the projection center
709    const osg::Matrix projectionCenter = osg::Matrix::translate( -C );
710
711    //calc d the perspective transform depth or light space y extents
712    const double d = osg::absolute(B_ls.zMax()-B_ls.zMin());
713
714    //the lispsm perspective transformation
715
716    //here done with a standard frustum call that maps P onto the unit cube with
717    //corner points [-1,-1,-1] and [1,1,1].
718    //in directX you can use the same mapping and do a mapping to the directX post-perspective cube
719    //with corner points [-1,-1,0] and [1,1,1] as the final step after all the shadow mapping.
720    osg::Matrix P = osg::Matrix::frustum( -1.0,1.0,-1.0,1.0, n, n+d );
721
722    //invert the transform from right handed into left handed coordinate system for the ndc
723    //done by the openGL style frustumGL call
724    //so we stay in a right handed system
725    P = P * osg::Matrix::scale( 1.0,1.0,-1.0 );
726    //return the lispsm frustum with the projection center
727    return projectionCenter * P;
728}
729
730osg::Vec3d LispSM::getProjViewDir_ls(const osg::Matrix& lightSpace ) const {
731    //get the point in the LVS volume that is nearest to the camera
732    const osg::Vec3d e = _E;
733    //construct edge to transform into light-space
734    const osg::Vec3d b = e+getEyeDir();
735    //transform to light-space
736    osg::Vec4d e_lp = osg::Vec4d( e, 1.0 ) * lightSpace;
737    osg::Vec4d b_lp = osg::Vec4d( b, 1.0 ) * lightSpace;
738
739    if( e_lp[3] <= 0 )
740    {
741        e_lp[3] = e_lp[3];
742    }
743
744    if( b_lp[3] <= 0 )
745    {
746        osg::Vec4d v = (e_lp - b_lp)/(e_lp[3]-b_lp[3]);
747
748        v = ( e_lp + v  ) * 0.5;
749
750        b_lp = v;
751    }
752
753    osg::Vec3d projDir( osg::Vec3( b_lp[0], b_lp[1], b_lp[2] ) / b_lp[3] -
754                        osg::Vec3( e_lp[0], e_lp[1], e_lp[2] ) / e_lp[3] );
755
756    projDir.normalize();
757
758    //project the view direction into the shadow map plane
759    projDir.y() = 0.0;
760    return projDir;
761}
762
763void LispSM::updateLightMtx
764    ( osg::Matrix& lightView, osg::Matrix& lightProj ) const
765{
766    //calculate standard light space for spot or directional lights
767    //this routine returns two matrices:
768    //lightview contains the rotated translated frame
769    //lightproj contains in the case of a spot light the spot light perspective transformation
770    //in the case of a directional light a identity matrix
771    // calcLightSpace(lightView,lightProj);
772
773    if( _hull._faces.empty() ) {
774        //debug() << "empty intersection body -> completely inside shadow\n";//debug output
775        return;
776    }
777
778    _E = getNearCameraPointE();
779
780    lightProj = lightProj * osg::Matrix::scale( 1, 1, -1 );
781
782    //coordinate system change for calculations in the article
783    osg::Matrix switchToArticle = osg::Matrix::identity();
784    switchToArticle(1,1) = 0.0;
785    switchToArticle(1,2) =-1.0; // y -> -z
786    switchToArticle(2,1) = 1.0; // z -> y
787    switchToArticle(2,2) = 0.0;
788    //switch to the lightspace used in the article
789    lightProj = lightProj * switchToArticle;
790
791    osg::Matrix L = lightView * lightProj;
792
793    osg::Vec3d projViewDir = getProjViewDir_ls(L);
794
795    if( getUseLiSPSM() /* && projViewDir.z() < 0*/ ) {
796        //do Light Space Perspective shadow mapping
797        //rotate the lightspace so that the proj light view always points upwards
798        //calculate a frame matrix that uses the projViewDir[light-space] as up vector
799        //look(from position, into the direction of the projected direction, with unchanged up-vector)
800        lightProj = lightProj *
801            osg::Matrix::lookAt( osg::Vec3d(0,0,0), projViewDir, osg::Vec3d(0,1,0) );
802
803        osg::Matrix lispsm = getLispSmMtx( lightView * lightProj );
804        lightProj = lightProj * lispsm;
805    }
806
807    const osg::Matrix PL = lightView * lightProj;
808
809    osg::BoundingBox bb = _hull.computeBoundingBox( PL );
810
811    osg::Matrix fitToUnitFrustum;
812    fitToUnitFrustum.makeOrtho( bb._min[0],  bb._max[0],
813                                bb._min[1],  bb._max[1],
814                               -bb._max[2], -bb._min[2] );
815
816    //map to unit cube
817    lightProj = lightProj * fitToUnitFrustum;
818
819    //coordinate system change for calculations in the article
820    osg::Matrix switchToGL = osg::Matrix::identity();
821    switchToGL(1,1) =  0.0;
822    switchToGL(1,2) =  1.0; // y -> z
823    switchToGL(2,1) = -1.0; // z -> -y
824    switchToGL(2,2) =  0.0;
825
826    //back to open gl coordinate system y <-> z
827    lightProj = lightProj * switchToGL;
828    //transform from right handed system into left handed ndc
829    lightProj = lightProj * osg::Matrix::scale(1.0,1.0,-1.0);
830}
831
832void LightSpacePerspectiveShadowMapAlgorithm::operator()
833    ( const osgShadow::ConvexPolyhedron* hullShadowedView,
834      const osg::Camera* cameraMain,
835      osg::Camera* cameraShadow ) const
836{
837    lispsm->setHull( *hullShadowedView );
838    lispsm->setViewMatrix( cameraMain->getViewMatrix() );
839    lispsm->setProjectionMatrix( cameraMain->getViewMatrix() );
840
841#if 1
842    osg::Vec3d lightDir = osg::Matrix::transform3x3( osg::Vec3d( 0, 0, -1 ), osg::Matrix::inverse( cameraShadow->getViewMatrix() ) );
843    osg::Vec3d eyeDir = osg::Matrix::transform3x3( osg::Vec3d( 0, 0, -1 ), osg::Matrix::inverse( cameraMain->getViewMatrix() ) );
844
845#else
846
847    osg::Vec3d lightDir = osg::Matrix::transform3x3( cameraShadow->getViewMatrix(), osg::Vec3d( 0.0, 0.0, -1.0 ) );
848    osg::Vec3d eyeDir = osg::Matrix::transform3x3( cameraMain->getViewMatrix(), osg::Vec3d( 0.0, 0.0, -1.0 ) );
849
850#endif
851
852    lightDir.normalize();
853    eyeDir.normalize();
854
855    lispsm->setLightDir(lightDir);
856
857
858    osg::Matrix &proj = cameraShadow->getProjectionMatrix();
859    double l,r,b,t,n,f;
860    if( proj.getOrtho( l,r,b,t,n,f ) )
861    {
862        osg::Vec3d camPosInLightSpace =
863            osg::Vec3d( 0, 0, 0 ) *
864            osg::Matrix::inverse( cameraMain->getViewMatrix() ) *
865            cameraShadow->getViewMatrix() *
866            cameraShadow->getProjectionMatrix();
867    }
868
869    eyeDir.normalize();
870
871    lispsm->setEyeDir( eyeDir );
872
873    osg::BoundingBox bb =
874        hullShadowedView->computeBoundingBox( cameraMain->getViewMatrix() );
875
876    lispsm->setNearDist( -bb.zMax() );
877    lispsm->setFarDist( -bb.zMin() );
878
879    lispsm->updateLightMtx
880        ( cameraShadow->getViewMatrix(), cameraShadow->getProjectionMatrix() );
881}
882
883LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
884{
885    lispsm = new LispSM;
886}
887
888LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
889{
890    delete lispsm;
891}
892
893
894#endif
Note: See TracBrowser for help on using the browser.