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
35 using 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 ////////////////////////////////////////////////////////////////////////////////
64 void 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
LightSpacePerspectiveShadowMapAlgorithm()151 LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
152 {
153 lispsm = NULL;
154 }
155
~LightSpacePerspectiveShadowMapAlgorithm()156 LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
157 {
158 }
159
operator ()(const osgShadow::ConvexPolyhedron * hullShadowedView,const osg::Camera * cameraMain,osg::Camera * cameraShadow) const160 void 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
LightSpacePerspectiveShadowMapAlgorithm()255 LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
256 {
257 lispsm = NULL;
258 }
259
~LightSpacePerspectiveShadowMapAlgorithm()260 LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
261 {
262 }
263
operator ()(const osgShadow::ConvexPolyhedron * hullShadowedView,const osg::Camera * cameraMain,osg::Camera * cameraShadow) const264 void 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>
418 static const float OSG_INFINITY = FLT_MAX;
419
420 namespace osgShadow {
421
422 class LispSM {
423 public:
424 typedef std::vector<osg::Vec3d> Vertices;
425
setProjectionMatrix(const osg::Matrix & projectionMatrix)426 void setProjectionMatrix( const osg::Matrix & projectionMatrix )
427 { _projectionMatrix = projectionMatrix; }
428
setViewMatrix(const osg::Matrix & viewMatrix)429 void setViewMatrix( const osg::Matrix & viewMatrix )
430 { _viewMatrix = viewMatrix; }
431
setHull(const ConvexPolyhedron & hull)432 void setHull( const ConvexPolyhedron & hull )
433 { _hull = hull; }
434
getHull() const435 const ConvexPolyhedron & getHull( ) const
436 { return _hull; }
437
getProjectionMatrix(void) const438 const osg::Matrix & getProjectionMatrix( void ) const
439 { return _projectionMatrix; }
440
getViewMatrix(void) const441 const osg::Matrix & getViewMatrix( void ) const
442 { return _viewMatrix; }
443
getUseLiSPSM() const444 bool getUseLiSPSM() const
445 { return _useLiSPSM; }
446
setUseLiSPSM(bool use)447 void setUseLiSPSM( bool use )
448 { _useLiSPSM = use; }
449
getUseFormula() const450 bool getUseFormula() const
451 { return _useFormula; }
452
setUseFormula(bool use)453 void setUseFormula( bool use )
454 { _useFormula = use; }
455
getUseOldFormula() const456 bool getUseOldFormula() const
457 { return _useOldFormula; }
458
setUseOldFormula(bool use)459 void setUseOldFormula( bool use )
460 { _useOldFormula = use; }
461
setN(const double & n)462 void setN(const double& n )
463 { _N = n; }
464
getN() const465 const double getN() const
466 { return _N; }
467
468 //for old LispSM formula from paper
getNearDist() const469 const double getNearDist() const
470 { return _nearDist; }
471
setNearDist(const double & nearDist)472 void setNearDist( const double & nearDist )
473 { _nearDist = nearDist; }
474
getFarDist() const475 const double getFarDist() const
476 { return _farDist; }
477
setFarDist(const double & farDist)478 void setFarDist( const double & farDist )
479 { _farDist = farDist; }
480
getEyeDir() const481 const osg::Vec3d & getEyeDir() const
482 { return _eyeDir; }
483
getLightDir() const484 const osg::Vec3d & getLightDir() const
485 { return _lightDir; }
486
setEyeDir(const osg::Vec3d eyeDir)487 void setEyeDir( const osg::Vec3d eyeDir )
488 { _eyeDir = eyeDir; }
489
setLightDir(const osg::Vec3d lightDir)490 void setLightDir( const osg::Vec3d lightDir )
491 { _lightDir = lightDir; }
492
493 protected:
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
533 public:
LispSM()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
getNearCameraPointE() const541 osg::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
getZ0_ls(const osg::Matrix & lightSpace,const osg::Vec3d & e,const double & b_lsZmax,const osg::Vec3d & eyeDir) const567 osg::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
calcNoptGeneral(const osg::Matrix lightSpace,const osg::BoundingBox & B_ls) const600 double 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
calcNoptOld(const double gamma_) const650 double 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
getN(const osg::Matrix lightSpace,const osg::BoundingBox & B_ls) const676 double 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
getLispSmMtx(const osg::Matrix & lightSpace) const689 osg::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
getProjViewDir_ls(const osg::Matrix & lightSpace) const730 osg::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
updateLightMtx(osg::Matrix & lightView,osg::Matrix & lightProj) const763 void 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
operator ()(const osgShadow::ConvexPolyhedron * hullShadowedView,const osg::Camera * cameraMain,osg::Camera * cameraShadow) const832 void 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
LightSpacePerspectiveShadowMapAlgorithm()883 LightSpacePerspectiveShadowMapAlgorithm::LightSpacePerspectiveShadowMapAlgorithm()
884 {
885 lispsm = new LispSM;
886 }
887
~LightSpacePerspectiveShadowMapAlgorithm()888 LightSpacePerspectiveShadowMapAlgorithm::~LightSpacePerspectiveShadowMapAlgorithm()
889 {
890 delete lispsm;
891 }
892
893
894 #endif
895