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

#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 Polyhedron *polyhedron = 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 CGprogram sCgShadowProgram;
static CGparameter sShadowParam;


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;
}



ShadowMap::ShadowMap(Light *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);
        // the diffuse color buffer
        mFbo->AddColorBuffer(ColorBufferObject::BUFFER_UBYTE, ColorBufferObject::WRAP_CLAMP_TO_EDGE, ColorBufferObject::FILTER_LINEAR, false);

        mShadowCam = new Camera(mSceneBox.Size().x * 0.5f, mSceneBox.Size().y * 0.5f);
        mShadowCam->SetOrtho(true);
}


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


void ShadowMap::DrawPolys()
{
        if (!polyhedron) return;

        for (size_t i = 0; i < polyhedron->NumPolygons(); ++ i)
        {
                float r = (float)i / polyhedron->NumPolygons();
                float g = 1;
                float b = 1;

                glColor3f(r, g, b);

                glBegin(GL_LINE_LOOP);

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

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

                glEnd();
        }
}


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: get point surely 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);                 
                }
        }
}


bool ShadowMap::CalcLightProjectionLispSM(Matrix4x4 &lightProj)
{
        Matrix4x4 lispMtx = IdentityMatrix();

        const float dotProd = DotProd(mCamera->GetDirection(), mShadowCam->GetDirection());
        const float sinGamma = sqrt(1.0f - dotProd * dotProd);


        ///////////////////
        //-- 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, mShadowCam->GetDirection(), mSceneBox);


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

        AxisAlignedBox3 extremalPoints;
        extremalPoints.Initialize();

        Matrix4x4 lightView;
        mShadowCam->GetModelViewMatrix(lightView);

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

                extremalPoints.Include(pt);
        }


        ///////////////
        //-- We apply the lispsm algorithm in order to calculate an optimal light projection matrix

        // use lispsm formulas
        const float fac = 1.0f / sinGamma;
        
        const float z_n = fac * mCamera->GetNear();

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

        const float z_f = z_n + d * sinGamma;

        // this is the famous n parameter
        const float n = (z_n + sqrt(z_f * z_n)) /sinGamma;
        const float f = n + d;

        const Vector3 nearPt = mShadowCam->GetNear() * mShadowCam->GetDirection() + mShadowCam->GetPosition();

        //get the coordinates of the near camera point in light space
        const Vector3 lsNear = lightView * nearPt;

        //c start has the x and y coordinate of e, the z coord of B.min()
        const Vector3 startPt = Vector3(lsNear.x, lsNear.y, extremalPoints.Max().z);

        // the new projection center
        Vector3 projCenter = startPt - Vector3(0, 0, 1) * n;

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

        lightProj = IdentityMatrix();

        //one possibility for a simple perspective transformation matrix
        //with the two parameters n(near) and f(far) in y direction
        
        // a = (f+n)/(f-n); b = -2*f*n/(f-n);
        // [ 1 0 0 0]
        // [ 0 a 0 b]
        // [ 0 0 1 0]
        // [ 0 1 0 0]

        lightProj.x[1][1] = (f + n) / (f - n);           
        lightProj.x[1][3] = 1.0f;

        lightProj.x[3][1] = -2.0f * f * n / (f - n);            
        lightProj.x[3][3] = 0.0f;                               

        cout << "here4\n" << lightProj << endl;

        // transform into OpenGL right handed system
        Matrix4x4 scale = ScaleMatrix(1.0f, 1.0f, -1.0f);

        lightProj = projectionCenter * scale * lightProj;

        Vector3 pmax = extremalPoints.Max();
        Vector3 pmin = extremalPoints.Min();

        cout << "min: " << lightProj * pmin << endl;
        cout << "max: " << lightProj * pmax << endl;
        

        return true;
}


bool ShadowMap::CalcLightProjection(Matrix4x4 &lightProj)
{
        DEL_PTR(polyhedron);


        ///////////////////
        //-- First step: calc frustum clipped by scene box

        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, mShadowCam->GetDirection(), mSceneBox);


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

        AxisAlignedBox3 extremalPoints;
        extremalPoints.Initialize();

        Matrix4x4 lightView;
        mShadowCam->GetModelViewMatrix(lightView);

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

                extremalPoints.Include(pt);
        }

        // focus projection matrix on the extremal points
        lightProj = GetFittingProjectionMatrix(extremalPoints);

        return true;
}


Polyhedron *ShadowMap::CalcClippedFrustum(const AxisAlignedBox3 &box) const
{
        Vector3 ftl, ftr, fbl, fbr;
        Vector3 ntl, ntr, nbl, nbr;

        VertexArray sides[6];

        mCamera->ComputePoints(ftl, ftr, fbl, fbr, ntl, ntr, nbl, nbr);

        for (int i = 0; i < 6; ++ i)
                sides[i].resize(4);
        
        // left, right
        sides[0][0] = ftl; sides[0][1] = fbl; sides[0][2] = nbl; sides[0][3] = ntl;
        sides[1][0] = fbr; sides[1][1] = ftr; sides[1][2] = ntr; sides[1][3] = nbr;
        // bottom, top
        sides[2][0] = fbl; sides[2][1] = fbr; sides[2][2] = nbr; sides[2][3] = nbl;
        sides[3][0] = ftr; sides[3][1] = ftl; sides[3][2] = ntl; sides[3][3] = ntr;
        // near, far
        sides[4][0] = ntr; sides[4][1] = ntl; sides[4][2] = nbl; sides[4][3] = nbr; 
        sides[5][0] = ftl; sides[5][1] = ftr; sides[5][2] = fbr; sides[5][3] = fbl; 

        //sides[4][0] = ntl; sides[4][1] = ntr; sides[4][2] = nbr; sides[4][3] = nbl; 
        //sides[5][0] = ftr; sides[5][1] = ftl; sides[5][2] = fbl; sides[5][3] = fbr;



        //////////
        //-- compute polyhedron

        PolygonContainer polygons;

        for (int i = 0; i < 6; ++ i)
        {
                Polygon3 *poly = new Polygon3(sides[i]);
                polygons.push_back(poly);
        }

        Polyhedron *p = new Polyhedron(polygons);
        Polyhedron *clippedPolyhedron = box.CalcIntersection(*p);
        
        DEL_PTR(p);
        

        return clippedPolyhedron;
}


void ShadowMap::ComputeShadowMap(RenderTraverser *renderer, const Matrix4x4 &projView)
{
        const float xlen = Magnitude(mSceneBox.Diagonal() * 0.5f);
        const float ylen = Magnitude(mSceneBox.Diagonal() * 0.5f);
        
        const Vector3 dir = mLight->GetDirection();
        //const Vector3 dir(0, 0, 1);

        mShadowCam->SetDirection(dir);
        
        //cout << "lightdir: " << mShadowCam->GetDirection() << endl;

        // set position so that we can see the whole scene
        Vector3 pos = mSceneBox.Center();

        //Matrix4x4 camView;
        //mCamera->GetModelViewMatrix(camView);

        pos -= dir * Magnitude(mSceneBox.Diagonal() * 0.5f);
        mShadowCam->SetPosition(pos);

        mFbo->Bind();
        
        glDrawBuffers(1, mrt);

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

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

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

        glPolygonOffset(1.0f, 2000.0f);
        glEnable(GL_POLYGON_OFFSET_FILL);

        glShadeModel(GL_FLAT);
        glEnable(GL_DEPTH_TEST);

        
        // setup projection
        /*glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
        glOrtho(+xlen, -xlen, +ylen, -ylen, 0.0f, Magnitude(mSceneBox.Diagonal())); 

        Matrix4x4 dummyMat;
        glGetFloatv(GL_PROJECTION_MATRIX, (float *)dummyMat.x);
        cout << "old:\n" << dummyMat << endl;
        */
        Matrix4x4 lightView, lightProj;

        mShadowCam->GetModelViewMatrix(lightView);

        //CalcLightProjection(lightProj);
        CalcLightProjectionLispSM(lightProj);

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

        cout << "new:\n" << lightProj << endl;

        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
        glLoadIdentity();

        mShadowCam->SetupCameraView();

        mLightProjView = lightView * lightProj;


        //////////////
        //-- 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; 




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

        renderer->RenderScene();

        
        glDisable(GL_POLYGON_OFFSET_FILL);
        glMatrixMode(GL_MODELVIEW);
        glPopMatrix();

        glMatrixMode(GL_PROJECTION);
        glPopMatrix();

        glPopAttrib();
#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::GetTextureMatrix(Matrix4x4 &m) const
{
        m = mTextureMatrix;
}

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



} // namespace