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

#include "ShadowMapping.h"
#include "FrameBufferObject.h"
#include "RenderState.h"
#include "RenderTraverser.h"
#include "Light.h"
#include "Polygon3.h"
#include "Polyhedron.h"

#include <IL/il.h>
#include <assert.h>


using namespace std;


namespace CHCDemoEngine
{

static CGprogram sCgShadowProgram;
static CGparameter sShadowParam;


static Polyhedron *polyhedron = NULL;
static Polyhedron *lightPoly = NULL;


static void PrintGLerror(char *msg)
{
        GLenum errCode;
        const GLubyte *errStr;
        
        if ((errCode = glGetError()) != GL_NO_ERROR) 
        {
                errStr = gluErrorString(errCode);
                fprintf(stderr,"OpenGL ERROR: %s: %s\n", errStr, msg);
        }
}


static Polyhedron *CreatePolyhedron(const Matrix4x4 &lightMatrix, 
                                                                        const AxisAlignedBox3 &sceneBox)
{
        Frustum frustum(lightMatrix);

        vector<Plane3> clipPlanes;

        for (int i = 0; i < 6; ++ i)
        {
                ////////////
                //-- normalize the coefficients

                // the clipping planes look outward the frustum,
                // so distances > 0 mean that a point is outside
                const float invLength = -1.0f / Magnitude(frustum.mClipPlanes[i].mNormal);

                frustum.mClipPlanes[i].mD *= invLength;
                frustum.mClipPlanes[i].mNormal *= invLength;
        }

        // first create near plane because of precision issues
        clipPlanes.push_back(frustum.mClipPlanes[4]);

        clipPlanes.push_back(frustum.mClipPlanes[0]);
        clipPlanes.push_back(frustum.mClipPlanes[1]);
        clipPlanes.push_back(frustum.mClipPlanes[2]);
        clipPlanes.push_back(frustum.mClipPlanes[3]);
        clipPlanes.push_back(frustum.mClipPlanes[5]);

        return Polyhedron::CreatePolyhedron(clipPlanes, sceneBox);
}


static void GrabDepthBuffer(float *data, GLuint depthTexture)
{
        glEnable(GL_TEXTURE_2D);
        glBindTexture(GL_TEXTURE_2D, depthTexture);

        glGetTexImage(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT, GL_FLOAT, data);

        glBindTexture(GL_TEXTURE_2D, 0);
        glDisable(GL_TEXTURE_2D);
}


static void ExportDepthBuffer(float *data, int size)
{
        ilInit();
        assert(ilGetError() == IL_NO_ERROR);

        ILstring filename = ILstring("shadow.tga");
        ilRegisterType(IL_FLOAT);

        const int depth = 1;
        const int bpp = 1;

        if (!ilTexImage(size, size, depth, bpp, IL_LUMINANCE, IL_FLOAT, data))
        {
                cerr << "IL error " << ilGetError() << endl;
        
                ilShutDown();
                assert(ilGetError() == IL_NO_ERROR);

                return;
        }

        if (!ilSaveImage(filename))
        {
                cerr << "TGA write error " << ilGetError() << endl;
        }

        ilShutDown();
        assert(ilGetError() == IL_NO_ERROR);

        cout << "exported depth buffer" << endl;
}



static AxisAlignedBox3 GetExtremalPoints(const Matrix4x4 &m,
                                                                                 const VertexArray &pts)
{
        AxisAlignedBox3 extremalPoints;
        extremalPoints.Initialize();

        VertexArray::const_iterator it, it_end = pts.end();
                
        for (it = pts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;
                pt = m * pt;

                extremalPoints.Include(pt);
        }

        return extremalPoints;
}


ShadowMap::ShadowMap(DirectionalLight *light
                                         , int size, 
                                         const AxisAlignedBox3 &sceneBox,
                                         Camera *cam):
mSceneBox(sceneBox), mSize(size), mCamera(cam), mLight(light)
{
        mFbo = new FrameBufferObject(size, size, FrameBufferObject::DEPTH_32, true);

        // need a color buffer to keep opengl happy
        mFbo->AddColorBuffer(ColorBufferObject::RGB_UBYTE,
                                 ColorBufferObject::WRAP_CLAMP_TO_EDGE, 
                                                 ColorBufferObject::FILTER_NEAREST);


        mShadowCam = new Camera(mSize, mSize);
        mShadowCam->SetOrtho(true);
}


ShadowMap::~ShadowMap()
{
        DEL_PTR(mFbo);
        DEL_PTR(mShadowCam);
}


static void DrawPolyhedron(Polyhedron *poly, const Vector3 &color)
{
        if (!poly) return;

        for (size_t i = 0; i < poly->NumPolygons(); ++ i)
        {
                glColor3f(color.x, color.y, color.z);

                glBegin(GL_LINE_LOOP);

                Polygon3 *p = poly->GetPolygons()[i];

                for (size_t j = 0; j < p->mVertices.size(); ++ j)
                {
                        Vector3 v = p->mVertices[j];
                        glVertex3d(v.x, v.y, v.z);
                }

                glEnd();
        }
}


void ShadowMap::VisualizeFrustra()
{
        DrawPolyhedron(lightPoly, Vector3(1, 0, 1));
        DrawPolyhedron(polyhedron, Vector3(0, 1, 0));
}


// z0 is the point that lies on the parallel plane to the near plane through e (A)
//and on the near plane of the C frustum (the plane z = bZmax) and on the line x = e.x
Vector3 ShadowMap::GetLightSpaceZ0(const Matrix4x4 &lightSpace, 
                                                                   const Vector3 &e, 
                                                                   const float maxZ, 
                                                                   const Vector3 &eyeDir) const 
{
        // to calculate the parallel plane to the near plane through e we 
        // calculate the plane A with the three points
        Plane3 planeA(e, eyeDir);

        planeA.Transform(lightSpace);
        
        // get the parameters of A from the plane equation n dot d = 0
        const float d = planeA.mD;
        const Vector3 n = planeA.mNormal;
        
        // transform to light space
        const Vector3 e_ls = lightSpace * e;

        Vector3 z0;

        z0.x = e_ls.x;
        z0.y = (d - n.z * maxZ - n.x * e_ls.x) / n.y;
        z0.z = maxZ;

        return z0;
        //return V3(e_ls.x(),(d-n.z()*b_lsZmax-n.x()*e_ls.x())/n.y(),b_lsZmax);
}


float ShadowMap::ComputeNOpt(const Matrix4x4 &lightSpace, 
                                                         const AxisAlignedBox3 &extremalPoints,
                                                         const VertexArray &body) const
{
        const Vector3 nearPt = GetNearCameraPointE(body);
        const Vector3 eyeDir = mCamera->GetDirection();

        Matrix4x4 eyeView;
        mCamera->GetModelViewMatrix(eyeView);

        const Matrix4x4 invLightSpace = Invert(lightSpace);

        const Vector3 z0_ls = GetLightSpaceZ0(lightSpace, nearPt, extremalPoints.Max().z, eyeDir);
        const Vector3 z1_ls = Vector3(z0_ls.x, z0_ls.y, extremalPoints.Min().z);
        
        // transform back to world space
        const Vector3 z0_ws = invLightSpace * z0_ls;
        const Vector3 z1_ws = invLightSpace * z1_ls;

        // transform to eye space
        const Vector3 z0_es = eyeView * z0_ws;
        const Vector3 z1_es = eyeView * z1_ws;

        const float z0 = z0_es.z;
        const float z1 = z1_es.z;

        cout << "z0 ls: " << z0_ls << " z1 ls: " << z1_ls << endl;
        cout << "z0: " << z0_es << " z1: " << z1_es << endl;

        const float d = fabs(extremalPoints.Max()[2] - extremalPoints.Min()[2]); 

        const float n = d / (sqrt(z1 / z0) - 1.0f);

        return n;
}


float ShadowMap::ComputeN(const AxisAlignedBox3 &extremalPoints) const 
{
        const float nearPlane = mCamera->GetNear();
        
        const float d = fabs(extremalPoints.Max()[2] - extremalPoints.Min()[2]); 
        
        const float dotProd = DotProd(mCamera->GetDirection(), mShadowCam->GetDirection());
        const float sinGamma = sin(fabs(acos(dotProd)));

        // test for values close to zero
        if (sinGamma < 1e-6f) return 1e6f;
        
        const float scale = 2.0f;
        return scale * (nearPlane + sqrt(nearPlane * (nearPlane + d * sinGamma))) /  sinGamma;
}


Matrix4x4 ShadowMap::CalcLispSMTransform(const Matrix4x4 &lightSpace, 
                                                                                 const AxisAlignedBox3 &extremalPoints,
                                                                                 const VertexArray &body
                                                                                 )
{
        AxisAlignedBox3 bounds_ls = GetExtremalPoints(lightSpace, body);

        ///////////////
        //-- We apply the lispsm algorithm in order to calculate an optimal light projection matrix
        //-- first find the free parameter values n, and P (the projection center), and the projection depth

        const float n = ComputeN(bounds_ls);
        //const float n = ComputeNOpt(lightSpace, extremalPoints, body); cout << "n: " << n << endl;

        if (n >= 1e6f) // light direction nearly parallel to view => switch to uniform
                return IdentityMatrix();

        const Vector3 nearPt = GetNearCameraPointE(body);
        
        //get the coordinates of the near camera point in light space
        const Vector3 lsNear = lightSpace * nearPt;

        // the start point has the x and y coordinate of e, the z coord of the near plane of the light volume
        const Vector3 startPt = Vector3(lsNear.x, lsNear.y, bounds_ls.Max().z);
        
        // the new projection center
        const Vector3 projCenter = startPt + Vector3::UNIT_Z() * n;

        //construct a translation that moves to the projection center
        const Matrix4x4 projectionCenter = TranslationMatrix(-projCenter);

        // light space y size
        const float d = fabs(bounds_ls.Max()[2] - bounds_ls.Min()[2]);

        const float dy = fabs(bounds_ls.Max()[1] - bounds_ls.Min()[1]);
        const float dx = fabs(bounds_ls.Max()[0] - bounds_ls.Min()[0]);

        

        //////////
        //-- now apply these values to construct the perspective lispsm matrix

        Matrix4x4 matLispSM;
        
        matLispSM = GetFrustum(-1.0, 1.0, -1.0, 1.0, n, n + d);

        // translate to the projection center
        matLispSM = projectionCenter * matLispSM;

        // transform into OpenGL right handed system
        Matrix4x4 refl = ScaleMatrix(1.0f, 1.0f, -1.0f);
        matLispSM *= refl;
        
        return matLispSM;
}

#if 0
Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
{
        float maxDist = -1e25f;
        Vector3 nearest = Vector3::ZERO();

        Matrix4x4 eyeView;
        mCamera->GetModelViewMatrix(eyeView);

        VertexArray newPts;
        polyhedron->CollectVertices(newPts);
        
        //the LVS volume is always in front of the camera
        VertexArray::const_iterator it, it_end = pts.end();     

        for (it = pts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;
                Vector3 ptE = eyeView * pt;
                
                if (ptE.z > 0) cerr <<"should not happen " << ptE.z << endl;
                else
                if (ptE.z > maxDist)
                {
                        cout << " d " << ptE.z;
        
                        maxDist = ptE.z;
                        nearest = pt;
                }
        }

        //      return Invert(eyeView) * nearest;
        return nearest;
}

#else

Vector3 ShadowMap::GetNearCameraPointE(const VertexArray &pts) const
{
        VertexArray newPts;
        polyhedron->CollectVertices(newPts);

        Vector3 nearest = Vector3::ZERO();
        float minDist = 1e25f;

        const Vector3 camPos = mCamera->GetPosition();

        VertexArray::const_iterator it, it_end = newPts.end();

        for (it = newPts.begin(); it != it_end; ++ it)
        {
                Vector3 pt = *it;

                const float dist = SqrDistance(pt, camPos);

                if (dist < minDist)
                {
                        minDist = dist;
                        nearest = pt;
                }
        }

        return nearest;
}

#endif

Vector3 ShadowMap::GetProjViewDir(const Matrix4x4 &lightSpace, 
                                                                  const VertexArray &pts) const 
{
        //get the point in the LVS volume that is nearest to the camera
        const Vector3 e = GetNearCameraPointE(pts);

        //construct edge to transform into light-space
        const Vector3 b = e + mCamera->GetDirection();
        //transform to light-space
        const Vector3 e_lp = lightSpace * e;
        const Vector3 b_lp = lightSpace * b;

        Vector3 projDir(b_lp - e_lp);

        //project the view direction into the shadow map plane
        projDir.y = .0f;

        return Normalize(projDir);
}


bool ShadowMap::CalcLightProjection(Matrix4x4 &lightProj)
{
        ///////////////////
        //-- First step: calc frustum clipped by scene box

        DEL_PTR(polyhedron);
        polyhedron = CalcClippedFrustum(mSceneBox);

        if (!polyhedron) return false; // something is wrong

        // include the part of the light volume that "sees" the frustum
        // we only require frustum vertices

        VertexArray frustumPoints;
        IncludeLightVolume(*polyhedron, frustumPoints, mLight->GetDirection(), mSceneBox);


        ///////////////
        //-- transform points from world view to light view and calculate extremal points

        Matrix4x4 lightView;
        mShadowCam->GetModelViewMatrix(lightView);

        const AxisAlignedBox3 extremalPoints = GetExtremalPoints(lightView, frustumPoints);

        // we use directional lights, so the projection can be set to identity
        lightProj = IdentityMatrix(); 

        // switch coordinate system to that used in the lispsm algorithm for calculations
        Matrix4x4 transform2LispSM = ZeroMatrix();

        transform2LispSM.x[0][0] =  1.0f;
        transform2LispSM.x[1][2] =  -1.0f; // y => -z
        transform2LispSM.x[2][1] =  1.0f; // z => y
        transform2LispSM.x[3][3] =  1.0f;


        //switch to the lightspace used in the article
        lightProj *= transform2LispSM;

        const Vector3 projViewDir = GetProjViewDir(lightView * lightProj, frustumPoints);

        //do DirectionalLight Space Perspective shadow mapping
        //rotate the lightspace so that the projected light view always points upwards
        //calculate a frame matrix that uses the projViewDir[lightspace] as up vector
        //look(from position, into the direction of the projected direction, with unchanged up-vector)
        //const Matrix4x4 frame = MyLookAt2(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());
        const Matrix4x4 frame = LookAt(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());

        lightProj *= frame;

        const Matrix4x4 matLispSM = 
                CalcLispSMTransform(lightView * lightProj, extremalPoints, frustumPoints);

        lightProj *= matLispSM;

        // change back to GL coordinate system
        Matrix4x4 transformToGL = ZeroMatrix();
        
        transformToGL.x[0][0] =   1.0f;
        transformToGL.x[1][2] =   1.0f; // z => y
        transformToGL.x[2][1] =  -1.0f; // y => -z
        transformToGL.x[3][3] =   1.0f;

        lightProj *= transformToGL;

        AxisAlignedBox3 lightPts = GetExtremalPoints(lightView * lightProj, frustumPoints);

        // focus projection matrix on the extremal points => scale to unit cube
        Matrix4x4 scaleTranslate = GetFittingProjectionMatrix(lightPts);

        lightProj = lightProj * scaleTranslate;

        Matrix4x4 mymat = lightView * lightProj;

        AxisAlignedBox3 lightPtsNew = GetExtremalPoints(mymat, frustumPoints);

        // we have to flip the signs in order to tranform to opengl right handed system
        Matrix4x4 refl = ScaleMatrix(1, 1, -1); 
        lightProj *= refl;
        
        return true;
}


Polyhedron *ShadowMap::CalcClippedFrustum(const AxisAlignedBox3 &box) const
{
        Polyhedron *p = mCamera->ComputeFrustum();
        
        Polyhedron *clippedPolyhedron = box.CalcIntersection(*p);

        DEL_PTR(p);
        
        return clippedPolyhedron;
}


//calculates the up vector for the light coordinate frame
static Vector3 CalcUpVec(const Vector3 viewDir, const Vector3 lightDir) 
{
        //we do what gluLookAt does...
        //left is the normalized vector perpendicular to lightDir and viewDir
        //this means left is the normalvector of the yz-plane from the paper
        Vector3 left = CrossProd(lightDir, viewDir);
        
        //we now can calculate the rotated(in the yz-plane) viewDir vector
        //and use it as up vector in further transformations
        Vector3 up = CrossProd(left, lightDir);

        return Normalize(up);
}


void ShadowMap::GetTextureMatrix(Matrix4x4 &m) const
{
        m = mTextureMatrix;
}

 
unsigned int ShadowMap::GetDepthTexture() const
{
        return mFbo->GetDepthTex();
}

unsigned int ShadowMap::GetShadowColorTexture() const
{
        return mFbo->GetColorBuffer(0)->GetTexture();
        
}


void ShadowMap::IncludeLightVolume(const Polyhedron &polyhedron, 
                                                                   VertexArray &frustumPoints, 
                                                                   const Vector3 lightDir,
                                                                   const AxisAlignedBox3 &sceneBox
                                                                   ) 
{
        // we don't need closed form anymore => just store vertices
        VertexArray vertices;
        polyhedron.CollectVertices(vertices);

        // we 'look' at each point and calculate intersections of rays with scene bounding box
        VertexArray::const_iterator it, it_end = vertices.end();

        for (it = vertices.begin(); it != it_end; ++ it)
        {
                Vector3 v  = *it;

                frustumPoints.push_back(v);
                
                // hack: start at point which is guaranteed to be outside of box
                v -= Magnitude(mSceneBox.Diagonal()) * lightDir;

                SimpleRay ray(v, lightDir);

                float tNear, tFar;

                if (sceneBox.Intersects(ray, tNear, tFar))
                {
                        Vector3 newpt = ray.Extrap(tNear);
                        frustumPoints.push_back(newpt);                 
                }
        }
}


void ShadowMap::ComputeShadowMap(RenderTraverser *renderer, const Matrix4x4 &projView)
{
        mFbo->Bind();
        
        glDrawBuffers(1, mrt);

        glPushAttrib(GL_VIEWPORT_BIT);
        glViewport(0, 0, mSize, mSize);

        glDisable(GL_LIGHTING);
        glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);

        glShadeModel(GL_FLAT);


        /////////////
        //-- render scene into shadow map

        _Render(renderer);


        //////////////
        //-- compute texture matrix

        static Matrix4x4 biasMatrix(0.5f, 0.0f, 0.0f, 0.5f,
                                                                0.0f, 0.5f, 0.0f, 0.5f,
                                                                0.0f, 0.0f, 0.5f, 0.5f,
                                                                0.0f, 0.0f, 0.0f, 1.0f);

        mTextureMatrix = mLightProjView * biasMatrix; 

        glPopAttrib();
        
        glShadeModel(GL_SMOOTH);
        glEnable(GL_LIGHTING);
        glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);

#if 0
        float *data = new float[mSize * mSize];

        GrabDepthBuffer(data, mFbo->GetDepthTex());
        ExportDepthBuffer(data, mSize);

        delete [] data;
        
        PrintGLerror("shadow map");
#endif
        
        FrameBufferObject::Release();
}


void ShadowMap::RenderShadowView(RenderTraverser *renderer, const Matrix4x4 &projView)
{
        glEnable(GL_LIGHTING);
        
        _Render(renderer);
        
        /*glDisable(GL_LIGHTING);
        glDisable(GL_DEPTH_TEST);

        //glLineWidth(2);
        Polyhedron *hpoly = CreatePolyhedron(projView, mSceneBox);
        DrawPoly(hpoly, Vector3(1, 1, 1));

        DEL_PTR(hpoly);

        glEnable(GL_LIGHTING);
        glEnable(GL_DEPTH_TEST);*/

        glDisable(GL_POLYGON_OFFSET_FILL);
}


void ShadowMap::_Render(RenderTraverser *renderer)
{
        const Vector3 dir = mLight->GetDirection();

        mShadowCam->SetDirection(dir);

        // set position so that we can see the whole scene
        Vector3 pos = mSceneBox.Center();
        pos -= dir * Magnitude(mSceneBox.Diagonal() * 0.5f);

        mShadowCam->SetPosition(mCamera->GetPosition());

        Vector3 upVec = CalcUpVec(mCamera->GetDirection(), dir);
        Matrix4x4 lightView = LookAt(mShadowCam->GetPosition(), dir, upVec);

        mShadowCam->mViewOrientation = lightView;

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        glPolygonOffset(5.0f, 100.0f);
        glEnable(GL_POLYGON_OFFSET_FILL);
        
        Matrix4x4 lightProj;
        CalcLightProjection(lightProj);

        glMatrixMode(GL_PROJECTION);
        glPushMatrix();
        glLoadMatrixf((float *)lightProj.x);

        mLightProjView = lightView * lightProj;

        DEL_PTR(lightPoly);
        lightPoly = CreatePolyhedron(mLightProjView, mSceneBox);

        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
        
        mShadowCam->SetupCameraView();

        
        /////////////
        //-- render scene into shadow map

        renderer->RenderScene();


        glMatrixMode(GL_PROJECTION);
        glPopMatrix();

        glMatrixMode(GL_MODELVIEW);
        glPopMatrix();

        
        glDisable(GL_POLYGON_OFFSET_FILL);
}
} // namespace