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

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



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(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);
        //mFbo->AddColorBuffer(ColorBufferObject::BUFFER_FLOAT_32, ColorBufferObject::WRAP_CLAMP_TO_EDGE, ColorBufferObject::FILTER_LINEAR, false);

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


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


void ShadowMap::DrawPoly(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::DrawPolys()
{
        DrawPoly(lightPoly, Vector3(1, 0, 1));
        DrawPoly(polyhedron, Vector3(0, 1, 0));
}


float ShadowMap::ComputeN(const AxisAlignedBox3 &extremalPoints) const 
{
        const float n = 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)));

        return (n + sqrt(n * (n + d * sinGamma))) /  sinGamma;
}


Matrix4x4 ShadowMap::CalcLispSMTransform(const Matrix4x4 &lightSpace, 
                                                                                 const AxisAlignedBox3 &extremalPoints,
                                                                                 const VertexArray &body
                                                                                 )
{
        //return IdentityMatrix();

        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 = 1e6f;
        const float n = 1e3f;
        //const float n = ComputeN(bounds_ls) * 100;

        cout << "n: " << n << endl;

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

        //c start 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);

        cout << "mx: " <<  bounds_ls.Max() << endl;
        cout << "mn: " << bounds_ls.Min() << endl;

        // the new projection center
        Vector3 projCenter = startPt + Vector3::UNIT_Z() * n;

        cout <<"start: " << startPt << " " << projCenter << " " << Distance(lightSpace * mCamera->GetPosition(), startPt) << endl;

        //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])+dummy;

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

        cout << "d: " << d << " dy: " << dy << " dx: " << dx << endl;

        

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

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

        //cout << "lispsm\n" << matLispSM << endl;

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

        //cout << "new\n" << matLispSM << endl;

        // 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;
//cout<<"i"<< pt.z;
                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;
                }
        }

        //      Invert(eyeView); return 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);

        Matrix4x4 dummy = lightSpace;
        Invert(dummy);
        Vector3 dummyVec = dummy * e_lp;
        Vector3 dummyVec2 = dummy * b_lp;

        //projDir.z = -projDir.z;

        cout << "dummy: " << Normalize(dummyVec2 - dummyVec) << endl;
        //project the view direction into the shadow map plane
        projDir.y = 0.0;

        return Normalize(projDir);
        //return 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, mShadowCam->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 = lightProj * transform2LispSM;

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


        cout << "projViewDir: " << projViewDir << " orig " << mCamera->GetDirection() << endl;

        //do Light 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 = LookAt(Vector3::ZERO(), projViewDir, Vector3::UNIT_Y());

        cout << "frame\n " << frame << endl;
        lightProj = lightProj * frame;

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

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

        Vector3 mydummy = projViewDir;//Vector3::UNIT_Z();
        //cout << "transformed unit z vector: " << Normalize(lightView * lightProj * mydummy) << endl;
        cout << "transformed unit z vector: " << Normalize(frame * mydummy) << endl;

        lightProj = 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 = lightProj * transformToGL;
        //cout << "here4 \n" << lightProj << endl;

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

        //cout << "ma2: " << lightPts.Max() << endl;
        //cout << "mi2: " << lightPts.Min() << endl;

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

        cout << "max: " << lightProj * extremalPoints.Max() << endl;
        cout << "min: " << lightProj * extremalPoints.Min() << endl;

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

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

        //glEnable(GL_LIGHTING);
        //glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
        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);

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

        mLightProjView = lightView * lightProj;

        Frustum frustum(mLightProjView);
        //Frustum frustum(projView);

        DEL_PTR(lightPoly);

        vector<Plane3> clipPlanes;

        for (int i = 0; i < 6; ++ i)
        {
                frustum.mClipPlanes[i].mNormal *= -1;
                frustum.mClipPlanes[i].mD *= -1;

                clipPlanes.push_back(frustum.mClipPlanes[i]);
        }

        DEL_PTR(lightPoly);
        lightPoly = Polyhedron::CreatePolyhedron(clipPlanes, mSceneBox);

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

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

        mShadowCam->SetupCameraView();


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

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


} // namespace